2.1. Formulation

We briefly revisit the gist of the linear stability analysis of wake vortices by Lee & Marcus (Reference Lee and Marcus2023) and then incorporate it into a transient growth analysis. Unless specified otherwise, we use a cylindrical coordinate system $(r, \phi , z)$ . All variables are non-dimensionalised with respect to the characteristic radial length scale, $R_0$ , the characteristic azimuthal velocity scale, $U_0$ , and the fluid density $\rho$ . Detailed definitions of $R_0$ and $U_0$ can be found from Lessen et al. (Reference Lessen, Singh and Paillet1974, p. 755). The base velocity profile $\overline {\boldsymbol{U}}$ , represented as the Batchelor vortex (Batchelor Reference Batchelor1964), or the $q$ -vortex in its non-dimensional form, is given by

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

where $q$ is the swirl parameter that determines the relative strength of the swirling motion. The vortex core region is defined as the radial location where the azimuthal velocity component is maximised, which is $r \leqslant 1.12$ .

The governing equations of fluid motion assume a Newtonian fluid with constant density, $\rho$ , and constant kinematic viscosity, $\nu$ . In terms of the total velocity, $\boldsymbol{u} \equiv u_r \hat {\boldsymbol{e}}_r + u_\phi \hat {\boldsymbol{e}}_\phi + u_z \hat {\boldsymbol{e}}_z$ , and the total specific energy, $\varphi \equiv (\boldsymbol{u} \cdot \boldsymbol{u}) / 2 + p$ , where $p$ denotes the total pressure (non-dimensionalised by $\rho U_0^2$ ), they are expressed as

(2.2) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} = - \boldsymbol{\nabla } \varphi + \boldsymbol{u} \times \boldsymbol{\omega } + \frac {1}{{\textit {Re}}} {\nabla }^2 \boldsymbol{u}\,\,\,\text{with}\,\boldsymbol{\nabla } \cdot \boldsymbol{u} = 0, \end{equation}

where $\boldsymbol{\omega } \equiv \boldsymbol{\nabla } \times \boldsymbol{u}$ is the total vorticity and ${\textit {Re}} \equiv U_{0}R_{0} / \nu$ is the Reynolds number. The equations are linearised around the base flow profile by decomposing $\boldsymbol{u}$ and $p$ into base terms (indicated by overbars $\overline {\ast }$ ) and perturbations (indicated by primes $\ast '$ ). Using the toroidal-poloidal decomposition with $\hat {\boldsymbol{e}}_z$ as a reference vector, the resulting form becomes

(2.3) \begin{equation} \frac {\partial }{\partial t} \begin{pmatrix}\psi ' \\ \chi ' \end{pmatrix} = \mathbb{P} \big ( \boldsymbol{\overline {U}}(r) \times \boldsymbol{\omega }' \big ) - \mathbb{P} \big ( \boldsymbol{\overline {\omega }}(r) \times \boldsymbol{u}' \big ) + \frac {1}{{\textit {Re}}} {\nabla }^2 \begin{pmatrix}\psi ' \\ \chi ' \end{pmatrix}, \end{equation}

where $\psi ' (r,\phi ,z,t)$ and $\chi ' (r,\phi ,z,t)$ are the toroidal and poloidal streamfunctions of $\boldsymbol{u}' (r,\phi ,z,t)$ , respectively. The operator $\mathbb{P}$ decomposes a smooth vector field into its toroidal and poloidal scalar components. If the input vector field is solenoidal, $\mathbb{P}$ is invertible; in other words, if $\boldsymbol{\nabla } \cdot \boldsymbol{u}' = 0$ and $\mathbb{P} ( \boldsymbol{u}' ) = (\psi ',\,\chi ')$ , then $\boldsymbol{u}'$ can be uniquely reconstructed from $(\psi ',\,\chi ')$ (for further details, see Lee & Marcus Reference Lee and Marcus2023, pp. 9–11). Thus, (2.3) only involves two state variables: $\psi '$ and $\chi '$ . This reduced form facilitates the imposition of the analyticity constraint at $r=0$ , as each state variable is treated independently of one another without coupling.

If we introduce the Fourier ansatz (indicated with tildes $\tilde {\ast }$ ) for the azimuthal and axial wavenumbers $m \in \mathbb{Z}$ , $\kappa \in \mathbb{R}\setminus \left \{ 0 \right \}$ , to represent perturbations of finite axial wavelengths, i.e.

(2.4) \begin{equation} \begin{pmatrix}\psi ' \\ \chi ' \end{pmatrix} = \begin{pmatrix} \tilde {\psi }(r,t) \\ \tilde {\chi }(r,t) \end{pmatrix} e^{\mathrm{i}(m\phi + \kappa z)} , \end{equation}

(2.3) further reduces to a spatially one-dimensional form expressed as

(2.5) \begin{equation} \frac {\partial \tilde {\boldsymbol{\wp }}}{\partial t} = \mathcal{L}_{m \kappa }^{\nu } ( \tilde {\boldsymbol{\wp }} ), \end{equation}

where $\tilde {\boldsymbol{\wp }} (r,t) \equiv ( \tilde {\psi } (r,t),\,\tilde {\chi } (r,t) )$ represents the toroidal-poloidal streamfunction set (equivalent to its corresponding velocity Fourier ansatz, $\tilde {\boldsymbol{u}} (r,t) e^{\mathrm{i}(m\phi + \kappa z)}$ ). Here, $\mathcal{L}_{m\kappa }^{\nu }$ is the linear operator with respect to $\tilde {\boldsymbol{\wp }}$ , representing the right-hand side of (2.3) with the inclusion of (2.4). Note that the operator varies with the wavenumbers $m$ and $\kappa$ and the Reynolds number $\textit {Re}$ , as indicated by the subscript and superscript.

To obtain physically meaningful solutions to (2.5) in an unbounded domain ( $0 \leqslant r \lt \infty$ ), the analyticity at the origin and rapid decay as $r \rightarrow \infty$ are necessary. The most prevalent discretisation schemes for computing these solutions have been Chebyshev spectral collocation methods (e.g. Ash & Khorrami Reference Ash, Khorrami and Green1995; Antkowiak & Brancher Reference Antkowiak and Brancher2004; Fontane, Brancher & Fabre Reference Fontane, Brancher and Fabre2008; Mao & Sherwin Reference Mao and Sherwin2011; Muthiah & Samanta Reference Muthiah and Samanta2018), which use a bounded domain requiring two closed ends and therefore demands approximations of the above constraints. In contrast, the mapped Legendre spectral collocation method, described by Lee & Marcus (Reference Lee and Marcus2023), is specifically designed for unbounded domains while accurately satisfying these conditions without additional treatments. When the problem is discretised using either method, we obtain

(2.6) \begin{equation} \frac {{\rm d} \tilde {\boldsymbol{\varrho }}}{{\rm d} t} = \unicode{x1D647}_{m\kappa }^{\nu } \tilde {\boldsymbol{\varrho }} , \end{equation}

where $\tilde {\boldsymbol{\varrho }} (t)$ is the discretised version of $\tilde {\boldsymbol{\wp }} (r,t)$ in spectral space, consisting of the spectral coefficients of $\tilde {\psi }$ and $\tilde {\chi }$ in order, and $\unicode{x1D647}_{m\kappa }^{\nu }$ is the matrix expression of $\mathcal{L}_{m\kappa }^{\nu }$ .

Similarly, we define the discretised version of $\tilde {\boldsymbol{u}} (r,t)$ in physical space as $\tilde {\boldsymbol{\upsilon }} (t)$ , which consists of the collocated values of $\tilde {u}_r$ , $\tilde {u}_\phi$ and $\tilde {u}_z$ in order. The conversion between $\tilde {\boldsymbol{\varrho }}$ and $\tilde {\boldsymbol{\upsilon }}$ can be achieved through the matrix expression of $\mathbb{P}$ , denoted as $\unicode{x1D64B}$ . The construction of $\mathbb{P}$ based on the mapped Legendre spectral collocation method is described by Matsushima & Marcus (Reference Matsushima and Marcus1997, pp. 330–333). We use the following notation: $\tilde {\boldsymbol{\varrho }} = \unicode{x1D64B}^{\dagger } \tilde {\boldsymbol{\upsilon }}$ and $\tilde {\boldsymbol{\upsilon }} = \unicode{x1D64B} \tilde {\boldsymbol{\varrho }}$ . Under the solenoidal velocity assumption, both $\unicode{x1D64B}^{\dagger }\unicode{x1D64B}$ and $\unicode{x1D64B}\unicode{x1D64B}^{\dagger }$ can be treated as identity maps. Now we define the ‘energy’ $E$ of a velocity of the form $\tilde {\boldsymbol{u}} (r,t) e^{\mathrm{i}(m\phi + \kappa z)}$ as

(2.7) \begin{equation} E(\tilde {\boldsymbol{u}}) \equiv \int _{0}^{\infty } \left ( \tilde {u}_r^* \tilde {u}_r + \tilde {u}_\phi ^* \tilde {u}_\phi + \tilde {u}_z^* \tilde {u}_z \right ) r{\rm d}r. \end{equation}

A similar usage can be found from Mao & Sherwin (Reference Mao and Sherwin2012). Using numerical integration, e.g. a quadrature rule, $E(\tilde {\boldsymbol{u}})$ is expressed as

(2.8) \begin{equation} E(\tilde {\boldsymbol{u}}) = \tilde {\boldsymbol{\upsilon }}^* \unicode{x1D648} \tilde {\boldsymbol{\upsilon }} = \tilde {\boldsymbol{\varrho }}^{*} \unicode{x1D64B}^{\kern1pt *} \unicode{x1D648} \unicode{x1D64B} \tilde {\boldsymbol{\varrho }} , \end{equation}

where $\unicode{x1D648}$ represents the numerical integration form of (2.7). A specific example of $\unicode{x1D648}$ using the Gauss–Legendre quadrature rule is given in Appendix A.

At last, we apply the transient growth formalism (Schmid & Henningson Reference Schmid and Henningson2001; Mao & Sherwin Reference Mao and Sherwin2012; Muthiah & Samanta Reference Muthiah and Samanta2018) to complete our formulation. Consider a set of eigenmodes of $\unicode{x1D647}_{m\kappa }^{\nu }$ containing $p$ elements $\left \{ \tilde {\boldsymbol{\varrho }}_1,\,\tilde {\boldsymbol{\varrho }}_2, \ldots ,\, \tilde {\boldsymbol{\varrho }}_p \right \}$ , corresponding to eigenvalues $\left \{ \sigma _1,\,\sigma _2,\ldots ,\,\sigma _p \right \} \subset \mathbb{C}$ , respectively. Assuming that $\tilde {\boldsymbol{\varrho }}$ belongs to the eigenspace spanned by these $p$ eigenmodes, i.e.

(2.9) \begin{equation} \tilde {\boldsymbol{\varrho }} (t) = \sum _{k=1}^{p} \tilde {\xi }_k e^{\sigma _k t} \tilde {\boldsymbol{\varrho }}_k, \end{equation}

we use a new vector $\tilde {\boldsymbol{\xi }}_0 \equiv (\tilde {\xi }_1, \tilde {\xi }_2, \ldots , \tilde {\xi }_p ) \in \mathbb{C}^p$ to represent $\tilde {\boldsymbol{\varrho }}$ . For instance, at time $t=\tau$ , $\tilde {\boldsymbol{\varrho }} (\tau )$ is expressed as $\exp (\tau \unicode{x1D64E}) \tilde {\boldsymbol{\xi }}_0$ , where $\unicode{x1D64E} \equiv \text{diag}(\sigma _1, \sigma _2, \ldots , \sigma _p)$ . Focusing on the transient growth process, the chosen eigenmodes are assumed to be asymptotically stable in time ( ${\textrm {Re}} (\sigma _k) \lt 0$ ), which mostly holds for strong swirling $q$ -vortices. By defining , (2.9) at time $t = \tau$ becomes

(2.10) \begin{equation} \tilde {\boldsymbol{\varrho }} (\tau ) = \unicode{x1D651} {\exp (\tau \unicode{x1D64E})} \tilde {\boldsymbol{\xi }}_0, \end{equation}

and applying this to (2.8) leads to the energy formula using the following $\ell ^2$ norm:

(2.11) \begin{equation} \begin{aligned} E(\tilde {\boldsymbol{u}}(\tau )) & = \tilde {\boldsymbol{\xi }}_0^* {\exp (\tau \unicode{x1D64E}^*)} \unicode{x1D651}^* \unicode{x1D64B}^* \unicode{x1D648} \unicode{x1D64B} \unicode{x1D651} {\exp (\tau \unicode{x1D64E})} \tilde {\boldsymbol{\xi }}_0 \\ & = \big \lVert \unicode{x1D641} \exp (\tau \unicode{x1D64E}) \tilde {\boldsymbol{\xi }}_0 \big \rVert _2^2, \end{aligned} \end{equation}

where the matrix $\unicode{x1D641}$ is defined such that $\unicode{x1D641}^* \unicode{x1D641} = \unicode{x1D651}^* \unicode{x1D64B}^* \unicode{x1D648} \unicode{x1D64B} \unicode{x1D651}$ . The maximum energy growth, $G$ , which determines the optimal perturbations under the transient growth formalism, at time $t=\tau$ is given by

(2.12) \begin{equation} \begin{aligned} G(\tau ) & \equiv \sup _{\tilde {\boldsymbol{u}}(0) \neq \boldsymbol{0}} \frac {E(\tilde {\boldsymbol{u}}(\tau ))}{E(\tilde {\boldsymbol{u}}(0))} = \sup _{\tilde {\boldsymbol{\xi }}_0 \neq \boldsymbol{0}} \frac { \big \lVert \unicode{x1D641} \exp (\tau \unicode{x1D64E}) \tilde {\boldsymbol{\xi }}_0 \big \rVert _2^2}{ \big \lVert \unicode{x1D641} \tilde {\boldsymbol{\xi }}_0 \big \rVert _2^2} \\ & = \big \lVert \unicode{x1D641} \exp (\tau \unicode{x1D64E}) \unicode{x1D641}^{\kern1pt -1} \big \rVert _2^2. & \end{aligned} \end{equation}

Using the fact that the $L^2$ norm of an arbitrary matrix is the same as its largest singular value, we finally reach the following. For the largest singular value $\varsigma _{1}$ (assumed to be non-zero) of $\unicode{x1D641} \exp (\tau \unicode{x1D64E}) \unicode{x1D641}^{\kern1pt -1}$ and its associated right and left singular vectors $\boldsymbol{r}_{1}$ and $\boldsymbol{l}_{1}$ , i.e.

(2.13) \begin{equation} \unicode{x1D641} \exp (\tau \unicode{x1D64E}\kern1pt) \unicode{x1D641}^{\kern1pt -1} \boldsymbol{r}_{1} = \varsigma _{1}\boldsymbol{l}_{1}, \end{equation}

which can result from the singular value decomposition (SVD), we get

(2.14) \begin{equation} G(\tau ) = \varsigma _{1}^2, \end{equation}

and the optimal perturbation velocity input and output at time $t = \tau$ are

(2.15) \begin{equation} \tilde {\boldsymbol{\upsilon }}_{\textit{opt}}(0) = \unicode{x1D64B} \unicode{x1D651} \unicode{x1D641}^{\kern1pt -1} \boldsymbol{r}_{1} \,\,\,\text{and}\,\,\, \tilde {\boldsymbol{\upsilon }}_{\textit{opt}}(\tau ) = \varsigma _{1}\unicode{x1D64B} \unicode{x1D651} \unicode{x1D641}^{\kern1pt -1} \boldsymbol{l}_{1} . \end{equation}

2.2. Numerical parameters

Discretisation is essential for the current problem formulation, and care must be taken to minimise non-physical errors arising from numerical parameters. In this study, we employ the mapped Legendre spectral collocation method, which is well-suited for analysing rotating flows in unbounded domains. The key numerical parameters in this scheme include the number of spectral basis elements $M$ , the number of radial collocation points $N$ and the map parameter $L$ – defined in (A2). Detailed procedures for accurately resolving the eigenmodes are discussed by Lee & Marcus (Reference Lee and Marcus2023), which interested readers are encouraged to read. In our set-up, we keep $M + |m| = N$ to ensure $N \geqslant M$ , and we vary $L$ to explore both potential and viscous critical-layer eigenmode families by adjusting the scheme’s characteristic numerical resolution.

For comparison, we secondarily consider the Chebyshev spectral collocation method with domain truncation at $r = R_\infty$ . The domain of the Chebyshev polynomials is linearly mapped from $[-1, 1]$ to $[0, R_\infty ]$ , as favoured in previous studies (e.g. Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989; Mao & Sherwin Reference Mao and Sherwin2011), where the primitive variables $\tilde {\boldsymbol{u}}$ and $\tilde {p}$ are taken as state variables. However, the use of this scheme in the present study is solely restricted to investigating how sensitively the domain truncation affects the transient growth analysis outcome, in comparison to the mapped Legendre spectral collocation method developed for radially unbounded domains (Matsushima & Marcus Reference Matsushima and Marcus1997; Lee & Marcus Reference Lee and Marcus2023), serving as a significant drawback despite its constructional convenience for computation.

2.3. Physical parameters

There are five physical parameters that influence the nature of the problem: $\tau$ , $q$ , $\textit {Re}$ , $m$ and $\kappa$ . Below, we clarify the range or value of each parameter that will be the focus of this study.

The maximum energy growth $G$ is explicitly dependent on the total time of growth $\tau$ , indicating the duration over which we allow linear transient dynamics of the wake vortex to develop. In the context of aircraft trailing vortices, an upper limit on $\tau$ is identifiable due to the dominance of the Crow instability mechanism after several hundred time units ( ${=}R_0/U_0$ ). For example, Matsushima & Marcus (Reference Matsushima and Marcus1997, pp. 341–343) reported the prevalence of long-wavelength instability around $t = 200$ in simulations of a counter-rotating vortex pair configuration. Under proper rescaling of units, the trailing vortex simulation by Han et al. (Reference Han, Lin, Schowalter, Arya and Proctor2000, pp. 295–297, also see figure 10) exhibited vortex linkage at $t = 229.12$ under moderate ambient turbulence. Based on these findings, we concentrate on the region where $\tau \lesssim O(10^2)$ . Our typical attention to the transient growth is in the time range of $10 \lt \tau \leqslant 100$ , where relatively fast transient growth is expected (Mao & Sherwin Reference Mao and Sherwin2012). However, we note that a longer range may be explored in case it is needed to verify long-term characteristics of transient growth.

Additionally, the analysis outcomes are also subject to physical parameters such as the swirl parameter $q$ , the Reynolds number $\textit {Re}$ , and the azimuthal and axial wavenumbers $m$ and $\kappa$ . For this study, we fix the first two parameters and, unless specified otherwise, use $q = 4$ and ${\textit {Re}} = 10^5$ . These values represent conditions where the swirling motion is sufficiently strong to exclude significant linear instabilities (e.g. $q \geqslant 2.31$ , see Heaton Reference Heaton2007a ), and where viscous diffusion is small enough to treat the base vortex profile as quasi-steady. It is remarked that, according to experiment-based estimation by Fabre and Jacquin (Reference Fabre and Jacquin2004, p. 259), this set-up aligns with the condition of actual trailing vortices behind large transport aircraft. As for the perturbation wavenumbers, we turn attention to axisymmetric or helical cases ( $m=0\,\text{or}\,1$ ) with small axial wavenumbers of order unity or less. This choice is driven not only by their prevalence in vortex transient growth literature (e.g. Antkowiak & Brancher Reference Antkowiak and Brancher2004; Pradeep & Hussain Reference Pradeep and Hussain2006; Mao & Sherwin Reference Mao and Sherwin2012; Navrose et al. Reference Navrose, Jonhson, Brion, Jacquin and Robinet2018), but also by the anticipation that such low-frequency perturbations better account for the principal perturbation structure that we later aim to initiate via particles near the vortex.

3.1. Numerical sensitivity and proper discretisation

We construct optimal perturbations by combining the eigenmodes of the wake vortex. To obtain accurate results, the chosen computation scheme should reliably capture each physically relevant eigenmode family in a well-resolved manner, while maintaining insensitivity to variations in numerical parameters. We evaluate the numerical sensitivity of the mapped Legendre spectral collocation method and the Chebyshev spectral collocation method to determine which is more appropriate for the present analysis.

When considering viscous eigenmodes that are regular across the entire radial domain, including their asymptotic behaviours near the origin and as $r$ approaches infinity, there are three important eigenmode families: the discrete family, the potential family and the viscous critical-layer family (Lee & Marcus Reference Lee and Marcus2023). The first family, as its name suggests, is associated with discrete spectra (i.e. sets of eigenvalues) and each eigenmode’s spatial structure is uniquely characterised by the number of ‘wiggles’ clustered in or around the vortex core. The other two families comprise continuous eigenmodes, whose spatial structures vary continuously and are associated with continuous spectra.

Figure 1 shows the numerically resolved spectra of the $q$ -vortex using the following physical parameters: $(m,\,\kappa ,\,q,\,{\textit {Re}}) = (1,\,1.0,\,4.0,\,10^5)$ , envisioning the families of eigenmodes. Both the Chebyshev spectral collocation method (with $M=800$ ) and the mapped Legendre spectral collocation method (with $M=400$ ) were employed for this analysis. To illustrate the continuous spectra, we collected all numerical eigenvalues obtained by varying the domain truncation radius $R_\infty$ , ranging from $12$ to $13$ for the Chebyshev spectral method, and by adjusting the map parameter $L$ , spanning from $3$ to $3.1$ for the mapped Legendre spectral method. These variations in $R_\infty$ or $L$ introduce small shifts in the numerically resolved continuous eigenvalues, making it possible to trace continuous spectra (see Lee & Marcus Reference Lee and Marcus2023, § 6.4.2).

Given perfect resolution, the free stream and potential spectra are expected to stretch out to ${\textrm {Re}}(\sigma ) \rightarrow -\infty$ . In figure 1, the spectra are shown in two panels with different aspect ratios. The left panel extends to large $|{\textrm {Re}}(\sigma )|$ , showcasing both the free stream and spurious spectra. The eigenmodes related to the free stream spectrum and the spurious spectrum exhibit non-regular characteristics: the free stream eigenmodes are singular since they do not decay to zero as $r \rightarrow \infty$ (Mao & Sherwin Reference Mao and Sherwin2011), and the spurious eigenmodes are non-physical, characterised by irregular oscillations near the origin (Lee & Marcus Reference Lee and Marcus2023). The right panel displays the discrete, potential and viscous critical-layer spectra, which correspond to the regular eigenmode families discussed earlier.

Figure 1. Numerical spectra of the $q$ -vortex with $(m,\,\kappa ,\,q,\,{\textit {Re}}) = (1,\,1.0,\,4.0,\,10^5)$ using the Chebyshev spectral collocation method (grey squares) and the mapped Legendre spectral collocation method (black dots). The consistent discrete spectrum demonstrates the robustness of both methods. The free stream and spurious spectra (shown in the left panel) are associated with either singular or non-physical eigenmodes, making them irrelevant to this study, with most generated by the Chebyshev method. The discrete, potential and viscous critical-layer spectra (shown in the right panel) are associated with regular eigenmodes, with the mapped Legendre method providing a clearer distinction of the viscous critical-layer spectrum.

A few issues arise when using the Chebyshev spectral method instead of the mapped Legendre spectral method for resolving the eigenmodes. As depicted in the left panel of figure 1, a significant portion of the numerically resolved spectra accounts for eigenmode families that are either singular or non-physical, making them irrelevant to the present problem. This issue likely arises from approximating the asymptotic constraints through subordinate boundary conditions at both ends of the computational domain. For the Chebyshev spectral method, as for $m = 1$ , the boundary conditions implemented are

(3.1) \begin{equation} \frac {{\rm d} \tilde {u}_r}{{\rm d}r} = \frac {{\rm d} \tilde {u}_\phi }{{\rm d}r} = \tilde {u}_z = \tilde {p} = 0\,\,\,\,\text{at}\,\,r=0 \end{equation}

and

(3.2) \begin{equation} \tilde {u}_r = \tilde {u}_\phi = \tilde {u}_z = \tilde {p} = 0\,\,\,\,\text{at}\,\,r=R_\infty , \end{equation}

which are proxies for analyticity at the origin and rapid decay as $r \rightarrow \infty$ , respectively (see Ash & Khorrami Reference Ash, Khorrami and Green1995). While (3.1) and (3.2) may serve as necessary conditions for what they are supposed to mimic, they cannot be considered formally equivalent. For instance, (3.2) does not prohibit solutions from oscillating in the far field as long as the oscillation is momentarily zeroed out at $r=R_\infty$ , which explains the emergence of the free stream spectrum.

The second issue comes from the unclear distinction between the viscous critical-layer spectrum and the potential spectrum. As noted by Lee & Marcus (Reference Lee and Marcus2023, pp. 41–42), the Chebyshev spectral method produces scattered traces of the viscous critical-layer spectrum curves, making it challenging to distinguish these curves from the surrounding continuous region. This scattering can be attributed to the high sensitivity of continuous spectra to minor errors. In the Chebyshev spectral method, domain truncation removes spatial information far from the origin, which, albeit diminutive, holds physical significance. Figure 2 illustrates the difference between the viscous critical-layer eigenmodes and the potential ones. Despite their structural resemblance on a large scale, the viscous critical-layer eigenmodes retain the structure of their inviscid counterparts beyond the region where viscosity effects dominate locally, scaled of the order of $Re^{-1/3}$ (Lin Reference Lin1955). In contrast, the potential eigenmodes turns into null outside this region, epitomising their ‘wave packet’ form (Mao & Sherwin Reference Mao and Sherwin2011), which conforms to the twist condition presented by Trefethen & Embree (Reference Trefethen and Embree2005, pp. 98–114). Further details of their comparison are omitted in this article; they are elucidated by Lee & Marcus (Reference Lee and Marcus2023).

A numerical sensitivity test evaluating the maximum energy growth $G$ at $\tau = 10$ for the entire eigenspace, using both methods, is presented in figure 3. The remaining physical parameters are kept the same: $(m,\,\kappa ,\,q,\,{\textit {Re}}) = (1,\,1.0,\,4.0,\,10^5)$ . As expected, increasing the number of spectral elements $M$ reduces sensitivity to changes in numerical parameters for both methods. At a fixed $M$ , the map parameter $L$ acts as a resolution-tuning parameter in the mapped Legendre spectral collocation method, while the domain truncation radius $R_\infty$ serves this role in the Chebyshev spectral collocation method. The parameter test ranges shown in figure 3 are based on the typical usage found by Mao & Sherwin (Reference Mao and Sherwin2011) and Lee & Marcus (Reference Lee and Marcus2023). Changes in $L$ have minimal impact on $G(\tau = 10)$ within the mapped Legendre spectral collocation method’s test range ( $1 \leqslant L \leqslant 10$ ), therefore allowing for arbitrary selection of $L$ within this interval. In contrast, $G(\tau = 10)$ is notably influenced by variations in $R_\infty$ within the Chebyshev spectral collocation method’s test range ( $10 \leqslant R_\infty \leqslant 20$ ), especially as $R_\infty$ increases. This presents a challenge, as using a large $R_\infty$ should be preferred to preserve the unbounded nature of the radial domain. We found that manually excluding the sub-eigenspace spanned by free stream eigenmodes mitigates this issue, which is, in fact, a step that is proactively taken in the mapped Legendre spectral collocation method. Exclusion of these eigenmodes is reasonable from a physical standpoint, as their non-decaying behaviour implies that they analytically possess infinite energy. This, in turn, renders them formally inapplicable in the current transient growth analysis context.

Figure 2. Schematic comparison between viscous critical-layer and potential eigenmodes, illustrating a velocity component. Both exhibit a similar large-scale structure within the region where viscosity effects are locally dominant, corresponding to a singularity at $r=r_c$ in the inviscid limit (for details of this inviscid singularity, see Lee & Marcus Reference Lee and Marcus2023). The width of this large-scale structure scales as $O(Re^{-1/3})$ (see Lin Reference Lin1955) for the viscous critical-layer eigenmode. In contrast, the width of the potential eigenmode can vary even at the same $Re$ , forming a ‘wave packet’ in compliance with the pseudomode analysis (see Trefethen & Embree Reference Trefethen and Embree2005). Aside from the location $r=r_c$ , the viscous critical-layer eigenmode retains the structure of its inviscid counterpart, while the potential eigenmode simply turns into null.

Figure 3. Numerical sensitivity test evaluating $G(\tau = 10)$ for the entire eigenspace associated with $(m,\,\kappa ,\,q,\,{\textit {Re}}) = (1,\,1.0,\,4.0,\,10^5)$ using $(a)$ the mapped Legendre spectral collocation method with varying $L$ and $(b)$ the Chebyshev spectral collocation method with varying $R_\infty$ , for $M=400,\,600$ and $800$ . The test ranges reflect typical parameter usage. $L$ can be arbitrarily adjusted without impacting domain unboundedness. In contrast, setting a finite $R_\infty$ compromises the unbounded nature of the problem, and therefore, larger values of $R_\infty$ should be preferred. However, increasing $R_\infty$ leads to undesirable numerical sensitivity.

To recapitulate, when it comes to resolving the eigenmodes of the $q$ -vortex in a radially unbounded domain, the Chebyshev spectral collocation method with domain truncation faces several challenges, which unfavourably influence numerical sensitivity in transient growth evaluation. In contrast, the mapped Legendre spectral collocation method effectively mitigates these numerical limitations, making it a more suitable choice for the present problem. Therefore, we adopt the mapped Legendre spectral collocation method for examining the transient growth of the wake vortex.

Figure 4. Maximum energy growth as a function of total growth time: $(a)$ for $m=0$ (axisymmetric) and $(b)$ for $m=1$ (helical). The values of $G$ are evaluated from the entire eigenspace, incorporating the discrete, potential and viscous critical-layer families as basis elements. Here, $q=4$ and ${\textit {Re}} = 10^5$ .

3.2. Maximum energy growth

Mao & Sherwin (Reference Mao and Sherwin2012) demonstrated that transient growth primarily results from the non-normality of continuous eigenmodes, while discrete eigenmodes play a less significant role. In their analysis, they used the term ‘continuous eigenmodes’ as a compilation of potential and free stream eigenmodes. However, as we pointed out earlier, free stream eigenmodes are unsuitable for evaluating maximum energy growth because their energy reaches infinity. Thus, the term ‘continuous eigenmodes’ in their argument should more specifically refer to potential eigenmodes. Not only that, but their argument also requires further refinement, as it did not account for viscous critical-layer eigenmodes. This eigenmode family was not distinguished from the potential family, presumably due to spectral overlap between the two (see Mao & Sherwin Reference Mao and Sherwin2011, p. 8) and their large-scale structural similarity, as shown in figure 2.

Accordingly, we believe that the argument put forth by Mao & Sherwin (Reference Mao and Sherwin2012) still necessitates clarification about which continuous eigenmode family predominantly contributes to optimal perturbations that maximise energy growth: the potential family or the viscous critical-layer family. To that end, we first evaluate $G (\tau )$ across the entire eigenspace and then compare the results with those from different sub-eigenspaces, each spanned by a distinct eigenmode family.

Figure 4 presents the numerically evaluated values of $G(\tau )$ from the entire eigenspace at various wavenumbers $\kappa$ for the $m=0$ and $m=1$ cases. In the $m=0$ cases, as shown in figure 4 $(a)$ , where perturbations are two-dimensional (i.e. functions of $r$ and $z$ only), the dependence of the $G$ curves on $\kappa$ is clear; in the short run, growth is stronger with larger $\kappa$ , while in the long run, the largest $G$ is achieved with smaller $\kappa$ . One may check in figure 4 $(a)$ that the upper envelope of the curves sequentially corresponds to decreasing $\kappa$ as $\tau$ increases. This trend aligns with previous observations in the literature (Pradeep & Hussain Reference Pradeep and Hussain2006; Mao & Sherwin Reference Mao and Sherwin2012), supporting the validity of our evaluation.

In the $m=1$ cases, depicted in figure 4 $(b)$ , involving three-dimensional perturbations, the $\kappa$ -dependence of the $G$ curves becomes complex, as previously noted by Antkowiak & Brancher (Reference Antkowiak and Brancher2004) and Pradeep & Hussain (Reference Pradeep and Hussain2006). To further clarify this trend, figure 5(a) presents supplementary slices of $G$ as a function of axial wavenumber $\kappa$ for $m=1$ at four growth periods within the time range of interest: $\tau =31.6$ , $50$ , $75$ and $100$ . As $\tau$ increases, a local energy growth peak becomes more pronounced, especially at $\tau =100$ , where the local peak around $\kappa = 2.0$ nearly stands at the largest $G$ at $\kappa = 0.1$ . This peak feature in the $G$ – $\kappa$ curves aligns with the findings of Antkowiak & Brancher (Reference Antkowiak and Brancher2004, p. L3), who ascribed the intricate (stretching and tilting) nature of three-dimensional perturbations to such irregularities.

In figure 5(b), the local maximum of $G$ around $\tau = 100$ , denoted $G_{{max}}$ , is plotted alongside its growth time, $\tau _{{max}}$ , as a function of $\kappa$ , provided this maximum is identifiable (e.g. see the $\kappa = 1.0$ or $2.0$ cases in figure 4 b). This approach follows Antkowiak & Brancher (Reference Antkowiak and Brancher2004) and Pradeep & Hussain (Reference Pradeep and Hussain2006), who used $G_{{max}}$ to examine local energy growth features. Similar to the $G{-}\kappa$ curve at $\tau =100$ , the local energy growth peak is observed , as $\tau _{{max}}$ consistently appears near $\tau = 100$ . Note that we avoid using the ‘global’ maximum of $G$ over the entire range of $\tau$ , as suggested in the literature. Due to the higher order of magnitude of ${\textit {Re}}(=10^5)$ in the present study compared with that in the literature ( ${\lesssim}10^4$ ), the global maximum of $G$ occurs at a much larger growth time, $\tau \sim O(10^3{-}10^4)$ . This time range far exceeds the intervals considered both in the literature and our study ( $\tau \lesssim 10^2$ ), thus placing it outside the current scope of analysis.

Figure 5. (a) Maximum energy growth as a function of axial wavenumber for $m=1$ , at specified growth periods of $\tau = 31.6,\,50,\,75$ and $100$ ; and (b) local maximum of $G$ across growth time around $\tau = 100$ and corresponding growth time. As in figure 4, the values of $G$ are evaluated from the entire eigenspace. Here, $q=4$ and ${\textit {Re}} = 10^5$ .

When comparing the $m=1$ cases with the $m=0$ cases, focusing on relatively short-term growth, we find that the largest $G$ for $m=1$ generally exceeds that for $m=0$ . For example, at $\tau = 10^2$ , the largest $G$ among the evaluated values is $7.2 \times 10^2$ for $\kappa = 0.1$ at $m=1$ , whereas it is $2.8 \times 10^2$ for $\kappa = 5.0$ at $m=0$ .

The maximum energy growth curves, evaluated from the entire eigenspace and compared with those from sub-eigenspaces respectively spanned by the discrete family, the viscous critical-layer family and the potential family, are presented in figure 6. It is evident that the curves derived from the entire eigenspace are primarily reproduced by those obtained from the sub-eigenspace of the viscous critical-layer family, highlighting its dominant contribution. However, the values of $G(\tau )$ from the rest of the continuous sub-eigenspace, for which the potential eigenmodes account, are of a similar magnitude to those from the discrete sub-eigenspace. Thus, the contribution of the potential family to transient growth is as minor as that of the discrete family.

Figure 6. Comparison of the curves of $G(\tau )$ evaluated from different sub-eigenspaces, each spanned by a distinct eigenmode family: $(a)$ $(m,\,\kappa ) = (0,\,0.1)$ ; $(b)$ $(m,\,\kappa ) = (0,\,5.0)$ ; $(c)$ $(m,\,\kappa ) = (1,\,0.1)$ and $(d)$ $(m,\,\kappa ) = (1,\,5.0)$ . Here, $q=4$ and ${\textit {Re}} = 10^5$ . The maximum energy growth curves from the entire eigenspace, identical to those plotted in figure 4, are nearly reproduced by those from the sub-eigenspace spanned by the viscous critical-layer family.

There are two minor exceptions worth noting. For $m=0$ , the discrete family contributes as significantly as the viscous critical-layer family to short-term optimal growth, particularly when $\kappa$ is small. We believe that this is relevant to the fact that, in the limit of $\kappa \rightarrow 0$ , critical layers vanish and so do the derived continuous eigenmodes, while the discrete ones persist. Additionally, for $m=1$ , the potential family’s contribution to optimal growth slightly supersedes that of the viscous critical-layer family during a brief period ( $\tau \lt 10$ ), accounting for the presence of a quirk in the $G$ curves around $\tau = 5.6$ . However, the exception clears quickly beyond this period, and the maximum $G$ attainable during this period never exceeds $10$ , thus not overturning the general dominance of the viscous critical-layer family in transient growth.

Based on these observations, we revisit the demonstration provided by Mao & Sherwin (Reference Mao and Sherwin2012) with the following clarification; the non-normality of the continuous eigenmodes induces significant transient growth, and it is specifically the viscous critical-layer family that predominantly contributes to this growth, rather than the potential family. The distinction is important, as it addresses the ‘true’ origin of transient growth of the wake vortex as critical layers. The potential eigenmodes have their theoretical root in the wave packet pseudomode analysis (Trefethen & Embree Reference Trefethen and Embree2005). As showcased in figure 2, they omit the asymptotic information of critical layers, making their birth irrelevant to phenomena that require asymptotic matching or equivalently, critical layer analysis (Lin Reference Lin1955; Le Dizès Reference Le Dizès2004). Although the wave packet pseudomode analysis is a powerful tool for exploring all possible forms of continually varying eigensolutions, it can divert attention too much from the genuine gems more worthy of our focus. Furthermore, this clarification better aligns with the argument made by Heaton (Reference Heaton2007 b), who suggested that inviscid continuous spectrum (CS) transients dominate growth over short time intervals. The viscous critical-layer family in our classfication corresponds to the viscous regularisation of the inviscid CS, which we denoted as the inviscid critical-layer spectrum (Lee & Marcus Reference Lee and Marcus2023), when ${\textit {Re}} \lt \infty$ .

3.3. Perturbation structures

The effects of changing $m$ and $\kappa$ on the perturbation structures that lead to optimal transient growth have been widely investigated and are well-established in the context of linear vortex dynamics (Antkowiak & Brancher Reference Antkowiak and Brancher2004; Pradeep & Hussain Reference Pradeep and Hussain2006; Mao & Sherwin Reference Mao and Sherwin2012). In this section, we examine whether our transient growth calculation complies with these established findings and then conduct a comparative analysis of the perturbation structures across different values of $m$ and $\kappa$ for further consideration.

Pradeep & Hussain (Reference Pradeep and Hussain2006) reported that axisymmetric perturbations ( $m=0$ ) generally produce the largest energy growth, as illustrated in figure 4. However, as the largest $G$ increases, the total duration of perturbation growth also lengthens (i.e. the $\tau$ needed to achieve $G$ increases), as the spatial structure shifts further from the vortex core, necessitating longer time for interactions to occur. For helical perturbations ( $m=1$ ), a common spatial structure emerges, with the main motion concentrated around a specific radius near the vortex core. In these cases, the growth can potentially trigger fluctuations within the vortex core, even though the initial perturbation originates outside it. This mechanism has occasionally been identified as a cause of erratic long-wavelength displacements in experimental vortices (e.g. Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta2016; Bölle et al. Reference Bölle, Brion, Couliou and Molton2023), often termed ‘vortex meandering’ (see Antkowiak & Brancher Reference Antkowiak and Brancher2004, p. L4). Mao & Sherwin (Reference Mao and Sherwin2012) affirmed that the vortex meandering phenomenon can be driven by the transient response of the vortex to an out-of-core perturbation.

Figure 7. Optimal perturbation inputs with unit energy $E$ (as defined in (2.7)) and their amplified outputs at $t=\tau$ , shown through absolute velocity components alongside their corresponding three-dimensional structures. Each structure is represented by the isosurface at 50 % of the maximum specific energy in space. Dark and light colours indicate counterclockwise and clockwise swirling of the flow, respectively. Here, four representative cases with the largest value of $G$ from figure 4 are displayed: $(a,\,b)$ the axisymmetric cases $(m,\,\kappa ) = (0,\,5.0)$ for $\tau = 31.6$ and $\tau = 100$ , respectively; and $(c,\,d)$ the helical cases $(m,\,\kappa ) = (1,\,1.0)$ for $\tau = 31.6$ and $\tau = 100$ , respectively. The initial dominance of azimuthal velocity components is found in all cases, while for $m=1$ , energy is transferred into the core region $r \leqslant 1.12$ (shaded in each plot).

As mentioned in § 2.3, our focus is on the relatively short time period of $10 \lt \tau \leqslant 100$ to study the transient growth process. In longer periods, classical linear instability mechanisms like the Crow instability may dominate under real conditions. Within this time range, the largest values of $G$ attained from our considerations (see figure 4) occur at $\kappa = 5.0$ for the $m=0$ cases and at $\kappa = 0.1$ for the $m=1$ cases. We consider these cases as representative. In figure 7, we depict the optimal perturbation velocity inputs and outputs for $\tau = 31.6$ and $\tau =100$ , all of which are visualised by the absolute velocity components. For clearer visualisation, we portray their corresponding three-dimensional structures alongside, represented by the iso-surface at 50 % of the maximum specific energy in physical space, i.e. $|(\tilde {\boldsymbol{u}}(r) e^{\mathrm{i}(m \phi + \kappa z)} + \mathrm{c.c.})|^2/2$ , where $\mathrm{c.c.}$ stands for the complex conjugate of the antecedent term. Dark and light surfaces express counterclockwise and clockwise swirling directions, respectively.

In all cases, the following characteristics are consistently observed. First, azimuthal velocity components are initially dominant in all optimal perturbations, while the other velocity components evolve significantly towards the end of the growth period. This clearly indicates that the azimuthal velocity component should be prioritised when inducing these optimal perturbations from an unperturbed state. Second, the most energetic part of the optimal perturbation inputs, coinciding with the peak of the absolute azimuthal velocity component, tends to be distant from the vortex core as $\tau$ increases. This tendency is found to be more evident in panels $(c)$ and $(d)$ , where the major perturbation structure overlaps the core region at $\tau = 31.6$ but moves out of the core at $\tau = 100$ .

For panels $(a)$ and $(b)$ , where $(m,\,\kappa )=(0,\,5.0)$ , the input perturbations generally form a ring structure owing to their azimuthal symmetry. As the perturbation evolves, the radius of the ring remain largely unchanged. This tendency for local confinement of the optimal perturbation structures becomes more pronounced with increasing $\kappa$ . This indicates that perturbations with shorter axial wavelengths have a more localised influence around the initially perturbed region.

For panels $(c)$ and $(d)$ , where $(m,\,\kappa )=(1,\,0.1)$ , a spiral structure develops in the most energetic region of the input perturbation due to alternating layers of oppositely swirling fluid motions at the periphery of the vortex core. Unlike the axisymmetric cases, the perturbation structure undergoes a drastic transformation from its input to output states. Notably, the most energetic region of the perturbation, initially located outside the vortex core, eventually penetrates into the vortex core. During this process, the transverse velocity ( $\tilde {u}_r$ and $\tilde {u}_\phi$ ) becomes maximal at the vortex centre. In other words, the principal response of the vortex to optimal perturbations with $m=1$ is characterised by the transverse motion of the vortex core, which is likely linked to the vortex meandering phenomenon (Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta2016; Bölle Reference Bölle2021).

It is important to clarify that the induction of vortex meandering by optimal helical perturbations with an axially long wavelength was formerly given by Mao & Sherwin (Reference Mao and Sherwin2012). They employed mesh-based direct numerical simulations rather than the matrix-based analysis (corresponding to (2.6)–(2.15) in our formulation), even though they used a matrix-based approach when $m=0$ . In a way, this choice seems to have been made to address challenges related to analyticity at the origin, which depends on the value of $m$ (see Lee & Marcus Reference Lee and Marcus2023, pp. 51–52). In contrast, our approach, using the mapped Legendre spectral collocation method, is fundamentally designed to be robust for any value of $m$ . Therefore, our contribution here lies in confirming the same phenomenon linked to $m=1$ perturbations using a computationally fast and formally consistent matrix-based transient growth analysis.

3.4. Nonlinear impacts on an optimally perturbed vortex

Given an optimally perturbed vortex, the linearised theory (see § 2.1) predicts that the perturbation will gradually amplify as time approaches $t=\tau$ , after which it decays in the absence of extrinsic factors capable of triggering secondary instabilities from the most perturbed state. This section is dedicated to verifying whether such transient behaviour remains significant in the original nonlinear system, governed by (2.2), despite the influence of higher-order energy transfer across different wavenumbers and other nonlinear effects, which could potentially cause early vortex growth to deviate from the linear prediction.

Although the optimal perturbation structures vary with the selection of $m$ , $\kappa$ and $\tau$ , our primary aim here is to investigate their general trend of evolution over time, as anticipated by the linearised theory. We focus on a specific case of optimal perturbation where $(m\,,\kappa ) = (1\,,0.1)$ and the optimal growth time is $\tau =50$ . The $z$ -component of the perturbation vorticity input, $\omega '_z (t=0)$ , on the $z=0$ plane is illustrated in figure 8. This case was chosen because the substantial shift of the most energetic portion of the perturbation from the periphery to the vortex core, as shown in figure 7 $(c,\,d)$ , offers a clear illustration of vortex growth. However, we emphasise that this particular behaviour at $m=1$ (potentially related to vortex meandering) is not the primary focus of this study. For readers interested in vortex meandering, we suggest referring to Edstrand et al. (Reference Edstrand, Davis, Schmid, Taira and Cattafesta2016) and Bölle (Reference Bölle2021).

According to the linearised theory, the optimal perturbation velocity input, expressed as $\tilde {\boldsymbol{\upsilon }}_{\textit{opt}}(0) = \unicode{x1D64B} \unicode{x1D651} \unicode{x1D641}^{\kern1pt -1} \boldsymbol{r}_{1}$ as in (2.15), evolves at $t=\eta$ as

(3.3) \begin{equation} \tilde {\boldsymbol{\upsilon }}_{\textit{opt}}(\eta ) = \unicode{x1D64B} \unicode{x1D651} \exp (\eta \unicode{x1D64E}) \unicode{x1D641}^{\kern1pt -1} \boldsymbol{r}_{1}. \end{equation}

One may check the consistency of the above equation when $\eta = \tau$ , with $\tilde {\boldsymbol{\upsilon }}_{\textit{opt}}(\tau )$ given in (2.15). We label this prediction from the linearised theory as ‘linear’. The ‘linear’ prediction, however, may be ideal as it strips off all higher-order interactions coming from the nonlinear convection term, i.e. $\boldsymbol{u} \times \boldsymbol{\omega }$ in (2.2), which facilitates energy transfer from the perturbation wavenumbers $(m,\,\kappa )$ to their multiples (e.g. $(2m,\,2\kappa )$ , $(3m,\,3\kappa ),\,\ldots$ ) or vice versa. We denote the growth of the optimal perturbation, accounting for higher-order interactions, as ‘nonlinear’. The extent of this nonlinearity substantially depends on the initial perturbation’s energy level. If the perturbation energy approaches zero (or the perturbation is infinitesimal), the ‘nonlinear’ evolution should follow the ‘linear’ prediction. We set aside the numerical details of our nonlinear simulations in Appendix B.

Figure 8. Axial perturbation vorticity contour on the $z=0$ plane of the optimal input for the case where $(m,\,\kappa ) = (1,\,0.1)$ with the optimal growth time $\tau =50$ . Solid contour lines represent positive levels (67 % and 33 % of the absolute maximum) and dashed lines indicate negative levels (–33 % and –67 % of the absolute maximum).

In the nonlinear simulations, the initial velocity field $\boldsymbol{u}(t=0)$ is defined as

(3.4) \begin{equation} \boldsymbol{u}(r,\phi ,z,t=0\,;\,\varepsilon ) \equiv \overline {\boldsymbol{U}}(r) + \varepsilon \tilde {\boldsymbol{u}}_{{opt}}(r) e^{i(m \phi + \kappa z)} + \mathrm{c.c.}, \end{equation}

where $\varepsilon$ determines how intense the initial perturbation is, adjusting the perturbation energy input. The base term representing the unperturbed $q$ -vortex, $\overline {U}(r)$ , was assumed to be unchanging in time in the linear analysis, as its radial viscous diffusion is negligible due to the high $Re$ in this problem set-up. In contrast, the nonlinear simulations take this small viscous diffusion of the base $q$ -vortex into account for enhanced accuracy. That is to say, even the unperturbed flow changes slowly over time, which can be calculated with $\varepsilon = 0$ . As a result, the perturbation velocity field at $t=\eta$ is assessed as the difference between two time-varying fields, i.e.

(3.5) \begin{equation} \boldsymbol{u}'(r,\phi ,z,\eta \,;\,\varepsilon ) \equiv \boldsymbol{u}(r,\phi ,z,\eta \,;\,\varepsilon ) - \boldsymbol{u}(r,\eta \,;\,0). \end{equation}

The perturbation energy at $t=\eta$ , denoted $E(\eta )$ , is evaluated as the volume integration of $\boldsymbol{u}' \cdot \boldsymbol{u}'$ divided by $2\pi$ times the axial wavelength ( $= 2\pi \cdot 2\pi / \kappa$ ), for consistency with the energy definition in (2.7).

In figure 9, three energy growth curves are plotted together for comparison. First, the energy growth curve in the linear evolution case peaks at $t=\tau =50$ with a maximum energy growth of $G = 3 \times 10^2$ . This curve serves as an index of the linear process’ prevalence during the early transient growth of vortices. Next, the energy growth curve in the nonlinear evolution case with $E(0) = 10^{-8}$ aligns with the ‘linear index’ curve. In this scenario, nonlinear effects arise but remain minimal, showing a slight debilitation of maximum energy growth at $t=\tau =50$ . However, it is unlikely that this change is entirely due to the nonlinearity introduced to the system because, in comparison with Mao & Sherwin (Reference Mao and Sherwin2012, p. 55), such a drop in energy growth at the peak may also be attributed to viscous diffusion of the base flow over time. Lastly, the energy growth curve in the nonlinear evolution case with $E(0) = 10^{-3}$ exhibits more pronounced debilitation at the peak. Compared with the previous case, this relatively high energy case demonstrates more intensification of the nonlinearity, represented by a secondary energy hump around $t=80$ that is unpredicted by the linear case. Nevertheless, the general trend of the curve does not drift away from the linear prediction. Overall, the coherence in trend holds particularly well up to the maximum vortex growth period ( $t \lt \tau = 50$ ), which suggests that the linearised transient growth framework remains effective in describing the early evolution of the original nonlinear system. Beyond this period, the nonlinearity intensifies, but the trend of decay remains persistent.

Figure 9. Comparison of transient energy growth curves for the ‘linear’ evolution case, calculated via (3.3), and for the ‘nonlinear’ evolution cases with initial perturbation energies of $10^{-8}$ and of $10^{-3}$ , obtained from three-dimensional nonlinear simulations. The initial perturbation is depicted in figure 8.

The prevalence of the linear process in the early-stage vortex growth becomes more evident when examining the evolution of the perturbation structure, as shown in figure 10. Using the same contour style across the three cases discussed above, we present three snapshots of axial vorticity perturbation contours on the $z=0$ plane at $t=25$ , $t=50$ and $t=100$ for each case. The structural coherence in vorticity perturbation between the linear and nonlinear cases is apparent at $t=25$ , representing the stage of rapid perturbation growth. At $t=50$ , the optimal growth time, the perturbation structures remain largely coherent. However, in the nonlinear evolution case with $E(0)=10^{-3}$ , a weak breakdown from the $m=1$ symmetry becomes observable, indicating that other azimuthal wavenumbers rather than $m=1$ begin to gain non-negligible energy through higher-order energy transfer across different wavenumbers. By $t=100$ , corresponding to the stage of asymptotic stabilisation, the perturbation structures no longer show strong resemblance. This is another indication of the nonlinearity intensification resulting from prolonged vortex growth over time.

Figure 10. Axial vorticity perturbation contours of the optimally perturbed vortex (refer to figure 8 for the initial perturbation) on the $z=0$ plane at $t=25$ , $t=50$ and $t=100$ : $(a)$ the ‘linear’ evolution case where maximum energy growth is known to occur at $t=\tau =50$ ; $(b)$ the ‘nonlinear’ evolution case with initial perturbation energy of $10^{-8}$ ; $(c)$ the ‘nonlinear’ evolution case with initial perturbation energy of $10^{-3}$ . The same contour style as figure 8 is used for all plots. Despite the presence and intensification of nonlinearity with increasing perturbation energy over time, the ‘linear’ process largely prevails the overall dynamics during early vortex growth in the original nonlinear system.

Last but not least, we note that the simulation with initial perturbation energy of $10^{-3}$ results in substantial displacement of the vortex core, as shown in figure 11, notwithstanding the seemingly small initial energy level. The $\lambda _2$ -isosurface at $\lambda _2 = -0.05$ is used to detect the vortex core (see Jeong & Hussain Reference Jeong and Hussain1995). The maximum displacement of the vortex centre in the simulation is nearly equal to the core radius, coinciding with experimentally observed meandering amplitudes of similar scale (see Devenport, Zsoldos & Vogel Reference Devenport, Zsoldos and Vogel1997; Bölle Reference Bölle2021). Based on the rough figures for a large transport aircraft from Fabre & Jacquin (Reference Fabre and Jacquin2004, p. 259), the characteristic scales in our formulation are $U_0 \approx 27\,\mathrm{m\, s^{-1}}$ and $R_0 \approx 0.5\,\mathrm{m}$ . Using the density of air $\rho \approx 1\,\mathrm {kg\, m^{-3}}$ , the ‘dimensionless’ energy of $10^{-3}$ translates to an ‘actual’ kinetic energy of $(10^{-3}) \times 2\pi \rho (U_0^2/2) R_0^2 \approx 0.6\,\mathrm{J\, m^{-1}}$ (‘per metre’ stands for axial unit length), which appears to be not exorbitant in practice. We believe this strengthens the practicability of the optimal transient growth process under consideration, along with the radially concentrated nature of the optimal perturbation structures (see figure 8).

Figure 11. Three-dimensional illustration of the $q$ -vortex with a helical perturbation (see figure 8), initiated with an energy level of $10^{-3}$ : $(a)$ the initial core structure at $t=0$ , and $(b)$ the most excited core structure at $t=\tau =50$ . To detect the vortex core, the $\lambda _2$ -isosurface at $\lambda _2 = -0.05$ is depicted (see Jeong & Hussain Reference Jeong and Hussain1995). The maximum displacement of the vortex centre comes up to the order of the core radius, comparable with experimental observations of vortex meandering amplitude (see Devenport et al. Reference Devenport, Zsoldos and Vogel1997; Bölle Reference Bölle2021).

4.1. On a means of initiating optimal transient growth

The transient evolution of the optimally perturbed $q$ -vortex, as analysed in the previous section, unveils a promising way for significantly disturbing the vortex, even when the base vortex is known to be linearly stable. For this growth process to hold practical significance, an important question remains: by what means can such perturbations be initiated? In our analyses and simulations so far, the presence of perturbations has been presumed to be initially present along with the vortex supposedly in an undisturbed state. However, from a practical standpoint, there should be $q$ -vortex, which by itself is quasi-steady. Without a plausible initiation process, the optimal transient vortex growth process may remain purely theoretical.

Ambient turbulence might seem a compelling means of initiation, as often addressed in linear instability contexts (e.g. Crow & Bate Reference Crow and Bate1976; Han et al. Reference Han, Lin, Schowalter, Arya and Proctor2000). However, unlike linear instability mechanisms – where perturbations are destined to be predominant in the limit of $t \rightarrow \infty$ due to the most unstable eigenmode’s exponential growth – the transient growth process generally necessitates a specific (optimal) form of perturbation as input. Whether such a specific perturbation can spontaneously emerge from ambient turbulence, which is fundamentally stochastic and uncontrolled, has led to recurrent criticisms of optimal transient growth (see Fontane et al. Reference Fontane, Brancher and Fabre2008, p. 235). In the work by Fontane et al. (Reference Fontane, Brancher and Fabre2008), where the transient dynamics of vortices with stochastic forcing was examined, optimal perturbations were shown to be activated by noise-like forcing that is random in both space and time. This mitigates the aforementioned criticisms of optimal transient growth. Nonetheless, as the authors stated, it remains questionable whether such random isotropic forcing effectively represents turbulence in real conditions. The absence of clear universality in modelling turbulence poses a significant challenge in integrating ambient turbulence into this current problem.

Instead, we take an initiative in considering a different, non-stochastic means of initiating optimal transient growth: ice crystals (or particles). The presence of such particles in real-world scenarios is ascertained through contrail observations. Contrail formation primarily starts with jet exhaust plumes produced by aircraft engines, which contain particulate matter that eventually acts as condensation nuclei (Kärcher Reference Kärcher2018). However, our interest does not lie in this very initial contrail stage during jet plume development, when the vortex roll-up process is still ongoing. According to early experimental findings by El-Ramly & Rainbird (Reference El-Ramly and Rainbird1977), there is no appreciable influence of the engine exhaust on altering the rolled-up structure at this stage. As a matter of course, our attention is directed towards the later stage involving a fully formed wake vortex, where interactions between the vortex and ice crystals become manifest.

During the stage when ice crystals interact with the wake vortex, individual particle sizes reach a few microns, or approximately 1–5 microns (Kärcher et al. Reference Kärcher, Peter, Biermann and Schumann1996; Paoli & Garnier Reference Paoli and Garnier2005; Naiman et al. Reference Naiman, Lele and Jacobson2011; Voigt et al. Reference Voigt, Schumann, Jessberger, Jurkat, Petzold, Gayet, Krämer, Thornberry and Fahey2011; Kärcher Reference Kärcher2018). It is noted that, according to Voigt et al. (Reference Voigt, Schumann, Jessberger, Jurkat, Petzold, Gayet, Krämer, Thornberry and Fahey2011), actual contrail samples exhibited particle sizes ranging from 0.39 to 17.7 microns, with an effective radius of 2.9 microns. Given the substantial size difference between the vortex scale (measured in metres) and the ice crystals (microns), these particles are often treated as flow tracers, i.e. with no backward influence on the carrier fluid (Paoli & Garnier Reference Paoli and Garnier2005; Naiman et al. Reference Naiman, Lele and Jacobson2011). This assumption can reasonably simplify the particle-flow dynamics. However, when considering a large particle number density, reportedly ranging from $10^{9}$ to $10^{11}$ per cubic metre (Paoli et al. Reference Paoli, Hélie and Poinsot2004; Paoli & Garnier Reference Paoli and Garnier2005), combined with the density ratio of ice to air (approximately $10^{3}$ ), their bulk impact may not be simply ignored. In the most optimistic estimation based on the given figures (a particle radius of 5 microns and a particle number density of $10^{11}$ per cubic metre), the upper limit of the particle mass fraction could fall between $10^{-2}$ and $10^{-1}$ , which, although representing one extreme end, is not preposterous. This mass fraction is substantial enough to initiate perturbations ultimately evolving into significant disturbances (recall § 3.4).

When considering ice crystals, a primary emphasis has typically been on their microphysical growth, typically using an ice microphysics model alongside a flow solver (e.g. Lewellen & Lewellen Reference Lewellen and Lewellen2001; Paoli et al. Reference Paoli, Hélie and Poinsot2004; Paoli & Garnier Reference Paoli and Garnier2005; Naiman et al. Reference Naiman, Lele and Jacobson2011). We exclude the context of microphysical growth in this study and, instead, direct our attention to two-way coupling, specifically through drag momentum exchange. This direction aligns with our essential goal of exploring whether particle-induced drag can initiate optimal transient growth.

Inspired by natural contrail formation, our analysis in this section seeks to shed light on how particles’ backward influence can act as a ‘controlled’ perturbation mechanism to trigger optimal transient growth. Unlike ambient turbulence, which is stochastic and uncontrollable, particles may offer the potential for a more deliberate approach. By strategically releasing particles around the vortex periphery, e.g. adjusting their distance from the vortex or timing their ejection period, it is possible to introduce targeted perturbations that invoke the optimal structures identified for transient growth.

Recent studies by Shuai & Kasbaoui (Reference Shuai and Kasbaoui2022) and Shuai et al. (Reference Shuai, Jeswin Dhas, Roy and Kasbaoui2022) indicate that weakly inertial particles within a vortex, under two-way coupled conditions, can trigger instabilities and expedite the process of vortex decay. While the initial particle distributions considered in these studies – where particles are loaded either throughout the entire domain or inside the vortex core – are not directly applicable to our case, where particles are at the periphery of the vortex core and interact with the vortex, these studies support the underlying concept that even a dilute amount of particles can meaningfully influence the surrounding vortex. By considering particles, our study aims at broadening the understanding of transient growth dynamics in a controlled manner, providing insights into potential applications of particle-driven perturbations.

4.2. Two-way coupled equations for vortex–particle interaction

To simulate the initiation process leading to transient vortex growth via particle drag, it is necessary to introduce additional parameters and variables to establish the equations governing particle motion. Additionally, a coupling term must be added to the momentum equation of fluid motion in (2.2) to complete a two-way coupled formulation. In this study, we consider dispersed ice crystals, whose density is approximately $10^3$ times that of the surrounding fluid (air). We define the ratio of particle density $\rho _{p}$ to fluid density $\rho$ as a new dimensionless parameter, denoted by $\vartheta$ ( ${\equiv}\rho _{p} / \rho$ ). For subsequent calculations, we set $\vartheta$ to a constant value of $10^3$ .

In this study, we employ the Eulerian approach adopting the fast equilibrium approximation proposed by Ferry & Balachandar (Reference Ferry and Balachandar2001). This method treats the set of particles as a continuum, allowing the flow–particle system to behave like a two-phase flow. Given the high particle number density, we opt against Lagrangian approaches that track all individual particles (e.g. Paoli et al. Reference Paoli, Hélie and Poinsot2004; Naiman et al. Reference Naiman, Lele and Jacobson2011; Shuai et al. Reference Shuai, Jeswin Dhas, Roy and Kasbaoui2022). Thanks to the relatively moderate computational cost, we believe that the Eulerian approach is favourable for future scale-up simulations involving two or more vortices to explore secondary vortex evolution. Moreover, our focus on the ‘bulk’ influence of particles on the surrounding vortex, rather than individual particle statistics, justifies the continuum treatment of particles.

The following two variables now represent the particles in the dispersed phase form: the particle velocity field, $\boldsymbol{u}_p (r,\phi ,z, t)$ , and the particle volume fraction, $c (r,\phi ,z, t)$ . The fast equilibrium approximation enables explicit evaluation of $\boldsymbol{u}_p$ in terms of the fluid velocity field, $\boldsymbol{u}$ . Using the Maxey–Riley equation with the added mass effect (Maxey & Riley Reference Maxey and Riley1983; Auton, Hunt & Prud’Homme Reference Auton, Hunt and Prud’Homme1988), $\boldsymbol{u}_p$ can be reduced and the resulting two-way coupled equations (see Ferry & Balachandar Reference Ferry and Balachandar2001, p. 1221) are

(4.1) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} = - \boldsymbol{\nabla } \varphi + \boldsymbol{u} \times \boldsymbol{\omega } + \frac {1}{{\textit {Re}}} {\nabla }^2 \boldsymbol{u} - (\vartheta - 1)c \frac {D\boldsymbol{u}}{Dt} \quad \text{with}\ \boldsymbol{\nabla } \cdot \boldsymbol{u}=0, \end{equation}

and

(4.2) \begin{equation} \frac {\partial c}{\partial t} = - \boldsymbol{u} \cdot \nabla c + \frac {2 {Stk} (\vartheta - 1)}{2\vartheta + 1} \nabla \cdot \left ( c \frac {D\boldsymbol{u}}{Dt} \right ), \end{equation}

where $D/Dt \equiv \partial /\partial t + \boldsymbol{u} \cdot \nabla$ is the material derivative with respect to the fluid phase and $Stk$ is the Stokes number, i.e. the dimensionless particle relaxation time normalised by $R_0/U_0$ . For calculations, $Stk$ is set to $10^{-5}$ to maintain compatibility with the fast equilibrium approximation as well as practical conditions (Kärcher et al. Reference Kärcher, Peter, Biermann and Schumann1996). The discretisation and time-integration procedures for (4.1) and (4.2) are not different from those used in the previous pure vortex cases, as detailed in Appendix B. Compared with (2.2), it is found that the last term in (4.1) accounts for the particle drag, with its magnitude dependent on the order of $c$ .

4.3. Initial particle distribution

As discussed in recent studies on vortex–particle interactions (Shuai & Kasbaoui Reference Shuai and Kasbaoui2022; Shuai et al. Reference Shuai, Jeswin Dhas, Roy and Kasbaoui2022), the observable effects of these interactions depend on the initial particle distribution. To explore whether particle drag can initiate optimal transient growth in the vortex, it is necessary to establish an effective initial distribution of the particles. In the following discussion, we first emphasise that the particles’ influence is limited to generating only a weak disturbance within the vortex system, under the assumption that $(\vartheta - 1)c$ remains significantly less than order unity, comparable to the perturbation order.

To comprehend how the particles induce perturbations in the carrier fluid, we reorganise the coupled momentum equation in (4.1), yielding

(4.3) \begin{equation} \left ( 1 + (\vartheta - 1 ) c \right ) \frac {\mathrm{D} \boldsymbol{u}}{\mathrm{D} t} = - \boldsymbol{\nabla } p + \frac {1}{{\textit {Re}}} {\nabla }^2 \boldsymbol{u}, \end{equation}

which is in fact the original formulation provided by Ferry & Balachandar (Reference Ferry and Balachandar2001). In this form, it is evident that the combined motion of the two phases resembles a single-phase flow with slight density variations, as characterised by the factor $(1 + (\vartheta - 1) c)$ . This effect is understood as a consequence of the dispersed phase absorbing momentum from the fluid phase. Loosely speaking, if the fluid accelerates, the particles’ presence retards the fluid’s acceleration, resulting in a negative perturbation velocity, and vice versa.

Given that the perturbation velocity is induced by non-zero $c$ , we derive the perturbation momentum equation from the decompositions of velocity and pressure in (4.3) (i.e. $\boldsymbol{u} = \overline {\boldsymbol{u}} + \boldsymbol{u}'$ and $p = \overline {p} + p'$ ), with the assumption that $(\vartheta - 1 ) c$ is comparable to the perturbation (prime) order, e.g. $O(\varepsilon )$ . Since zero particle concentration implies no perturbation, the mean quantity (overscore) equation is obtained by setting $c=0$ in (4.3), yielding ${\mathrm{D} \boldsymbol{\overline {u}}}/{\mathrm{D}t} = - \boldsymbol{\nabla }\, \overline {p} + \nabla ^2{\overline {\boldsymbol{u}}}/{{\textit {Re}}}$ . Subtracting this from (4.3) and neglecting terms of order higher than $O(\varepsilon )$ , we obtain

(4.4) \begin{equation} (\vartheta - 1 ) c \left [ \frac {\partial \overline {\boldsymbol{u}}}{\partial t} + (\overline {\boldsymbol{u}} \cdot \nabla )\,\overline {\boldsymbol{u}} \right ] + \frac {\partial \boldsymbol{u}'}{\partial t} + (\overline {\boldsymbol{u}} \cdot \nabla )\,\boldsymbol{u}'= - \boldsymbol{\nabla } p' + \frac {1}{{\textit {Re}}} {\nabla }^2 \boldsymbol{u}'. \end{equation}

Our goal is to estimate the initial particle distribution $c(t=0)$ , denoted $c_0$ , that effectively perturbs the ‘undisturbed’ vortex towards optimal transient growth, i.e. $\boldsymbol{u}'(t=0) = 0$ . Under this condition, (4.4) reduces to the initial perturbation growth formula:

(4.5) \begin{equation} \frac {\partial \boldsymbol{u}'}{\partial t}\Big \rvert _{\, t=0} = - \boldsymbol{\nabla } p' (t=0) -(\vartheta - 1 ) c_0 \frac {\mathrm{D}\overline {\boldsymbol{u}}}{\mathrm{D}t}\Big \rvert _{\, t=0}, \end{equation}

which describes the growth of $\boldsymbol{u}'$ over a brief initial period, e.g. $0 \leqslant t \lt \delta t \ll 1$ . After a brief advance of time, the spatial structure of $\boldsymbol{u}'$ is primarily attributed to that of $-(\vartheta - 1 ) c_0\mathrm{D}\overline {\boldsymbol{u}}/\mathrm{D}t\rvert _{\, t=0}$ (initial particle drag) while $- \boldsymbol{\nabla } p' (t=0)$ serves only to maintain the divergence-free constraint. In other words, we expect $\boldsymbol{u'}(t=\delta t) \approx \delta t [ - (\vartheta - 1 ) c_0\mathrm{D}\overline {\boldsymbol{u}}/\mathrm{D}t\rvert _{\, t=0} ]$ with a minor adjustment for the field to be solenoidal (for details of the adjustment, see Appendix C).

Suppose that we aim to initiate a known optimal velocity perturbation $\breve {\boldsymbol{u}}'$ . If we identify $c_0$ such that $-(\vartheta - 1 ) c_0\mathrm{D}\overline {\boldsymbol{u}}/\mathrm{D}t\rvert _{\, t=0}$ matches (a positive constant multiple of) $\breve {\boldsymbol{u}}'$ , then we can expect $\boldsymbol{u}'$ to exhibit this perturbation form. However, the existence of such $c_0$ is rare, since this matching requires solving three component equations whereas $c_0$ is the only unknown. To circumvent this overdetermination issue, we inevitably focus on the most critical component among them. By choosing the azimuthal component, the problem reduces to finding $c_0$ such that

(4.6) \begin{equation} -(\vartheta - 1) c_0 { \left ( \frac {\mathrm{D}\overline {\boldsymbol{u}}}{\mathrm{D}t} \right )_\phi \, \Big \rvert _{\, t=0}} = { C \breve {u}'_\phi }, \end{equation}

where $C$ is an arbitrary positive constant, which can be used to scale $c_0$ when solving (4.1) and (4.2) for different levels of particle volumetric loading. Arranging the terms based on the fact that $(\mathrm{D}\overline {\boldsymbol{u}}/\mathrm{D}t)_\phi \rvert _{t=0} = -4r e^{-r^2} / {\textit {Re}}$ for the ‘undisturbed’ $q$ -vortex profile, we can further simplify this relation to $c_0 \propto \breve {u}'_\phi /(r e^{-r^2})$ .

The proposed $c_0$ has two serious limitations that restrict its practical applicability. Nevertheless, we affirm that it remains sufficiently effective for initiating optimal transient growth. First, the radial and axial components of perturbation velocity are excluded in this derivation. This is justifiable since the azimuthal component of velocity perturbation is generally dominant during the initial stages of transient growth (see the left panels in figure 7). Second, the particle volume fraction cannot take on negative values. The continuity of the fluid may alleviate this issue; in a local sense, a deficit (or surplus) in speed within the particle-laden fluid must be offset by a corresponding gain (or loss) of speed in the circumferential particle-free fluid. As a result, we set $c_0$ to locally zero when its calculated value from (4.6) is locally negative.

For comparison’s sake, we revisit the optimal perturbation case discussed in § 3.4, where $(m,\,\kappa ) = (1,\,0.1)$ with $\tau = 50$ . Figure 12 shows the initial particle volume fraction calculated using the proposed estimation, along with the axial vorticity after a brief advancement in time ( $t=0.01$ ), resulting from the two-way interactions between the particles and the vortex with an initial particle volumetric loading level of $c_{max } = 10^{-6}$ . The computation of perturbation velocity fields in two-way coupled vortex–particle simulations follows the same approach as described in (3.5), except that what determines the perturbation intensity is now the particle volumetric loading level $c_{max }$ . Despite the previously mentioned limitations, the resulting perturbation effectively resembles the intended perturbation input (see figure 8) needed for optimal transient growth. This observation provides a posteriori validation of (4.6).

Figure 12. Initiation of optimal transient growth via inertial particles: $(a)$ the initial particle volume fraction contour on the $z=0$ plane, in the pursuit of initiating the perturbation discussed in § 3.4 (see figure 8); and $(b)$ the axial vorticity perturbation contour on the $z=0$ plane after a brief advancement in time ( $t=0.01$ ) in the two-way coupled vortex–particle simulation, solving (4.1) and (4.2) with an initial particle volumetric loading level of $c_{max } = 10^{-6}$ . In panel (b), the same contour style as figure 8 is used and $|\omega _z' (t=0.01)|_{max } = 5.80 \times 10^{-6}$ .

4.4. Particle-initiated transient growth

Now that the vortex evolution is governed by (4.1), where the term $(\vartheta - 1) c (\mathrm{D} {\boldsymbol{u}}/\mathrm{D}t)$ serves as an additional forcing term, the overall perturbation dynamics is influenced not only by the transient growth process but also by the continual interaction between the particles and the vortex flow over time (see also a simliar discussion by Fontane et al. Reference Fontane, Brancher and Fabre2008, p. 249). In what follows, we substantiate particle-initiated transient growth by identifying some notable indications of transient growth over short time intervals in the vortex–particle system, using the initial particle distribution shown in figure 12.

Temporal changes in perturbation energy are displayed in figure13 for four different levels of particle volumetric loading: $c_{max } = 10^{-4}$ , $10^{-5}$ , $10^{-6}$ and $10^{-7}$ . The case of $c_{max } = 10^{-4}$ is considered to be the upper limit where $(\vartheta - 1) c$ remains significantly less than order unity (n.b. $\vartheta = 10^3$ ). All energy curves exhibit a nearly identical trend, suggesting that similar dynamics govern every case. In panel (b), the data are normalised by the perturbation energy at $t=10$ for each case (note that $E(0)$ is zero and the energy growth used, $E(t)/E(0)$ , becomes undefined here), allowing for a comparison of energy amplification across cases. Overall, energy amplification levels off at around $10^2$ times $E(10)$ , due to the long-term response of the vortex to the particle interactions. Arguably, amplification beyond this level should be attributable to the transient growth process, as evidenced by the energy amplification ‘hump’ up to $t = 80$ . Notably, the peak of this hump at $t=50$ coincides with the optimal transient growth period, $\tau = 50$ , as intended to induce, which further supports our argument.

Figure 13. (a) Temporal changes in perturbation energy for vortex–particle interactions (see figure 12 for the initial set-up), evaluated at various levels of particle volumetric loading: $c_{max } = 10^{-4}$ , $10^{-5}$ , $10^{-6}$ and $10^{-7}$ ; and (b) the same data, normalised by $E(10)$ , to facilitate comparison of energy amplification across cases. Note that $E(10)$ is used for normalisation because $E(t)/E(0)$ is undefined here ( $E(0) = 0$ ).

Figure 14. Axial vorticity perturbation contours of the vortex interacting with particles initially distributed around the periphery (see figure 12 for the initial particle distribution) on the $z=0$ plane at $t=25$ , $t=50$ and $t=100$ . Here, depicted is the case with $c_{max } = 10^{-4}$ .

Another indication of the particle-initiated transient growth is observed in the evolution of the perturbation structure. In figure 14, axial vorticity perturbation contours on the $z=0$ plane at $t=25$ , $t=50$ and $t=100$ in the vortex–particle simulation with $c_{max } = 10^{-4}$ are shown. These snapshots are compared with those of the optimally perturbed nonlinear vortex growth in figure 10 $(c)$ . The continual vortex–particle interactions produce structural discrepancies in the perturbation; this is clearly discernible at $t=100$ , where strong spiralling arms form at the periphery as a result of prolonged drag momentum exchange. Nonetheless, during the early-stage perturbation growth up to $t=50$ , key features characterising the optimal transient growth process are observed, such as the appearance of two weak spiralling arms at the periphery of the core at $t=25$ . Most notably, at the time of maximum energy growth ( $t=50$ ), perturbation energy transfer from the periphery to the core – the iconic feature of the optimal transient growth process for $m=1$ – is identifiable. We believe that this provides plausible evidence that near-optimal transient growth takes place through vortex–particle interactions.

Lastly, we report the transient development of the particle distribution in association with the vortex’s transient growth. In figure 15, the interaction between the vortex and particles for the case $c_{max } = 10^{-4}$ is visualised, using the $\lambda _2$ -isosurface at $\lambda _2 = -0.05$ for the vortex core and the $c$ -isosurface representing 20 % of $c_{max }$ for the particles. As the vortex evolves from its unperturbed state at $t=0$ to its most excited state at $t=50$ , the particles show reduced dispersion, forming a coherent helical structure that envelops the vortex core. However, this coherence is evanescent and dissipates rapidly after the peak perturbation growth at $t=50$ . The physical significance of this temporary coherence increase, as well as whether this is exclusive to the $m=1$ case, warrants further investigation. For now, it is noteworthy that the physical phenomenon linked to the current example, vortex meandering, is known to increase system ‘orderliness’ (i.e. reduce the number of dynamically active proper orthogonal decomposition (POD) modes), as reported by Bölle (Reference Bölle2021); the coherence of the particles during transient vortex growth may be associated with this tendency.

Figure 15. Three-dimensional illustration of the $q$ -vortex interacting with the peripherally located particles (see figure 12) of the initial particle volumetric loading level of $c_{max } = 10^{-4}$ : $(a)$ $t=0$ and $(b)$ $t=50$ . The vortex core is detected using the $\lambda _2$ -isosurface where $\lambda _2 = -0.05$ , as with figure 11. The green isosurface of $c = 2 \times 10^{-5}$ (20 % of $c_{max }$ ) is drawn together to visualise the particle distribution at each time.