# Mesh-Free Hydrodynamic Stability

Tianyi Chu<EMAIL_ADDRESS>Oliver T. Schmidt<EMAIL_ADDRESS>Department of
Mechanical and Aerospace Engineering, University of California San Diego, 9500
Gilman Drive, La Jolla, CA 92093, USA.

###### Abstract

A specialized mesh-free radial basis function-based finite difference (RBF-FD)
discretization is used to solve the large eigenvalue problems arising in
hydrodynamic stability analyses of flows in complex domains. Polyharmonic
spline functions with polynomial augmentation (PHS+poly) are used to construct
the discrete linearized incompressible and compressible Navier-Stokes
operators on scattered nodes. Rigorous global and local eigenvalue stability
studies of these global operators and their constituent RBF stencils provide a
set of parameters that guarantee stability while balancing accuracy and
computational efficiency. Specialized elliptical stencils to compute boundary-
normal derivatives are introduced and the treatment of the pole singularity in
cylindrical coordinates is discussed. The numerical framework is demonstrated
and validated on a number of hydrodynamic stability methods ranging from
classical linear theory of laminar flows to state-of-the-art non-modal
approaches that are applicable to turbulent mean flows. The examples include
linear stability, resolvent, and wavemaker analyses of cylinder flow at
Reynolds numbers ranging from 47 to 180, and resolvent and wavemaker analyses
of the self-similar flat-plate boundary layer at a Reynolds number as well as
the turbulent mean of a high-Reynolds-number transonic jet at Mach number
$0.9$. All previously-known results are found in close agreement with the
literature. Finally, the resolvent-based wavemaker analyses of the Blasius
boundary layer and turbulent jet flows offer new physical insight into the
modal and non-modal growth in these flows.

###### keywords:

RBF-FD , polyharmonic splines , polynomial augmentation , linear stability
theory, Navier-Stokes , wavemaker

##  1 Introduction

Flow instabilities and large-scale coherent structures are ubiquitous
phenomena in fluid mechanics that have been the focus of extensive research.
Linear stability (LST) analysis is specifically designed to investigate the
growth of small perturbations exclusively around laminar base flows, which are
the steady-state solution to the Navier-Stokes equations. One-dimensional LST
analysis was widely used in the past century, e.g., [20, 30, 52, 79]. Eriksson
and Rizzi [36] and Tuckerman and Marcus [146] were the pioneers in conducting
LST analysis in two-dimensional (2D). Subsequently, Jackson [57] and Zebib
[154] examined the 2D nature of vortex shedding in the wakes of bluff bodies.
Readers are referred to Huerre and Monkewitz [53] and Theofilis [136] for
comprehensive reviews of the concept of 2D LST modes. The implementation of
the 2D LST framework has enabled improved identification of flow instability
in non-parallel flows, including cylinder wakes [46, 76, 91], aerofoil wakes
[34, 151], boundary layers [1, 35], and jets in cross-flow [10, 95, 107].
Studies have demonstrated that an open flow can possess marginal stability
despite exhibiting local convective instability. LST analysis around a steady
laminar base flow, by its inherent nature, is not applicable to predict
finite-amplitude flow instabilities arising from nonlinear interactions. The
use of mean flow for LST analysis, despite violating the basic assumption of
linear theory, has been used successfully to predict coherent flow features in
diverse types of flows, including cylinder wakes [11, 100], open cavity flows
[126, 127], mixing layers [45, 87, 149], and turbulent or transitional jets
[48, 92, 116]. Theoretical conditions required for the validity of mean flow
stability analysis have been explored by Beneddine _et al._ [18]. Although
beyond the scope of this work, it is noteworthy that the weakly-nonlinear
extension of LST analysis has been successfully employed for studying the
dynamics of non-parallel flows [27, 127]. The LST-based semi-empirical $e^{N}$
method [55, 131] has succeeded in transition prediction for certain flows such
as boundary layers. However, the prediction of disturbance behavior in more
complex scenarios, such as crossflows or bypass transitions, falls outside the
scope of LST theory.

Despite these limitations, resolvent, or input-output analysis, has recently
emerged as a linear tool for accurately predicting large-scale coherent
structures in fully turbulent flows. Resolvent analysis (RA) originally stems
from the studies of transient growth [37, 104, 105, 144] and seeks the optimal
pairs of inputs and corresponding outputs through the linearized system that
maximizes the energy gain. Within the laminar regime, RA has been utilized to
investigate the linear response to external body forces or perturbations for
channel flows [64, 84, 111], boundary layers [4, 21, 24, 89, 97, 128, 142],
and jets [44, 59, 118]. In contrast to classical LST analysis, the input-
output perspective offers a mathematically rigorous framework for studying
turbulent mean-flows by identifying the forcing as the Reynolds stresses in
the perturbation-interaction terms in the Reynolds-decomposed Navier–Stokes
equations [81, 124]. Applications include near-wall flows [54, 1, 124],
boundary layers [6, 58, 110], incompressible [44, 69] or compressible jets
[59, 117, 143], and airfoil wakes [108, 153]. The validation of fundamental
relationships between RA and other modal decomposition techniques was
facilitated by Towne _et al._ [143], establishing RA as a well-suited tool for
turbulence modeling [98, 99]. Both classical LST and RA require, in their most
basic form, the construction of large matrices that have to be decomposed into
their singular- or eigen-components. The construction and decomposition of
these matrices are particularly challenging for flows in complex geometries.
Previous studies [33, 109, 125, 153] have leveraged the flexibility of the
finite-volume (FV) methods on unstructured meshes. However, a downside of
unstructured FV methods is that the accuracy is usually restricted to 2nd-
order. Alternatively, the utilization of finite-element (FE) methods for
spatial discretization provides high-order accuracy for flow instability
analysis [76, 77, 110, 127]. The commonly employed weak formulation of
governing equations in FE methods raises concerns regarding stability and
convergence. Finite-element (FE) methods offer the same flexibility, and the
FreeFEM+ toolbox by Hecht [51] has been employed in a number of recent studies
[76, 77, 110, 127]. Matrix-free methods, such as the time-stepper [7, 10, 145]
and other related techniques [32, 78, 89, 93], provide an alternative approach
where the decomposition of large matrix operations can be completely
circumvented. Iterative Krylov subspace methods are commonly employed to
obtain a partial eigendecomposition. However, their major limitation lies in
their capacity to extract only a limited portion of the spectra at a time.
Randomized approaches have also been explored as potential solutions to
decrease the computational cost of the singular- or eigendecomposition of
large stability matrices [85, 108]. In this study, we demonstrate the
capability and accuracy of radial basis functions (RBF) in effectively
addressing the aforementioned challenges.

The use of RBF-based methods provides high flexibility in meshing complex
geometries, allowing for local grid refinement and local adjustments to the
order of accuracy of the discretization. The RBF methodology is rooted in
scattered data interpolation, first introduced by Hardy [50], and offers a
systematic means of approximating multivariate functions on arbitrarily
scattered nodes. By generalizing the classical finite-difference (FD) methods,
the so-called RBF-FD methods enable the approximation of spatial derivatives
to arbitrary node layouts based on RBFs [139]. In fluid flow problems,
Gaussian (GA), multiquadric (MQ), and inverse multiquadric (IMQ) are among the
most commonly used radial basis functions (RBFs). Refer to Fornberg and Flyer
[42] for a detailed overview. The computation of these RBFs requires
specifying a shape parameter, a free parameter that can significantly impact
both the numerical stability and accuracy. Recently, a new class of RBF-FD
approximations has emerged that utilizes polyharmonic splines (PHS) augmented
with polynomials (PHS+poly). This technique was first introduced by Flyer _et
al._ [39] and has demonstrated great potential in achieving high-order
accuracy while avoiding the difficulties associated with the shape parameter
selection. Bayona _et al._ [15] and Flyer _et al._ [40] highlighted the
advantages of local PHS+poly-generated RBF-FD stencils in achieving high-order
accuracy for derivative approximations and solving elliptic partial
differential equations (PDEs), while also eliminating stagnation errors under
node refinement. The feasibility of implementing a larger stencil near domain
boundaries to avoid the Runge phenomenon was analytically verified by Bayona
[13] using RBFs in a closed-form. Additional numerical examples in 2-D and 3-D
were presented in Bayona _et al._ [16] to demonstrate the effectiveness of the
PHS+poly approach. The accuracy and stability of the PHS+poly RBF-FD method
depend on the combination of the stencil size, PHS exponent, and the degree of
polynomials. Several studies, including Chu and Schmidt [29], Le Borne and
Leinen [67], Shankar and Fogelson [123], have investigated and discussed the
optimal combinations of these parameters to achieve the ‘spectral-like’
accuracy of higher-order Padé-type compact FDs [68]. In addition, other
approaches have been explored to improve the computational accuracy and
efficiency of PHS+poly RBF-FDs, including overlapping stencils [122, 123], the
partition of unity (PU) method [83], and oversampling [140]. Comparative
analyses have been carried out on various higher-order mesh-free
discretizations, consistently demonstrating that PHS+poly-based RBF-FDs
outperform other methods in accuracy, robustness, and computational
efficiency, e.g., [14, 40, 47, 94, 112]. Recent successful applications of
PHS+poly RBF-FDs in flow simulations on scattered nodes [119, 120, 121, 147]
highlight the potential of RBF methods for investigating flow instabilities in
complex geometries. In our recent work [29], we introduced a staggered-node
RBF-FD-based semi-implicit fractional-step solver for the incompressible
Navier-Stokes equations, which will be utilized to compute the base state in
the present study.

This work introduces a high-order, mesh-free framework for hydrodynamic
stability analysis based on the PHS+poly RBF-FD discretizations. The efficient
construction of global Jacobians on scattered nodes using RBF-based methods is
demonstrated. The trade-off between computational efficiency and accuracy is
addressed, and best practices are provided. Furthermore, we discuss the
practical treatment of boundary conditions, particularly the boundary-normal
derivatives and the pole singularity in cylindrical coordinates. The numerical
framework is demonstrated on various hydrodynamic stability theoretical
methods and flows. Refer to table 3 for an overview. In particular, LST, RA,
and wavemaker (WM) analyses are shown for the canonical cylinder flow at
Reynolds numbers ranging from 47 to 180, RA and WM analyses for a laminar
zero-pressure-gradient (ZPG) flat-plate Blasius boundary layer at a Reynolds
number of $6\times 10^{5}$, RA and WM analyses for a turbulent transonic jet
at Mach number $0.9$ and a Reynolds number of approximately $10^{6}$. The
identified flow phenomena include vortex shedding in cylinder wakes [11, 100],
the Orr and Tollmien-Schlichting waves in boundary layers [1, 22, 89, 128],
and Orr and Kelvin-Helmholtz wavepackets [26, 44, 69, 117, 134] as well as
trapped acoustic waves [116, 141] in turbulent jets. The comparisons of these
benchmark problems with the literature allow us to highlight the viability,
accuracy, and robustness. The study furthermore encapsulates the first
application of RA-based WM analysis to the ZPG Blasius boundary layer and
turbulent jets, providing a new perspective on modal and non-modal growth in
these flows.

This paper is organized as follows: §2 introduces the PHS+poly RBF-FD method,
§3 presents the mesh-free hydrodynamic stability analysis, and §4 discusses
the numerical implementations. The performance of the mesh-free framework is
demonstrated in §5. Finally, §6 concludes and summarizes the paper.

##  2 Radial basis function-based finite differences (RBF-FD)

The goal of the RBF-FD method is to compute the discrete representation of any
linear differentiation operator $\mathcal{D}$ at a given location $x_{0}$ as a
linear combination of the function values $g({\bf\it x}_{j})$, such that

$\mathcal{D}g({\bf\it x}_{0})=\sum_{j=1}^{n}w_{j}g({\bf\it x}_{j}).$ (1)

To obtain the unknown differentiation weights $w_{j}$, we use the RBF
interpolant,

$s({\bf\it x})=\sum_{j=1}^{n}\gamma_{j}\phi(\|{\bf\it x}-{\bf\it x}_{j}\|),$
(2)

to approximate the given function $g({\bf\it x})$ by satisfying

$s({\bf\it x}_{j})=g({\bf\it x}_{j}),\qquad j=1,2,\cdots n.$ (3)

Here, $\phi(r)$ is a smooth radial function, $\\{{\bf\it x}_{j}\\}_{j=1}^{n}$
is a set of scattered nodes, and $\|\cdot\|$ denotes the standard Euclidean
norm. Combining equations (1-3) leads to the linear system

$\underbrace{\matrixquantity[\phi(\|{\bf\it x}_{1}-{\bf\it
x}_{1}\|)&\phi(\|{\bf\it x}_{1}-{\bf\it x}_{2}\|)&\cdots&\phi(\|{\bf\it
x}_{1}-{\bf\it x}_{n}\|)\\\ \phi(\|{\bf\it x}_{2}-{\bf\it
x}_{1}\|)&\phi(\|{\bf\it x}_{2}-{\bf\it x}_{2}\|)&\cdots&\phi(\|{\bf\it
x}_{2}-{\bf\it x}_{n}\|)\\\ \vdots&\vdots&&\vdots\\\ \phi(\|{\bf\it
x}_{n}-{\bf\it x}_{1}\|)&\phi(\|{\bf\it x}_{n}-{\bf\it
x}_{2}\|)&\cdots&\phi(\|{\bf\it x}_{n}-{\bf\it x}_{n}\|)]}_{{\bf\it
A}}\matrixquantity[w_{1}\\\ w_{2}\\\ \vdots\\\
w_{n}]=\matrixquantity[\evaluated{\mathcal{D}\phi(\|{\bf\it x}-{\bf\it
x}_{1}\|)}_{{\bf\it x}={\bf\it x}_{0}}\\\ \evaluated{\mathcal{D}\phi(\|{\bf\it
x}-{\bf\it x}_{2}\|)}_{{\bf\it x}={\bf\it x}_{0}}\\\ \vdots\\\
\evaluated{\mathcal{D}\phi(\|{\bf\it x}-{\bf\it x}_{n}\|)}_{{\bf\it x}={\bf\it
x}_{0}}],$ (4)

which can be solved directly to obtain the weight vector ${\bf\it
w}=[w_{1},\cdots,w_{n}]^{T}$. An implicit assumption for equation (4) is that
the derivative of the basis function, $\mathcal{D}\phi$, is continuous.

Polynomial augmentation is commonly applied to the RBF-FD method to enforce
consistency with Taylor expansion-based FD approximations [39, 42, 43, 66,
152]. The two-dimensional augmented RBF-FD method can be expressed as

$\mathcal{D}g({\bf\it x}_{0})=\sum_{j=1}^{n}w_{j}g({\bf\it
x}_{j})+\sum_{i=1}^{(q+1)(q+2)/2}c_{i}P_{i}({\bf\it x}_{0}),$ (5)

where $P_{i}({\bf\it x})$ are multivariate polynomials up to degree $q$. To
match the local Taylor series, additional constraints for the differentiation
weights,

$\displaystyle\sum_{j=1}^{n}w_{j}P_{i}({\bf\it
x}_{j})=\mathcal{D}P_{i}({\bf\it x}_{0})\qquad\text{for}\,\,1\leq
i\leq\frac{(q+1)(q+2)}{2},$ (6)

are included in the computation. These constraints, also referred to as the
vanishing momentum conditions [56], ensure that the RBF approximations locally
reproduce polynomial behavior up to degree $q$ [40] and decay in the far-field
[41]. A more general and compact representation of equations (5-6) are

$\underbrace{\left[\begin{array}[]{c c }{\bf\it A}&{\bf\it P}\\\ {\bf\it
P}^{\small{T}}&{\bf\it 0}\end{array}\right]}_{{\bf\it
A}_{\text{aug}}}\matrixquantity[{\bf\it w}\\\ {\bf\it
c}]=\matrixquantity[{\mathcal{D}{\bf\it\phi}}\\\ {\mathcal{D}{\bf\it P}}],$
(7)

where ${\bf\it A}$ represents the same interpolation matrix as defined in
equation (4). In equation (1), only the weight vector ${\bf\it w}$ is used to
approximate the differentiation operator $\mathcal{D}$. Motivated by recent
studies (see, e.g., [15, 16, 39, 40]), we use polyharmonic splines (PHS),

$\displaystyle\phi(r)=r^{m},$ (8)

as the basis functions, where $m$ is an odd positive integer. The accuracy and
stability of this PHS+poly RBF-FD method depend on the combination of the
stencil size, PHS exponent, and the degree of polynomials. The selection of
these parameters will be discussed later in §4.1.

RBF-FD methods provide a mesh-free approach for discretizing the domain
without the need for a pre-defined mesh or explicit node connectivity. This
feature enables the investigation of flow instabilities for complex
geometries, thereby offering advantages over traditional discretizations that
rely on structured meshes. The subsequent sections of this paper focus on the
construction of the discrete linearized Navier-Stokes (LNS) and then the
resolvent operators using PHS+poly RBF-FDs.

##  3 Mesh-free hydrodynamic stability analysis

###  3.1 Governing equations

The motion of a general incompressible Newtonian fluid is governed by the
Navier-Stokes equations,

$\displaystyle\partialderivative{{\bf\it u}}{t}+\left({\bf\it
u}\cdot{\nabla}\right){\bf\it u}$
$\displaystyle=-{\nabla}p+\mathrm{Re}^{-1}{\nabla}^{2}{\bf\it u},$ (9a)
$\displaystyle{\nabla}\cdot{\bf\it u}$ $\displaystyle=0.$ (9b)

Here, all variables are nondimensionalized by the velocity scale ${\bf\it
U}_{\infty}$ and the length scale $L$, and $\mathrm{Re}$ denotes the Reynolds
number.

In the case of laminar flows, we can decompose the flow state around the
steady-state solution of the Navier-Stokes equations (9) as

$\displaystyle{\bf\it u}={\bf\it U}+{\bf\it u}^{\prime},\quad
p={P}+p^{\prime},$ (10)

where $({\bf\it U},\,P)$ represents the base flow that satisfies

$\displaystyle\left({\bf\it U}\cdot{\nabla}\right){\bf\it U}$
$\displaystyle=-{\nabla}P+\mathrm{Re}^{-1}{\nabla}^{2}{\bf\it U},$ (11a)
$\displaystyle{\nabla}\cdot{\bf\it U}$ $\displaystyle=0,$ (11b)

and $(\cdot)^{\prime}$ denotes the small fluctuating components. In turbulent
flows, we take the Reynolds decomposition of the flow state into the temporal
mean, $\overline{\left(\cdot\right)}$, and fluctuating components, given by

$\displaystyle{\bf\it u}=\overline{{\bf\it u}}+{\bf\it u}^{\prime},\quad
p=\overline{p}+p^{\prime}.$ (12)

By generalizing the notation of the base state as $({\bf\it u}_{0},\,p_{0})$,
the resulting governing equations for the fluctuations take the form of

$\displaystyle\partialderivative{{\bf\it u}^{\prime}}{t}+\left({{\bf\it
u}_{0}}\cdot{\nabla}\right){\bf\it u}^{\prime}+\left({\bf\it
u}^{\prime}\cdot{\nabla}\right){{\bf\it u}_{0}}$
$\displaystyle=-{\nabla}p^{\prime}+\mathrm{Re}^{-1}{\nabla}^{2}{\bf\it
u}^{\prime}+{\bf\it f}^{\prime},$ (13a) $\displaystyle{\nabla}\cdot{\bf\it
u}^{\prime}$ $\displaystyle=0.$ (13b)

Here, the term ${\bf\it f}^{\prime}$ represents the remaining nonlinear
interactions between the fluctuation components. Table 1 summarizes the two
aforementioned decompositions.

Base state | Description | Notation | Obtained from | Remaning forcing ${\bf\it f}^{\prime}$
---|---|---|---|---
$({\bf\it u}_{0},\,p_{0})$ | Base flow | $({\bf\it U},\,P)$ | equation (11) | $-\left({\bf\it u}^{\prime}\cdot{\nabla}\right){\bf\it u}^{\prime}$
Mean flow | $(\overline{{\bf\it u}},\overline{p})$ | long-time average | $-\left({\bf\it u}^{\prime}\cdot{\nabla}\right){\bf\it u}^{\prime}+\overline{\left({\bf\it u}^{\prime}\cdot{\nabla}\right){\bf\it u}^{\prime}}$
Table 1: Summary of flow state decompositions.

Equation (13) can be written compactly in terms of the fluctuating state,
${\bf\it q}^{\prime}=[u^{\prime},\,v^{\prime},\,p^{\prime}]^{T}$, as

$\displaystyle\mathcal{P}\mathcal{P}^{T}\left(\partialderivative{t}{\bf\it
q}^{\prime}\right)$ $\displaystyle=\mathcal{L}{\bf\it
q}^{\prime}+\mathcal{P}{\bf\it f}^{\prime},$ (14)

where $\mathcal{P}$ is the prolongation operator that extends the velocity
vector $[u,\,v]^{T}$ into $[u,\,v,\,0]^{T}$, and its transpose is the
restriction operator that extracts the velocity vector from the extended state
vector [128]. The incompressible linearized Navier-Stokes (LNS) operator takes
the form of

The governing equations for compressible flows will be discussed in §5.3 and
A. Beyond the linear dynamics, the remaining forcing ${\bf\it f}^{\prime}$ in
equation (14) comprises products of fluctuating quantities, as outlined in
table 1. These terms will be either neglected or modeled.

Classical (temporal) linear stability (LST) analysis investigates fluctuations
with complex frequency $\lambda=\lambda_{r}+\mathrm{i}\lambda_{i}$, where
$\lambda_{r}$ is the exponential growth rate and $\lambda_{i}$ the oscillation
frequency. The fluctuations are assumed to be infinitesimally small, and the
forcing term, ${\bf\it f}^{\prime}$, is, therefore, negligible at $O(1)$, see,
e.g., Schmid and Henningson [115]. Substituting perturbations of the form
${\bf\it q}^{\prime}({\bf\it x},t)=\tilde{{\bf\it q}}({\bf\it
x})\mathrm{e}^{\lambda t}$ into the governing equations (13) yields the LST
equation,

$\displaystyle\lambda\mathcal{P}\mathcal{P}^{T}\tilde{{\bf\it
q}}=\mathcal{L}\tilde{{\bf\it q}}.$ (15)

Equation (15) is a generalized eigenvalue problem, and the eigenvector
associated with the largest growth rate ought to predict the dominant flow
instability mechanism.

The nonlinear interactions in equation (14) are no longer negligible for
general cases of finite amplitude fluctuations. Within the resolvent framework
of turbulent flows, the nonlinear interactions, along with the background
turbulence, can be interpreted as external forcing, ${\bf\it f}^{\prime}$, to
the otherwise linear dynamics. This interpretation was first proposed by
McKeon and Sharma [81]. By assuming a normal mode form for the fluctuating
components, $[{\bf\it q}^{\prime},\,{\bf\it f}^{\prime}]({\bf\it
x},t)=[\hat{{\bf\it q}},\,{\bf\it\hat{f}}]({\bf\it
x})\mathrm{e}^{\mathrm{i}\omega t}+c.c.$, where $\omega$ is the angular
frequency, or equivalently by taking the Fourier transform, we obtain the
linear time-invariant (LTI) representation of the governing equation (13) in
the frequency domain,

$\displaystyle\left(\mathrm{i}\omega\mathcal{P}\mathcal{P}^{T}-\mathcal{L}\right)\hat{{\bf\it
q}}$ $\displaystyle=\mathcal{P}\left(\mathcal{B}{\bf\it\hat{f}}\right),$ (16a)
$\displaystyle{\bf\it\hat{u}}$
$\displaystyle=\mathcal{P}^{T}\left(\mathcal{C}{\bf\it\hat{q}}\right).$ (16b)

The linear operators $\mathcal{B}$ and $\mathcal{C}$ are used to select
spatial regions of particular interest. We write equations (16) in a compact
form as

$\displaystyle\hat{{\bf\it u}}=\mathcal{H}(\omega){\bf\it\hat{f}},$ (17)

where
$\mathcal{H}(\omega)=\mathcal{P}^{T}\mathcal{C}\left(\mathrm{i}\omega\mathcal{P}\mathcal{P}^{T}-\mathcal{L}\right)^{-1}\mathcal{P}\mathcal{B}$
is known as the resolvent operator. The numerical discretization of the global
operator, $\mathcal{L}$, is at the core of both LST and resolvent analyses. In
this work, the PHS+poly RBF-FDs are utilized to construct differentiation
matrices and global Jacobians on scattered nodes, taking advantage of their
flexibility and accuracy in numerical discretization.

###  3.2 Global Jacobians and mesh-free stability analysis

For differentiation operators
$\mathcal{D}=\frac{\partial^{(\alpha+\beta)}}{\partial x_{1}^{\alpha}\partial
x_{2}^{\beta}}$ with $\alpha+\beta\leq 2$, we seek differentiation matrices
${\bf\it D}_{x_{1}^{\alpha}x_{2}^{\beta}}$ that satisfy

Here, $(x_{1},x_{2})$ represents the coordinate system, and $\\{{\bf\it
x}_{i}\\}_{i=1}^{N}$ represents the global computational grid. The $j$th row
of the sparse matrix ${\bf\it D}$ contains the $n\ll N$ weights that
approximate the derivative at node ${\bf\it x}_{j}$. The matrix hence has
$N\times n$ nonzero elements. In a slight change of notation, to enhance
readability, we now denote by ${\bf\it u}={\bf\it u}({\bf\it x})$ and ${\bf\it
v}={\bf\it v}({\bf\it x})$ the global velocity fields, and by ${\bf\it
p}={\bf\it p}({\bf\it x})$ the global pressure field in the computational
domain $\Omega$.

Upon the use of these global RBF-FD-based differentiation matrices, we
assemble the discrete global LNS operator from equation (3.1), taking the two-
dimensional Cartesian coordinate system as an example, as

where ${\bf\it S}=-\left({{\bf\it u}_{0}}\circ{\bf\it D}_{x}+{{\bf\it
v}_{0}}\circ{\bf\it D}_{y}\right)$, and $\circ$ denotes the Hadamard product.

For a weighted inner product $\left<{{\bf\it q}}_{1},{{\bf\it
q}}_{2}\right>={{\bf\it q}}_{2}^{*}{\bf\it W}{{\bf\it q}}_{1}$ that accounts
for the non-uniformly distributed nodes, the discrete LST and its adjoint
eigenvalue problems take the form of

$\displaystyle\lambda{\bf\it P}{\bf\it P}^{T}\tilde{{\bf\it q}}$
$\displaystyle={\bf\it L}\tilde{{\bf\it q}},$ (18a)
$\displaystyle\lambda^{+}{\bf\it P}{\bf\it P}^{T}\tilde{{\bf\it q}}^{+}$
$\displaystyle={\bf\it L}^{+}\tilde{{\bf\it q}}^{+},$ (18b)

respectively, where

$\displaystyle{\bf\it L}^{+}={\bf\it W}^{-1}{\bf\it L}^{H}{\bf\it W}$ (19)

is the discrete adjoint LNS operator. Here, the restriction matrix for two-
dimensional problems takes the form of ${\bf\it P}^{T}=\matrixquantity[{\bf\it
I}&{\bf\it 0}&{\bf\it 0}\\\ {\bf\it 0}&{\bf\it I}&{\bf\it 0}\\\ ]$, and
$\left(\cdot\right)^{H}$ denotes the Hermitian transpose. The eigenvectors of
these two generalized eigenvalue problems, $\tilde{{\bf\it q}}=[\tilde{{\bf\it
u}},\tilde{{\bf\it v}},\tilde{{\bf\it p}}]^{T}$ and $\tilde{{\bf\it
q}}^{+}=[\tilde{{\bf\it u}}^{+},\tilde{{\bf\it v}}^{+},\tilde{{\bf\it
p}}^{+}]^{T}$, are referred to as the LST and adjoint modes, respectively.

The wavemaker (WM), introduced by Giannetti and Luchini [46], identifies the
flow region with the strongest localized feedback, where the dominant
instability mechanisms act. Based on the leading LST and adjoint modes, the WM
is locally defined as

$\displaystyle\zeta_{\text{LST}}(x_{1},x_{2})=\frac{\|{\bf\it
P}^{T}{\bf\it\tilde{q}^{+}}(x_{1},x_{2})\|\,\|{\bf\it
P}^{T}{\bf\it\tilde{q}}(x_{1},x_{2})\|}{|\left<{\bf\it\tilde{q}^{+}},{\bf\it\tilde{q}}\right>|},$
(20)

where the norm $\|\cdot\|$ measures the localized energy. Referred to Luchini
and Bottaro [71] for a comprehensive review. In addition to quantifying the
receptivity, the wavemaker also indicates the non-normality level of the flow
field [27]. Initially introduced for base flows, Meliga _et al._ [82] extended
the application of wavemaker as a sensitivity analysis technique for mean
flows. We utilize this methodology to gauge the structural sensitivity of the
RBF-FD-based global LNS operator, ${\bf\it L}$.

###  3.3 Mesh-free resolvent analysis (RA)

The construction of the global Jacobian, ${\bf\it L}$, leads to the direct
discretization of the input-output system in equation (17), yielding

$\displaystyle{\bf\it\hat{u}}={\bf\it H}(\omega){\bf\it\hat{f}},$ (21)

where

$\displaystyle{\bf\it H}(\omega)={\bf\it P}^{T}{\bf\it
C}\left(\mathrm{i}\omega{\bf\it PP}^{T}-{\bf\it L}\right)^{-1}{\bf\it
P}{\bf\it B}$ (22)

is referred to as the discrete resolvent operator.

In the absence of nonlinear interactions described in table (1), the
application of input-output, or resolvent analysis (RA) provides a means to
model them as optimal forcing inputs to the linear system in equation (21).
The objective of resolvent analysis is to identify pairs of optimal forcings
and their corresponding responses that maximize the gain, $\sigma^{2}$,
defined as the ratio of the energy of the response to the energy of the
forcing,

$\displaystyle\sigma^{2}({\bf\it\hat{f}};\omega)=\frac{\|\hat{{\bf\it
u}}\|^{2}_{u}}{\|{\bf\it\hat{f}}\|^{2}_{f}}=\frac{\left<{\bf\it
H}(\omega){\bf\it\hat{f}},{\bf\it
H}(\omega){\bf\it\hat{f}}\right>_{u}}{\left<{\bf\it\hat{f}},{\bf\it\hat{f}}\right>_{f}}.$
(23)

Refer to Schmid and Henningson [115] for a detailed discussion. The energy of
the response and the forcing are measured in the norms $\|\cdot\|_{u}$ and
$\|\cdot\|_{f}$, induced by the inner products

$\displaystyle\left<\hat{{\bf\it u}}_{1},\hat{{\bf\it
u}}_{2}\right>_{u}=\hat{{\bf\it u}}_{2}^{*}{\bf\it W}_{u}\hat{{\bf\it
u}}_{1}\qquad\text{and}\qquad\left<{\bf\it\hat{f}}_{1},{\bf\it\hat{f}}_{2}\right>_{f}={\bf\it\hat{f}}_{2}^{*}{\bf\it
W}_{f}{\bf\it\hat{f}}_{1}$ (24)

on the output and input spaces, respectively. Here, ${\bf\it W}_{u}$ and
${\bf\it W}_{f}$ are weight matrices containing both the numerical quadrature
weights and weights associated with these inner products.

We further define the modified, or weighted, resolvent operator

$\displaystyle{\bf\it R}(\omega)\equiv{\bf\it W}_{u}^{\frac{1}{2}}{\bf\it
H}(\omega){\bf\it W}_{f}^{-\frac{1}{2}}=\hat{{\bf\it
U}}{\bf\it\Sigma}\hat{{\bf\it F}}^{*}$ (25)

to account for the energy in the input and output spaces. The optimal
responses ${\bf\it\hat{U}}=\matrixquantity[{\hat{{\bf\it
u}}}_{1},\cdots,{\hat{{\bf\it u}}}_{N}]$ and the corresponding forcings
${\bf\it\hat{F}}=\matrixquantity[{\bf\it\hat{f}}_{1},\cdots,{\bf\it\hat{f}}_{N}]$
are obtained from the singular value decomposition (SVD) of the modified
resolvent operator and ranked by the energy gains
$\Sigma=\matrixquantity[\sigma_{1},\cdots\sigma_{N}]$. The resulting modes are
orthogonal in their respective inner products, that is, $\left<\hat{{\bf\it
u}}_{j},\hat{{\bf\it
u}}_{k}\right>_{u}=\left<{\bf\it\hat{f}}_{j},{\bf\it\hat{f}}_{k}\right>_{f}=\delta_{jk}$.
The optimal input and output modes are related through

$\displaystyle{\bf\it
R}(\omega){\bf\it\hat{f}}_{j}=\sigma_{j}(\omega)\hat{{\bf\it u}}_{j},$ (26)

which provides a physical interpretation of the singular values and vectors.
In practice, the optimal input forcings, ${\bf\it\hat{f}}_{j}$, are determined
as the solutions of the eigenvalue problem

$\displaystyle{\bf\it W}_{f}^{-1}{\bf\it H}(\omega)^{*}{\bf\it W}_{u}{\bf\it
H}(\omega){\bf\it\hat{f}}_{j}=\sigma^{2}_{j}{\bf\it\hat{f}}_{j},$ (27)

where ${\bf\it W}_{f}^{-1}{\bf\it H}(\omega)^{*}{\bf\it W}_{u}{\bf\it
H}(\omega)$ is Hermitian. The matrix inversion of
$\left(\mathrm{i}\omega{\bf\it P}{\bf\it P}^{T}-{\bf\it L}\right)$ is solved
using LU-factorization [128].

The RA-based wavemaker,

$\displaystyle\zeta_{\text{RA}}(x_{1},x_{2})=\frac{\|{\bf\it\hat{u}}_{1}(x_{1},x_{2})\|\,\|{\bf\it\hat{f}}_{1}(x_{1},x_{2})\|}{|{\bf\it\hat{u}}_{1}^{*}{\bf\it
W}_{u}{\bf\it\hat{f}}_{1}|},$ (28)

is defined analogous to equation (20). Within the framework of RA, it provides
a quantitative measure of the effect of the localized feedback at any given
frequency [103, 109, 130, 135].

##  4 Numerical implementations

###  4.1 Node generation and parameter selection

We employ the unstructured triangular mesh generator DistMesh developed by
Persson [96] to efficiently generate scattered nodes with localized refinement
in regions of interest. We highlight that the local connection information is
not used in the computation. The accuracy and stability of the PHS+poly RBF-FD
method are contingent upon three parameters: the stencil size, $n$, the PHS
exponent, $m$, and the polynomial degree, $q$. we investigate the order of
accuracy ranging from the standard choice of $q=2$ to the desired high-order
accuracy for hydrodynamic stability analysis, here $q=4$. Increasing the order
of accuracy generally leads to a decrease in the size of the linear operator.

Figure 1: Maximum local $\mathcal{L}$-Lebesgue function for different
combinations of PHS exponents, $m$, and polynomial orders, $q$ for the first
(a-c) and second derivatives (d-f) in the $x$\- (solid line), $y$\- (’+’) and
$xy$-directions (circle).

The local $\mathcal{D}$-Lebesgue function proposed by Shankar and Fogelson
[123] provides a measure of the eigenvalue stability of the RBF-FD method. For
a given linear operator $\mathcal{D}$, the local Lebesgue function takes the
form of the 1-norm of the RBF-FD weight in equation (1), that is

$\displaystyle\Lambda_{\mathcal{D}}({\bf\it x}_{0};\\{{\bf\it
x}_{j}\\}_{j=1}^{n})=\|{\bf\it w}\|_{1}=\sum_{j=1}^{n}|w_{j}|.$ (29)

Larger values of the local Lebesgue function indicate an increased
susceptibility to numerical instability in the assembled global
differentiation matrix. Figure 1 shows the maximum value of the local Lebesgue
function across various parameter combinations based on the test grid shown in
figure 4. Refer to the accompanying context of figure 4 for more details of
the grid. For all considered stencil sizes, polynomial degrees, and relevant
operators, the minimax value is consistently achieved at $m=3$, suggesting its
suitability for hydrodynamic stability analysis.

For two-dimensional problems, it has been suggested that the stencil size
should satisfy $n\gtrsim(q+1)(q+2)$ [39, 40]. In practice, including a few
additional nodes beyond the minimum requirement is often beneficial to improve
performance. The specific choice of additional nodes may depend on the type of
node sets used. For instance, a recommended formula for Halton nodes is
$n=(q+1)(q+2)+\lfloor{\ln{\left[(q+1)(q+2)\right]}}\rfloor$ [67, 123]. We here
select the stencil size $n$ based on the assumption of a perfectly arranged
hexagonal node distribution comprising $q/2+1$ layers, despite the actual
distribution of nodes being heterogeneous. For practical examples, refer to
figure 3 and the surrounding context. The recommended parameter combinations
are summarized in table 2. In particular, we use $(n,m,q)=(37,3,4)$ for
interior nodes and $(n,m,q)=(19,3,2)$ for nodes near boundaries in the
construction of the discrete LNS operator, ${\bf\it L}$. The RBF stencil for
each node is determined as the nearest $n$ nodes acquired through a k-nearest
neighbor (kNN) search.

| Shankar and Fogelson [123] | Le Borne and Leinen [67] | present
---|---|---|---
$q=2$ | $(n,m)=(14,3)$ | $(n,m)=(14,5)$ | $(n,m)=(19,3)$
$q=3$ | $(n,m)=(22,3)$ | $(n,m)=(22,7)$ | $(n,m)=(28,3)$
$q=4$ | $(n,m)=(33,3)$ | $(n,m)=(33,7)$ | $(n,m)=(37,3)$
Table 2: Summary of parameter selections. Figure 2: Maximum condition number
of ${\bf\it A}_{\text{aug}}$ as a function of local grid spacing
$\Delta\tilde{r}$: (a) $n=19,\,q=2$; (b) $n=28,\,q=3$; (c) $n=37,\,q=4$. The
recommended local spacing of $\Delta\tilde{r}=0.25$ is highlighted as dashed
lines.

The condition number of the augmented matrix ${\bf\it A}_{\text{aug}}$,
defined in equation (7), affects the numerical stability in a similar way and
hence needs to be considered separately. The condition number for various
parameter combinations has been investigated and discussed in Le Borne and
Leinen [67]. Independent of the parameter selection, Shahane _et al._ [119]
demonstrated that scaling the local stencil to a unit length scale can improve
the condition number. We further reveal the existence of an optimal range for
the averaged local grid spacing, wherein the condition number of the augmented
matrix ${\bf\it A}_{\text{aug}}$ reaches its minimum value.

For a given stencil $\\{{{\bf\it x}}_{j}\\}_{j=1}^{n}$, we perform scaling
with respect to the location of ${\bf\it x}_{0}$ with a scaling factor
$\kappa$ such that

$\displaystyle\\{\tilde{{\bf\it x}}_{j}\\}_{j=1}^{n}=\\{\kappa({{\bf\it
x}}_{j}-{\bf\it x}_{0})\\}_{j=1}^{n},$ (34)

and the averaged local grid spacing becomes $\Delta\tilde{r}=\kappa\Delta r$.
Figure 2 shows the maximum condition number within the test grid as a function
of $\Delta\tilde{r}$ for the considered parameters. Optimal performance is
consistently achieved within the range of $0.15\lesssim\Delta\tilde{r}\lesssim
0.3$ across all examined parameter combinations. Based on this, we perform
spatial scaling to ensure an averaged local spacing of $\Delta\tilde{r}=0.25$
before solving the RBF-FD weights in equation (1). The chain rule is applied
to transform the RBF weights back to the original grid.

###  4.2 Boundary condition treatments

####  4.2.1 Homogeneous Dirichlet boundary conditions

Homogeneous Dirichlet boundary conditions, ${\bf\it q}^{\prime}=0$, are widely
employed in hydrodynamic stability analysis based on the assumption that
perturbation variables either vanish at the solid wall or in the far field. To
approximate spatial derivatives for a node near such a boundary, we follow the
Kansa method [65] and compute the RBF-FD weights using the local stencil
consisting of both the interior nodes, $\\{{\bf\it
x}^{(i)}\\}_{j=1}^{n^{(i)}}$, and the boundary nodes, $\\{{\bf\it
x}^{(b)}\\}_{j=1}^{n^{(b)}}$. Letting $g({\bf\it x}_{j}^{(b)})=0$, equation
(1) becomes

$\displaystyle\mathcal{D}g({\bf\it x}_{0})$
$\displaystyle=\sum_{j=1}^{n^{(i)}}w_{j}^{(i)}g({\bf\it
x}_{j}^{(i)})+\sum_{j=1}^{n^{(b)}}w_{j}^{(b)}g({\bf\it
x}_{j}^{(b)})=\sum_{j=1}^{n^{(i)}}w_{j}^{(i)}g({\bf\it x}_{j}^{(i)}),$ (35)

and only the weights for the interior nodes will be used in the computation.

####  4.2.2 Boundary-normal derivatives

Neumann conditions, inflow/outflow conditions, or the enforcement of
continuity necessitate the evaluation of boundary-normal derivatives and
require special treatment for scattered nodes that are not structured near
boundaries. To accurately approximate derivatives in the boundary-normal
direction, we propose an elliptical stencil with its major axis perpendicular
to the boundary. The elliptical form is achieved by utilizing a weighted
Euclidean distance metric

$\displaystyle\tilde{r}=\sqrt{a^{2}+b^{2}}\sqrt{\frac{\left(({\bf\it
x}-{\bf\it x}_{0})\cdot{\bf\it n}\right)^{2}}{a^{2}}+\frac{\left(({\bf\it
x}-{\bf\it x}_{0})\cdot{\bf\it t}\right)^{2}}{b^{2}}},$ (36)

where ${\bf\it n}$ and ${\bf\it t}$ are the normal and tangential directions,
respectively, and $a$ and $b$ denote the corresponding scaling factors. An
eccentricity of $\frac{\sqrt{a^{2}-b^{2}}}{a}={\frac{2\sqrt{6}}{5}}$ is used
in the computation. We use the $q=2$ polynomial augmentation for the boundary-
tangential direction and a higher-order of accuracy of $q=4$ for the boundary-
normal direction. An illustrative example is shown in figure 3(a). This new
strategy circumvents the need to increase the local stencil near domain
boundaries to prevent the Runge phenomenon, as demonstrated in Bayona _et al._
[15], Bayona [13]. Without this special treatment, stable solutions for
hydrodynamic stability analysis cannot be attained.

####  4.2.3 Symmetric and anti-symmetric boundary conditions

Figure 3:  RBF stencils (red-shaded circles) for a given node (blue star): (a)
boundary-normal derivatives; (b) interior node near the symmetric centerline
(dot-dashed). The latter stencil includes the interior nodes ($\\{{\bf\it
x}^{\bullet}\\}_{j=1}^{n^{\bullet}}$, dot) and image ghost nodes ($\\{{\bf\it
x}^{+}\\}_{j=1}^{n^{+}}$, ’+’). The counterparts of the image nodes,
$\\{{\bf\it x}^{\Delta}\\}_{j=1}^{n^{+}}$, and the remaining nodes,
$\\{{\bf\it x}_{j}^{\circ}\\}_{j=1}^{n^{\bullet}-n^{+}}$, are shown as
triangles and circles, respectively. The area of each red-shaded circle
represents the corresponding local radial control volume, $\differential
V_{i}$.

To address pole singularities that arise in the polar coordinates and
symmetric or anti-symmetric boundary conditions, we propose a generalization
of the pole treatment method by Mohseni and Colonius [86] for scattered nodes.
We introduce a set of ghost nodes, denoted as $\mathcal{X}^{+}$, which are
symmetric to the interior nodes, $\mathcal{X}^{\bullet}$, with respect to the
centerline. For each interior node, we divide its local stencil into two
disjoint sets such that $\\{{\bf\it x}\\}_{j=1}^{n}=\\{{\bf\it
x}_{j}^{\bullet}\\}_{j=1}^{n^{\bullet}}\cup\\{{\bf\it
x}_{j}^{+}\\}_{j=1}^{n^{+}}$, where ${\bf\it
x}^{\bullet}\in\mathcal{X}^{\bullet}$ and ${\bf\it x}^{+}\in\mathcal{X}^{+}$,
respectively, and $n=n^{\bullet}+n^{+}$. Additionally, we define $\\{{\bf\it
x}_{j}^{\Delta}\\}_{j=1}^{n^{+}}\in\\{{\bf\it
x}_{j}^{\bullet}\\}_{j=1}^{n^{\bullet}}$ as the counterparts of the image
nodes, and $\\{{\bf\it x}_{j}^{\circ}\\}_{j=1}^{n^{\bullet}-n^{+}}=\\{{\bf\it
x}_{j}^{\bullet}\\}_{j=1}^{n^{\bullet}}\setminus\\{{\bf\it
x}_{j}^{\Delta}\\}_{j=1}^{n^{+}}$ represents the remaining nodes within the
interior nodes, Refer to figure 3(b) for detailed symbol explanations. The
function values at the image nodes are determined by their corresponding
counterparts, given by

$\displaystyle g({\bf\it x}_{j}^{+})=\eta g({\bf\it x}_{j}^{\Delta}),$ (37)

where $\eta$ depends on the type of boundary conditions being imposed.
Specifically, we use $\eta=1$ for symmetric boundary conditions and $\eta=-1$
for anti-symmetric boundary conditions. The RBF-FD weights in equation (1) can
be written as

$\displaystyle\mathcal{D}g({\bf\it x}_{0})$
$\displaystyle=\sum_{j=1}^{n^{\bullet}}w_{j}^{\bullet}g({\bf\it
x}_{j}^{\bullet})+\sum_{j=1}^{n^{+}}w_{j}^{+}g({\bf\it
x}_{j}^{+})=\sum_{j=1}^{n^{\bullet}-n^{+}}w_{j}^{\circ}g({\bf\it
x}_{j}^{\circ})+\sum_{j=1}^{n^{+}}\left(\eta
w_{j}^{+}+w_{j}^{\Delta}\right)g({\bf\it x}_{j}^{\Delta}).$ (38)

This treatment allows us to approximate derivatives solely based on the
function values at the interior nodes, effectively handling pole singularities
and addressing the challenges associated with scattered nodes.

##  5 Applications

Flow | $\mathrm{Re}$ | $M$ | Flow type | Base state | Analysis | Sec.
---|---|---|---|---|---|---
Cylinder wake [29] | $47-180$ | - | 2D laminar | $\overline{u},\overline{v}$ | LST/RA/WM | §5.1
Boundary layer (ZPG) | $6\times 10^{5}$ | - | 2D laminar | $u_{\text{ZPG}},v_{\text{ZPG}}$ | RA/WM | §5.2
Jet [23] | $\approx 10^{6}$ | 0.9 | 3D turbulent | $\overline{\rho},\overline{u}_{x},\overline{u}_{r},\overline{u}_{\theta},\overline{T}$ | RA/WM | §5.3
Table 3: Overview of datasets and analyses. The columns from left to right
indicate the flow description, Reynolds number, Mach number, flow type, base
state, analysis type, and section number. The zero-pressure-gradient (ZPG)
Blasius solution is used for analyzing the boundary layer. Analyses include
linear stability (LST) and resolvent analyses (RA), along with wavemaker (WM).

We demonstrate the mesh-free RBF-FD-based hydrodynamic stability framework
outlined in §3 using three representative examples: canonical steady and
unsteady cylinder wakes, a self-similar non-parallel steady laminar boundary-
layer flow, and the turbulent mean of a transonic jet, as summarized in table
3. These three examples are benchmark problems for open flows and are
appropriate for validating mesh-free hydrodynamic stability.

###  5.1 Cylinder wake

We first consider the incompressible cylinder flow at diameter-based Reynolds
numbers, $\mathrm{Re}=\frac{U_{\infty}D}{\nu}$, ranging from 47 to 180 and
investigate the mean-flow stability within the two-dimensional laminar regime.
The occurrence of the periodic von Kármán vortex shedding in the cylinder wake
beyond the critical Reynolds number of $\mathrm{Re_{c}}\simeq 47$ is a well-
known phenomenon, owing to a Hopf bifurcation that results in flow
instability, see, e.g., [17, 148]. Beyond the limit of $\mathrm{Re}\simeq
188$, the cylinder flow becomes three-dimensional [12, 150].

Classical LST analysis of the cylinder base flow ought to predict the onset of
unsteadiness [57, 154] but fails to capture the vortex-shedding frequency
beyond $\mathrm{Re}_{c}$ [11, 127]. Previous studies by Hammond and Redekopp
[49] and Pier [100] show that LST analysis around the cylinder mean flow
accurately identifies the vortex-shedding frequency compared to experimental
measurements. Barkley [11] supported these findings and showed that the
cylinder mean flow is marginally stable in the 2D regime. Sipp and Marquet
[127] subsequently provided theoretical underpinning by conducting a weakly
nonlinear analysis and establishing criteria for utilizing mean flows in LST
analysis. The vortex-shedding dynamics beyond the critical Reynolds number
were later investigated using a self-consistent model [74, 75] and RA [61,
135] based on the mean flow stability. We here conduct both LST and resolvent
analyses of the cylinder wake mean flow.

Figure 4: Computational grid for cylinder flows with $N=118225$ nodes, colored
using the mean vorticity, $\overline{{\bf\it\omega}}={\bf\it
D}_{x}\overline{{\bf\it v}}-{\bf\it D}_{y}\overline{{\bf\it u}}$, at
$\mathrm{Re}=100$.

We define the the computational domain $\Omega$ as the exterior of the
cylinder $r\geq D/2=0.5$ and within the rectangle $-15\leq x\leq 30,\,-15\leq
y\leq 15$. The computational domain is discretized using $N=118225$ scattered
nodes. Local grid refinement is employed near the cylinder with a
characteristic distance of $\Delta r=0.03$ and around the wake centerline with
$\Delta r=0.04$ to better resolve the flow structures. The unsteady cylinder
flow is simulated using the PHS+poly RBF-FD version of the fractional-step,
staggered-grid incompressible Navier-Stokes solver by Chu and Schmidt [29].
The mean-flow profiles are obtained as the time average of flow over 20
vortex-shedding cycles. Figure 4 shows the mean vorticity at $\mathrm{Re}=100$
and the computational grid. Homogeneous boundary conditions,
$u^{\prime}=v^{\prime}=0$, are prescribed at the inlet and the cylinder
surface. Symmetric boundary conditions with $v^{\prime}=\partial
u^{\prime}/\partial y=0$ are applied at the transverse boundaries. A stress-
free outflow condition, $-p^{\prime}{\bf\it
n}+\frac{1}{\mathrm{Re}}{\nabla}{\bf\it u}^{\prime}\cdot{\bf\it n}=0$, where
${\bf\it n}=[1,0]^{T}$ is the outflow direction, is enforced at the outflow.

The local wavemaker sensitivity in equation (20) and the resolvent gain in
equation (23) and are both quantified in terms of the perturbed kinetic
energy. To this end, we define the integration matrices

in equation (24) to approximate the integral within the computational domain,
where $\otimes$ denotes the Kronecker product. Here, $\differential
V_{i}=\xi\pi(\Delta r_{i})^{2}$ is the local radial control volume for each
grid. The constant $\xi$ ensures consistency with the total control volume and
is defined by letting $\sum_{i=1}^{N}\differential
V_{i}=\xi\pi\sum_{i=1}^{N}(\Delta
r_{i})^{2}=\int_{\Omega}1\differential{\bf\it x}$. Refer to figure 3 for
practical illustrations.

Figure 5: Vortex shedding frequency predicted by the leading eigenvalue of the
mean-flow stability problem as a function of Reynolds number. The frequency is
given as the non-dimensional Strouhal number $\mathrm{St}=\lambda_{i}/2\pi$.
Shown for comparison are results from Pier [100] (green square) and Barkley
[11] (blue circle). Two representative Reynolds numbers for the following
analysis, $\mathrm{Re}_{c}\approx 47$ and $\mathrm{Re}=100$, are highlighted
as dashed lines.

Figure 5 shows the Strouhal number, $\mathrm{St}=\lambda_{i}/2\pi$, associated
with the leading eigenvalues at varying Reynolds numbers. Starting from the
critical Reynolds number of $\mathrm{Re}_{c}\approx 47$, the frequency-
Reynolds number dependence exhibits the typical features of a Hopf
bifurcation. Our results are in good agreement with Barkley [11], and
similarly, deviate no more than $3\%$ from those reported by Pier [100]. The
maximum growth rates are almost identical to zero (see also figure 6 below),
confirming that the mean flow is marginally stable.

Figure 6: Resolvent singular values (blue curve) and stability eigenvalue (red
circle) spectra for $\mathrm{Re}=47$ (a) and $\mathrm{Re}=100$ (b). The 20
eigenvalues closest to $\lambda_{c}=0+0.753\mathrm{i}$ or
$\mathrm{St}_{c}=0.1199$ and $\lambda=0+1.038\mathrm{i}$ or
$\mathrm{St}=0.1652$ (dashed lines) were found within the regions outlined by
the black dot-dashed lines using the shift-and-invert Arnoldi algorithm for
$\mathrm{Re}=47$ and $100$, respectively.

We next conduct a comparative study of LST and RA at two representative
Reynolds numbers, the critical Reynolds number of $\mathrm{Re}=47$ and
$\mathrm{Re}=100$, as an example of the unsteady regime. Figure 6 shows the
resulting resolvent singular value and stability eigenvalue spectra. At
$\mathrm{Re}=\mathrm{Re}_{c}$, both the peak of the resolvent gain and the
leading eigenvalue identify the same frequency, $\mathrm{St}_{c}=0.1199$, as
the vortex shedding frequency. This value is in good agreement with the
results of the LST analysis around the base flow, see, e.g., Giannetti and
Luchini [46] ($\mathrm{St}_{c}\simeq 0.118$ for $\mathrm{Re}_{c}\simeq 46.7$),
Marquet _et al._ [76] ($\mathrm{St}_{c}\simeq 0.116$ for
$\mathrm{Re}_{c}\simeq 46.8$), and Sipp and Marquet [127]
($\mathrm{St}_{c}\simeq 0.118$ for $\mathrm{Re}_{c}\simeq 46.6$). Similarly,
the resolvent singular value spectrum at $\mathrm{Re}=100$ in panel 6(b)
displays a clear peak at the vortex shedding frequency, which is now at
$\mathrm{St}=0.1652$, and again coinciding with the least stable LST
eigenvalue. This result matches closely with the findings of earlier Strouhal-
Reynolds number relationship [38, 60, 74, 101, 150] and RA [61, 135] studies.

Figure 7: Leading modes for $\mathrm{Re}=47$ at $\mathrm{St}=0.1199$ (a-d) and
$\mathrm{Re}=100$ at $\mathrm{St}=0.1652$ (e-h): (a,e) LST modes; (b,f)
response modes; (c,g) adjoint modes; (d,h) forcing modes.

Figure 7 shows the leading LST and RA modes for the cylinder mean flow at
$\mathrm{Re}=\mathrm{Re}_{c}$ and $\mathrm{Re}=100$. The leading LST and
optimal response modes, along with their corresponding adjoint and forcing
modes, are near-identical at both Reynolds numbers. The resemblance observed
between the LST and resolvent response modes is to be anticipated, as the
singular value of the resolvent attains its peak at the vortex-shedding
frequency [18, 135], thus exhibiting the characteristic vortex-shedding
structure. This implies that the optimal forcing leverages the global
instability mode to achieve maximum gain. Differing from the optimal response
and LST modes, the optimal forcing and adjoint modes are active far upstream
of the cylinder but peak downstream in close vicinity to the cylinder. This
hallmark of convective instability was similarly observed in previous works
[61, 76, 135].

Figure 8:  Wavemakers $\zeta_{\text{LST}}$ for mean flow at (a)
$\mathrm{Re}=47$ and (c) $\mathrm{Re}=100$. Respective results for
$\zeta_{\text{RA}}$ are shown in (b) and (d). The blue curve represents the
mean flow streamlines.

We finally investigate the wavemaker $\zeta_{\text{LST}}$, defined in equation
(20), in figure 8 to quantify the sensitivity of spatially localized feedback.
Wavemakers obtained from LST and RA look similar for both Reynolds numbers and
reach their maxima in two symmetrically positioned lobes located across the
separation bubble. This result signifies the promising applicability of RA-
based wavemakers in accurately identifying the region where the flow
instability mechanisms happen. The wavemaker patterns at the critical Reynolds
number shown in panels 8 (a,c) compare well to those previously reported by
Giannetti and Luchini [46] and Marquet _et al._ [76]. The shrinking of
wavemaker region at $\mathrm{Re}=100$ in panels 8 (b,d) closely matches with
the findings in Meliga _et al._ [82] and Symon _et al._ [135].

###  5.2 Blasius boundary layer

We now examine the incompressible two-dimensional non-normal flat-plate
boundary layer, a classic example of convectively unstable flows characterized
by the amplification of disturbances during downstream advection. The
convective instability of boundary-layer flows has been extensively studied
through 1D LST analysis over the past century, e.g., [63, 72, 73, 106, 132].
Above the critical displacement-thickness Reynolds number of
$\mathrm{Re}_{\delta,c}\approx 520$, the Tollmien–Schlichting (TS) waves are
known to arise as unstable eigenmodes of the Orr-Sommerfeld equation [25].
While these locally unstable waves are damped in 2D LST [1, 2, 7, 35],
investigations into the non-normality of the linearized Navier-Stokes
equations for open flows have revealed an alternative pathway for disturbance
amplification [27, 77, 129]. To analyze the non-modal behavior of boundary-
layer flows, input-output analysis has been conducted to determine the optimal
harmonic forcing that results in the largest asymptotic response [3, 8, 9, 89,
115, 128].

Following the setup by Sipp and Marquet [128], we compute the resolvent gain
in the restricted domain
$\Omega_{\text{RA}}=[x,\,y/\delta]\in[0.02,\,1]\times[0,\,22.52]$ within the
computational domain $\Omega=[x,\,y/\delta]\in[0.02,\,1.27]\times[0,\,22.52]$
such that the forcing optimizes the ratio between the restricted kinetic
energy,
$\iint_{\Omega_{\text{RA}}}\left(|\hat{u}|^{2}+|\hat{v}|^{2}\right)\,\differential
x\differential y$, and the integral of forcing,
$\iint_{\Omega}\left(|\hat{f}_{u}|^{2}+|\hat{f}_{v}|^{2}\right)\,\differential
x\differential y$. The zero-pressure-gradient (ZPG) asymptotic Blasius
solution, characterized by a local Reynolds number of
$\mathrm{Re}_{x}=\frac{U_{\infty}x}{\nu}=6\times 10^{5}$ or
$\mathrm{Re}_{\delta}=\frac{U_{\infty}\delta(x)}{\nu}=1332$ at the end of the
restricted domain $\Omega_{\text{RA}}$, is used as the base flow to
investigate the flow instabilities, where
$\delta(x)=1.72\sqrt{x/\mathrm{Re}_{x}}$ is the displacement thickness. The
leading edge ($x<0.02$) is removed to avoid the singularity of the self-
similar solution. Hereafter, we use the notations
$\mathrm{Re}=\mathrm{Re}_{x=1}$ and $\delta=\delta(x=1)$ for simplicity.

The computational domain $\Omega$ is discretized using $N=726143$ scattered
nodes, corresponding to $2.2\times 10^{6}$ degrees of freedom. In comparison,
Sipp and Marquet [128] employ $13.7\times 10^{6}$ degrees of freedom. The
characteristic distances of the grid are $\Delta r=0.069\delta$ near the flat
plate ($y/\delta<4$), $\Delta r=0.077\delta$ near the inlet ($x<0.025$), and
average at $0.081\delta$ over the whole domain. At the inlet and the flat
plate, we prescribe $u^{\prime}=v^{\prime}=0$. A symmetric boundary condition
with $v^{\prime}=\partial u^{\prime}/\partial y=0$ is applied at the far-field
($y/\delta=22.52$). A stress-free outflow condition is enforced at the
outflow. The weight matrix ${\bf\it W}_{u}$ is defined in equation (5.1), but
with zero weights in the region $\Omega\setminus\Omega_{\text{RA}}$.

Figure 9:  Resolvent gains (solid lines) and peak frequencies (dashed lines)
for the flat-plate boundary layer and $\mathrm{Re}=6\cdot 10^{5}$.

We first examine the leading resolvent singular value as a function of the
normalized frequency, $F=10^{6}\cdot\omega/\mathrm{Re}$, in figure 9. Overall,
our results agree well with those reported by Sipp and Marquet [128]. The
slight deviation of the peak and the resolvent gains beyond the peak
($F\gtrsim 88$) are most likely attributed to differences in the base flows,
that is, the ZPG self-similar Blasius solution in the present study and the
fully non-parallel numerical data of Sipp and Marquet [128] that includes the
leading edge. The discrepancy at higher frequencies is potentially also
related to the truncation of the leading edge, where the asymptotic Blasius
solution becomes singular.

Figure 10: Optimal and suboptimal resolvent forcings (b,d) and corresponding
responses (a,c) for the flat-plate boundary layer at
$F=10^{6}\cdot\omega/\mathrm{Re}=100$. The normalized stream-wise velocity
components have been interpolated onto a stretched Cartesian mesh for
visualization. Panels (e) and (f) show the local regions of the optimal
forcing (magenta box) and corresponding response modes (red box) with the
largest magnitudes, respectively, on the scattered nodes used for the
computation. The Blasius displacement thickness, $\delta(x)$, is highlighted
as the white dashed line.

Figure 10 shows the optimal and suboptimal forcings and corresponding
responses at $F=100$. The optimal response exhibits clear Tollmien–Schlichting
(TS) wavepackets in the downstream region, while the upstream tilted
structures in the optimal forcing highlight the active role of the Orr
mechanism in extracting energy from the mean shear via the Reynolds stress
[25]. The clear spatial separation between leading resolvent forcing and
response modes indicates the stream-wise non-normality of the system [27, 76,
129]. The suboptimal forcing and response are similar to the leading modes
except for two local maxima. This modulation confirms the orthogonality in
their respective inner products. Qualitative comparisons with previous studies
by Åkervik _et al._ [1], Monokrousos _et al._ [89], Brandt _et al._ [22] and
Sipp and Marquet [128] verify the capability of the present framework in
identifying the convective instability of the boundary-layer flow.

Figure 11: Spatial distributions of energy density for optimal (dashed) and
suboptimal (dotted) forcings and corresponding responses (solid and dot-
dashed, respectively) at $F=100$. The results (red) are compared to those
reported by Sipp and Marquet [128] (blue). The two vertical solid lines
represent the upstream neutral point (branch I) and the downstream neutral
point (branch II) from a local stability analysis.

As a quantitative assessment of the flow structures, we examine the energy
density functions,

$\displaystyle
d_{f}(x)=\int_{0}^{y_{\text{max}}}\left(|\hat{f}_{u}|^{2}+|\hat{f}_{v}|^{2}\right)\,\differential
y,\qquad\text{and}\qquad
d_{u}(x)=\int_{0}^{y_{\text{max}}}\left(|\hat{u}|^{2}+|\hat{v}|^{2}\right)\,\differential
y,$ (43)

in figure 11 for the modes shown in figure 10. Our results are almost
identical to Sipp and Marquet [128] The spatial distribution of the optimal
forcing unambiguously identifies the location of the upstream neutral point
(branch I) from a local stability analysis at $x=0.3$ and the corresponding
response is localized at $x=0.89$, which is in close proximity to the
downstream neutral point (branch II).

Figure 12: Energy density distributions for the optimal forcings (a) and
responses (b) as a function of frequency. The locations of maximum energy
densities are marked as red circles. The results reported by Sipp and Marquet
[128] are shown as blue stars for comparison. The neutral curve obtained from
a local LST analysis is shown as the magenta line.

Figure 12 shows the energy density distributions of the optimal forcing and
response as a function of frequency, with the stream-wise coordinate given in
terms of the local Reynolds number, $\mathrm{Re}_{\delta}$. The optimal
forcings agree well with Sipp and Marquet [128], and their maxima effectively
delineate the convectively stable/unstable boundary (branch I) obtained using
local stability theory. The Orr and TS mechanisms coexist and compete while
both contribute to the overall energy gain of RA. As frequency increases, the
spatial support of the TS-like optimal responses decreases, suggesting that
the TS mechanism is only supported in a limited region of the shear layer at
high frequencies. Additionally, the energy density of the forcing shows an
increasing trend, indicating that the Orr mechanism becomes dominant at high
frequencies. Similar to the deviation in the resolvent gain spectra shown in
figure 9, we observe a slight downstream shift in the peak of optimal
responses compared to the results reported Sipp and Marquet [128] for
increasing frequency ($F\gtrsim 130$).

Figure 13: RA-based wavemaker $\zeta_{\text{RA}}$ for the flat-plate boundary
layer as a function of frequency.

Finally, we show the RA-based wavemaker, defined in equation (28), in figure
13 to examine the instability mechanisms at different frequencies. The maximum
magnitude of $\zeta_{\text{RA}}$ shows an increasing trend with frequency. The
large magnitudes of the wavemaker are attributed to the fact that the value of
the term $|{\bf\it\hat{u}}_{1}^{*}{\bf\it W}_{u}{\bf\it\hat{f}}_{1}|^{-1}$ is
considerably greater than the value of the self-adjoint modes, which is equal
to 1 [77]. This again confirms that the boundary-layer flow exhibits high non-
normality. For $F\lesssim 100$ (panels 13 (a-c)), the RA-based wavemaker
exhibits an elongated shape and is artificially restricted within the
optimization domain $\Omega_{\text{RA}}$. The wavemaker becomes more
concentrated towards the upstream region for higher frequencies ($F\gtrsim
140$, panels 13 (d-g)). This observation suggests a transition from the TS-
dominated to the Orr-dominated mechanism as the frequency increases.

###  5.3 Turbulent jet

The above two examples demonstrate the capability of the proposed numerical
framework for analyzing incompressible flows within the laminar regime. We now
focus on the mean flow analysis of a turbulent iso-thermal jet at Mach number,
based on the jet velocity and the far-field speed of sound, of $M=0.9$ and
Reynolds number, based on the nozzle diameter and the jet velocity, of
$\real\approx 10^{6}$. The large eddy simulation (LES) data, generated using
the unstructured flow solver ‘Charles’, is used for analysis. Further details
about the dataset can be found in Brès _et al._ [23]. Previous studies have
demonstrated that the transonic turbulent jet under consideration displays a
variety of coherent features, including the well-known Kelvin–Helmholtz
instabilities of the shear layer [23], downstream non-modal Orr-type waves
[99, 117], and trapped acoustic waves in the potential core [116, 141].

The state vector,

$\displaystyle{\bf\it q}=[\rho,u_{x},u_{r},u_{\theta},T]^{T},$ (44)

comprises primitive variables: density $\rho$, temperature $T$, and
cylindrical velocity components $u_{x}$, $u_{r}$ and $u_{\theta}$ in the
streamwise, $x$, radial, $r$, and circumferential, $\theta$, directions,
respectively. The fluctuating state is defined as
${\bf\it{q^{\prime}}}=[{\rho^{\prime}},{u}_{x}^{\prime},{u}_{r}^{\prime},u_{\theta}^{\prime},{T^{\prime}}]^{T}$.
Owing to the rotational symmetry of the jet, we may decouple the governing
equations and construct the compressible linear operator $\mathcal{L}$ in
cylindrical coordinates, without loss of generality, for each azimuthal
wavenumber $m_{\theta}$ independently. For a round jet, the mean azimuthal
velocity component is zero. Upon linearization of the compressible Navier-
Stokes equations shown in (A.1), we obtain the linearized Navier-Stokes
operator around the azimuthally averaged long-time mean of the primitive
state, $\overline{{\bf\it
q}}=[\overline{\rho},\overline{u}_{x},\overline{u}_{r},0,\overline{T}]^{T}$,
such that

$\displaystyle\partialderivative{t}{\bf\it{q^{\prime}}}=\mathcal{L}{\bf\it{q^{\prime}}}+{\bf\it
f}.$ (45)

The general setup, including boundary conditions, sponge regions, and a
molecular Reynolds number of $\mathrm{Re}=3\times 10^{4}$, follows Schmidt _et
al._ [116] and Schmidt _et al._ [117].

Figure 14: Computational grid for jet with $N=210817\approx 1013\times 208$
nodes, colored using the mean streamwise velocity $\overline{u}_{x}$ at
$\mathrm{Re}\approx 10^{6}$. The potential core (white solid) and the jet
width (white dashed) are delineated as isolines corresponding to $99\%$ and
$5\%$ of the jet velocity, respectively

The computational domain $\Omega=[-1,31.01]\times[0,9.65]$ includes the
physical domain $x,r\in[0,30]\times[0,6]$ and the surrounding sponge regions
that prevent waves from being reflected. Local refinement is used near the
nozzle with $\Delta r=0.005$ and near the jet centerline with $\Delta
r=0.008$, resulting in $N=210817$ scattered nodes for the construction of the
RBF-FD-based global Jacobian ${\bf\it L}$, see figure 14. A structured mesh
with a comparable number of $N=185250$ was used in Schmidt _et al._ [117].

The resolvent gain is quantified in terms of the compressible energy norm [28]
through the weight matrices

$\displaystyle{\bf\it W}_{q}={\bf\it
W}_{f}\equiv\mathrm{diag}\left(\frac{\overline{T}}{\gamma\overline{\rho}M^{2}},\overline{\rho},\overline{\rho},\overline{\rho},\frac{\overline{\rho}}{\gamma(\gamma-1)\overline{T}M^{2}}\right)\otimes\mathrm{diag}(2\pi
r_{1}\differential V_{1},2\pi r_{2}\differential V_{2},\cdots,2\pi
r_{N}\differential V_{N}).$ (46)

For compressible flows, the discrete weighted resolvent operator takes the
form of

$\displaystyle{\bf\it R}(\omega)={\bf\it W}_{q}^{\frac{1}{2}}{\bf\it
C}\left(\mathrm{i}\omega{\bf\it I}-{\bf\it L}\right)^{-1}{\bf\it B}{\bf\it
W}_{f}^{-\frac{1}{2}},$ (47)

where the input and output matrices, ${\bf\it B}$ and ${\bf\it C}$, are used
to are used to focus the analysis exclusively on the physical domain. As a
first demonstration, we conduct the mesh-free RA for the symmetric component
of the jet with $m_{\theta}=0$ in the following.

Figure 15:  Leading resolvent singular value spectra for the transonic jet and
$m_{\theta}=0$. Results reported by Schmidt _et al._ [117] are shown as
comparisons (blue dashed).

Figure 15 compares the leading resolvent singular value spectra to those
reported by Schmidt _et al._ [117], obtained using a 4th-order finite
difference [80] and a tenth-order filter for discretization. Very good
agreement is observed within the frequency range
$0.15\lesssim\mathrm{St}\lesssim 1.2$, with only a minor deviation at the peak
value. This specific frequency range has been previously identified as the
regime where different physical mechanisms are active in turbulent jets, see,
e.g., [44, 98, 116, 117, 134, 138, 141].

Figure 16:  Streamwise velocity component of the optimal response (a,c,e) and
corresponding forcing modes (b,d,f) at three representative frequencies: (a-b)
$\mathrm{St}=0.2$; (c-d) $\mathrm{St}=0.6$; (e-f) $\mathrm{St}=1$. The modes
are interpolated onto a stretched Cartesian mesh with contours corresponding
to $\pm 0.6\|\cdot\|_{\infty}$ for visualization. The potential core and the
jet width, shown in figure 14, are visualized for comparisons.

To verify that large-scale coherent structures in the turbulent jet are
accurately captured, we investigate the leading resolvent modes at three
representative frequencies, as shown in figure 16. At $\mathrm{St}=0.2$, the
optimal response exhibits a clear downstream ($x\gtrsim 10$) Orr-type
wavepacket. The elongated structures tilted against the mean shear in the
corresponding forcing are a clear indicator of the Orr mechanism associated
with the optimal non-modal spatial growth [117, 137, 138]. Similar flow
structures have been reported in Garnaud _et al._ [44], Lesshafft _et al._
[69] and Pickering _et al._ [99]. For higher frequencies ($\mathrm{St}\gtrsim
0.6$), the optimal response take the form of compact Kelvin–Helmholtz (KH)
wavepackets localized upstream in the initial shear layer region of the jet
($x\lesssim 10$), which can be identified as the modal KH-type shear-layer
instability of the turbulent mean flow [48, 62, 134]. The corresponding
forcing distributions near the lip line ($r\simeq 0.5$) remain indicative of
the Orr mechanism, but this time with the KH instability [44, 59, 103, 118,
138]. Within the potential core, the optimal response exhibits duct modes at
$\mathrm{St}=0.6$ and trapped acoustic waves at $\mathrm{St}=1$. The presence
of the latter is a general feature of resonance mechanisms between propagating
waves associated with the isothermal and transonic jet, as previously
described in [116, 141]. Notably, comparable patterns can be observed in the
corresponding forcings, indicating the nearly self-adjoint nature of the
trapped instability mechanism. We confirmed the modal shapes are almost
identical to those obtained using the numerical scheme outlined in Schmidt _et
al._ [117] (not shown).

Figure 17:  RA-based wavemaker $\zeta_{RA}$ for: (a) $\mathrm{St}=0.2$; (b)
$\mathrm{St}=0.6$; (c) $\mathrm{St}=1$. The region with the strongest feedback
is marked as green ’+’.

Figure 17 shows how the RA-based wavemaker, $\zeta_{RA}$, unveils different
physical mechanisms that are active in the turbulent jet. Most notably, the
wavemaker region and its overall peak move upstream as the frequency increases
from $\mathrm{St}=0.2$ to 1. At $\mathrm{St}=0.2$, the optimal forcing and
response modes are spatially separated (see panels 16(a,b)), resulting in a
very weak wavemaker signature in the downstream self-similar region. This
result is anticipated, as the responses at low frequencies are purely
triggered by the Orr mechanism without being associated with a single global
mode or a local feedback mechanism that leads to the creation of waves in a
localized region. This is indicative of a non-modal convective instability. At
$\mathrm{St}=0.6$, the wavemaker peaks at the centerline at $x\simeq 5$. This
is the end of the potential core, where upstream traveling acoustic modes are
generated, also known as duct modes, as previously described by Towne _et al._
[141]. Parallel-flow models accurately predict the occurrence of duct modes,
but the wavemaker potentially reveals the location where duct modes originate,
which is not predicted by the theory. The wavemaker at $\mathrm{St}=1$ peaks
within the shear layer, which is associated with the KH instability, but it
also exhibits a comparable magnitude at the upstream region of the potential
core, which identifies the resonance mechanisms that trigger trapped acoustic
modes [116, 141]. That is, the wavemaker analysis confirms the phase-linking
between downstream KH waves and upstream traveling waves at $\mathrm{St}=1$.
Note that this is not the case for $\mathrm{St}=0.6$, where the wavemaker is
solely associated with the duct modes.

##  6 Summary and discussion

In this study, a novel higher-order mesh-free framework for hydrodynamic
stability analysis is developed. The framework has been demonstrated and
validated on three benchmark problems for open flows: a canonical laminar
cylinder wake, a non-parallel flat-plate Blasius boundary layer, and a
transonic turbulent jet. It also provides new insights into well-known
physics. The proposed framework incorporates PHS-type RBFs with polynomial
augmentations to discretize large hydrodynamic stability matrix problems on
scattered nodes with high accuracy, stability, and computational efficiency.
The resulting differentiation matrices are employed to construct the discrete
linearized Navier-Stokes operator for conducting LST and RA. We propose a set
of parameters to address the trade-off between accuracy and computational
efficiency for PHS+Poly RBF-FD discretizations that provide best practices. In
the context of scattered nodes, the practical implementations of various
boundary conditions arising in hydrodynamic stability analysis are discussed
and addressed.

The present mesh-free approach accurately predicts flow instabilities in wakes
behind a 2D cylinder, including the vortex-shedding frequency, corresponding
coherent structures, and structural sensitivity. The results are in good
agreement with previous studies [11, 46, 76, 100]. In the case of the Blasius
boundary layer, the present framework accurately identifies the TS wavepackets
and yields favorable quantitative and qualitative comparisons with those
previously reported in the literature [1, 22, 89, 128]. Furthermore, the RA-
based wavemaker confirms the high non-normality of the boundary layer and the
transition in dominant mechanisms from TS to Orr as the frequency increases.
When applied to the mean flow of a turbulent jet, the mesh-free RA yields
results that are almost identical to those reported by Schmidt _et al._ [117].
The application of RA-based WM analysis identifies distinct dominant physical
mechanisms operating at different frequencies. These mechanisms include the
Orr mechanism in the downstream region with a very weak wavemaker signature at
$\mathrm{St}=0.2$, prominent duct-like structures near the end of the
potential core at $\mathrm{St}=0.6$, and the phase coupling between the KH and
trapped acoustic mechanisms at $\mathrm{St}=1$. These results closely align
with the findings documented in prior literature, e.g., [62, 69, 116, 117,
141], and provide a new, global perspective on modal and non-modal growth in
jets, as well as on the origins of these instabilities.

Building upon the pioneering work of Hardy [50], RBF discretizations have been
widely employed in canonical problems, exhibiting a continual evolution
through advancements like RBF-FDs [139] and PHS+poly RBFs [39]. More recently,
this mesh-free approach has garnered attention in computational fluid dynamics
(CFD) for physical exploration or engineering applications, see, e.g., [29,
119, 120, 121, 147]. While RBF-FDs demonstrate the capability to construct
sparse differentiation operators, the potential applications in hydrodynamic
stability analysis have not been thoroughly explored. The current work is the
first demonstration of such applications. Future extensions encompass
applications involving 3D and complex geophysical or engineering fluid flows.

## CRediT authorship contribution statement

Tianyi Chu: Formal analysis, Methodology, Software, Validation, Visualization,
Writing – original draft. Oliver T. Schmidt: Conceptualization, Funding
acquisition, Methodology, Project administration, Resources, Supervision,
Writing – review $\&$ editing.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships
which may be considered as potential competing interests: Oliver T. Schmidt
reports financial support was provided by National Science Foundation under
Grant No. CBET-1953999. Tianyi Chu reports financial support was provided by
National Science Foundation under Grant No. CBET1953999.

## Data availability

Data will be made available on request.

## Acknowledgments

We gratefully acknowledge support by the National Science Foundation under
Grant No. CBET-1953999 (PM Ron Joslin).

##  Appendix A Governing equations for compressible flows

The compressible Navier-Stokes equations govern the motion of a general,
compressible Newtonian fluid,

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\nabla\cdot\rho\mathbf{u},$ (A.1)
$\displaystyle\frac{\partial\rho\mathbf{u}}{\partial t}$
$\displaystyle=-\frac{1}{2}\nabla\cdot(\mathbf{u}:\rho\mathbf{u}+\rho\mathbf{u}:\mathbf{u})-\nabla
p+\frac{1}{\mathrm{Re}}\nabla\cdot{\bf\it\tau},$
$\displaystyle\frac{\partial\rho e}{\partial t}$
$\displaystyle=-\nabla\cdot\rho
e\mathbf{u}+\frac{1}{(\gamma-1)\mathrm{Re}\mathrm{Pr}\mathrm{Ma}^{2}_{\infty}}\nabla\cdot
k\nabla T-\nabla\cdot
p\mathbf{u}+\frac{1}{\mathrm{Re}}\nabla\cdot{\bf\it\tau}\mathbf{u},$

where, $e$ is the total energy. For a Newtonian fluid, the viscous stress
tensor is
${\bf\it\tau}=\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)-\frac{2}{3}\mu(\nabla\cdot\mathbf{u})\mathrm{I}$.
All flow quantities are non-dimensionalized by their dimensional free-stream
values, denoted by $(\cdot)_{\infty}^{*}$, and the coordinates by the jet
diameter $D$. The dimensionless Reynolds number
$\mathrm{Re}={\rho_{\infty}^{*}u_{\infty}^{*}D}/{\mu_{\infty}^{*}}$, Prandtl
number $\mathrm{Pr}={\mathrm{c_{p}^{*}}\mu_{\infty}^{*}}/{k_{\infty}^{*}}$,
and Mach number $M={u_{\infty}^{*}}/{a_{\infty}^{*}}$ then fully describe the
flow. Here, $\mu_{\infty}^{*},k_{\infty}^{*},c_{p}^{*},\gamma,a_{\infty}^{*}$
are the free-stream values of the dynamic viscosity, heat conductivity, heat
capacity at constant pressure, heat capacity ratio, and speed of sound,
respectively. Closure of the equations is achieved under the assumption of an
ideal gas, and using Sutherland’s law to compute the dynamic viscosity from
the local temperature.

## References

  * Åkervik _et al._ [2008] E. Åkervik, U. Ehrenstein, F. Gallaire, D. S. Henningson, Global two-dimensional stability measures of the flat plate boundary-layer flow, Eur. J. Mech. B-Fluids. 27 (5)(2008) 501–513.
  * Alizard and Robinet [2007] F. Alizard, J.-C. Robinet, Spatially convective global modes in a boundary layer, Phys. Fluids. 19 (11)(2007) 114105.
  * Alizard _et al._ [2009] F. Alizard, S. Cherubini, J.-C. Robinet, Sensitivity and optimal forcing response in separated boundary layer flows, Phys. Fluids. 21 (6)(2009) 064108.
  * Alizard _et al._ [2022] F. Alizard, T. Gibis, B. Selent, U. Rist, C. Wenzel, Stochastic receptivity of laminar compressible boundary layers: An input-output analysis, Phys. Rev. Fluids. 7 (7)(2022) 073902.
  * Andersson _et al._ [1999] P. Andersson, M. Berggren, D. S. Henningson, Optimal disturbances and bypass transition in boundary layers, Phys. Fluids. 11 (1)(1999) 134–150.
  * Bae _et al._ [2020] H. J. Bae, S. T. M. Dawson, B. J. McKeon, Resolvent-based study of compressibility effects on supersonic turbulent boundary layers, J. Fluid Mech. 883 (2020) A29.
  * Bagheri _et al._ [2009] S. Bagheri, E. Åkervik, L. Brandt, D. S. Henningson, Matrix-free methods for the stability and control of boundary layers, AIAA J. 47 (5)(2009) 1057–1068 .
  * Bagheri _et al._ [2009] S. Bagheri, L. Brandt, D. S. Henningson, Input–output analysis, model reduction and control of the flat-plate boundary layer, J. Fluid Mech. 620 (2009) 263–298.
  * Bagheri _et al._ [2009] S. Bagheri, D. S. Henningson, J. Hoepffner, P. J. Schmid, Input-output analysis and control design applied to a linear model of spatially developing flows, Appl. Mech. Rev. 62 (2)(2009).
  * Bagheri _et al._ [2009] S. Bagheri, P. Schlatter, P. J. Schmid, D. S. Henningson, J. Hoepffner, Global stability of a jet in crossflow, J. Fluid Mech. 624 (2009) 33–44. J. Fluid Mech. 624 (2009) 33-44.
  * Barkley [2006] D. Barkley, Linear analysis of the cylinder wake mean flow, Europhys. Lett. 75 (5)(2006) 750.
  * Barkley and Henderson [1996] D. Barkley, R. D. Henderson, Three-dimensional Floquet stability analysis of the wake of a circular cylinder, J. Fluid Mech. 322 (1996) 215–241.
  * Bayona [2019a] V. Bayona, An insight into RBF-FD approximations augmented with polynomials, Comput. Math. Appl. 77 (9)(2019) 2337–2353.
  * Bayona [2019b] V. Bayona, Comparison of moving least squares and RBF+ poly for interpolation and derivative approximation, J. Sci. Comput. 81 (1)(2019) 486–512.
  * Bayona _et al._ [2017] V. Bayona, N. Flyer, B. Fornberg, G.A. Barnett, On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs, J. Comput. Phys. 332 (2017) 257–273.
  * Bayona _et al._ [2019] V. Bayona, N. Flyer, B. Fornberg, On the role of polynomials in RBF-FD approximations: III. Behavior near domain boundaries, J. Comput. Phys. 380 (2019) 378–399.
  * Bénard [1908] H. Bénard, The formation of gyration centres at the back of a moving obstacle, C. R. Acad. Sci. 147 (1908) 839–842.
  * Beneddine _et al._ [2016] S. Beneddine, D. Sipp, A. Arnault, J. Dandois, L. Lesshafft, Conditions for validity of mean flow stability analysis, J. Fluid Mech. 798 (2016) 485–504.
  * Berlin and Henningson [1999] S. Berlin, D. S. Henningson, A nonlinear mechanism for receptivity of free-stream disturbances, Phys. Fluids 11 (12)(1999) 3749–3760.
  * Betchov and Szewczyk [1963] R. Betchov, A. Szewczyk, Stability of a shear layer between parallel streams, Phys. Fluids 6 (10)(1963) 1391–1396.
  * Bonne _et al._ [2019] N. Bonne, V. Brion, E. Garnier, R. Bur, P. Molton, D. Sipp, L. Jacquin, Analysis of the two-dimensional dynamics of a Mach 1.6 shock wave/transitional boundary layer interaction using a RANS based resolvent approach, J. Fluid Mech. 862 (2019) 1166–1202.
  * Brandt _et al._ [2011] L. Brandt, D. Sipp, J. O. Pralits, O. Marquet, Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer, J. Fluid Mech. 687 (2011) 503–528.
  * Brès _et al._ [2018] G. A. Brès, and P. Jordan, V. Jaunet, M. Le Rallic, A. V. G. Cavalieri, A. Towne, S. K. Lele, T. Colonius, O. T. Schmidt, Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets, J. Fluid Mech. 851 (2018) 83-124.
  * Bugeat _et al._ [2019] B. Bugeat, J.-C. Chassaing, J.-C. Robinet, P. Sagaut, 3D global optimal forcing and response of the supersonic boundary layer, J. Comput. Phys. 398 (2019) 108888.
  * Butler and Farrell [1992] K. M. Butler, B. F. Farrell, Three-dimensional optimal perturbations in viscous shear flow, Phys. Fluids A. 4 (8)(1992) 1637–1650.
  * Cavalieri _et al._ [2013] A. V. G. Cavalieri, D. Rodríguez, P. Jordan, T. Colonius, Y. Gervais , Wavepackets in the velocity field of turbulent jets, J. Fluid Mech. 730 (2013) 559–592.
  * Chomaz [2005] J.-M. Chomaz, Global instabilities in spatially developing flows: non-normality and nonlinearity, Annu. Rev. Fluid Mech. 37 (2005) 357–392.
  * Chu [1965] B.-T. Chu, On the energy transfer to small disturbances in fluid flow (Part 1), Acta Mech. 1 (3)(1965) 215–234.
  * Chu and Schmidt [2023] T. Chu, O. Schmidt, RBF-FD discretization of the Navier-Stokes equations on scattered but staggered nodes, J. Comput. Phys. 474 (2023) 111756.
  * Crighton and Gaster [1976] D. G. Crighton, M. Gaster, Stability of slowly diverging jet flow, J. Fluid Mech. 77 (2)(1976) 397–413.
  * Dawson and McKeon [2020] S. T. M. Dawson, B. J. McKeon, Prediction of resolvent mode shapes in supersonic turbulent boundary layers, Int. J. Heat Fluid Flow. 85 (2020) 108677.
  * De _et al._ [2021] P. De, F. Miguel, D. Sipp, P. J. Schmid, Efficient evaluation of the direct and adjoint linearized dynamics from compressible flow solvers, J. Comput. Phys. 231 (23)(2012) 7739–7755.
  * Dwivedi _et al._ [2022] A. Dwivedi, G. S. Sidharth, M. R. Jovanović, Oblique transition in hypersonic double-wedge flow, J. Fluid Mech. 948 (2022) A37.
  * Edstrand _et al._ [2018] A. D. Edstrand, P. J. Schmid, K. Taira, L. N. Cattafesta, A parallel stability analysis of a trailing vortex wake, J. Fluid Mech. 837 (2018) 858–895.
  * Ehrenstein and Gallaire [2005] U. Ehrenstein, F. Gallaire, On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer, J. Fluid Mech. 536 (2005) 209–218.
  * Eriksson and Rizzi [1985] L. E. Eriksson, A. Rizzi, Computer-aided analysis of the convergence to steady state of discrete approximations to the Euler equations, J. Comput. Phys. 57 (1)(1985) 90–128.
  * Farrell and Ioannou [2005] B. F. Farrell, P. J. Ioannou, Stochastic forcing of the linearized Navier–Stokes equations, Phys. Fluids A 5 (11)(1993) 2600–2609.
  * Fey _et al._ [1998] U. Fey, M. König, H. Eckelmann, A new Strouhal–Reynolds-number relationship for the circular cylinder in the range $47<\mathrm{Re}<2\times 10^{5}$, Phys. Fluids 10 (7)(1998) 1547–1549.
  * Flyer _et al._ [2016] N. Flyer, G.A. Barnett, L.J. Wicker, Enhancing finite differences with radial basis functions: experiments on the Navier–Stokes equations, J. Comput. Phys. 316 (2016) 39–62.
  * Flyer _et al._ [2016] N. Flyer, B. Fornberg, V. Bayona, G.A. Barnett, On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy, J. Comput. Phys. 321 (2016) 21–38.
  * Fornberg _et al._ [2002] B. Fornberg, T. A. Driscoll, G. Wright, R. Charles, Observations on the behavior of radial basis function approximations near boundaries, Comput. Math. Appl. 43 (3-5)(2002) 473–490.
  * Fornberg and Flyer [2015] B. Fornberg, N. Flyer, Solving PDEs with radial basis functions, Acta Numer. 24 (2015) 215–258.
  * Fornberg and Letho [2011] B. Fornberg, E. Lehto, Stabilization of RBF-generated finite difference methods for convective PDEs, J. Comput. Phys. 230 (6)(2011) 2270–2285.
  * Garnaud _et al._ [2013] X. Garnaud, L. Lesshafft, P. J. Schmid, P. Huerre, The preferred mode of incompressible jets: linear frequency response analysis, J. Fluid Mech. 716 (2013) 189–202.
  * Gaster _et al._ [1985] M. Gaster, E. Kit, I. Wygnanski, Large-scale structures in a forced turbulent mixing layer, J. Fluid Mech. 150 (1985) 23–39.
  * Giannetti and Luchini [2007] F. Giannetti, P. Luchini, Structural sensitivity of the first instability of the cylinder wake, J. Fluid Mech. 581 (2007) 167–197.
  * Gunderman _et al._ [2020] D. Gunderman, N. Flyer, B. Fornberg, Transport schemes in spherical geometries using spline-based RBF-FD with polynomials, J. Comput. Phys. 408 (2020) 109256.
  * Gudmundsson and Colonius [2011] K. Gudmundsson, T. Colonius, Instability wave models for the near-field fluctuations of turbulent jets, J. Fluid Mech. 689 (2011) 97–128.
  * Hammond and Redekopp [1997] D. A. Hammond, L. G. Redekopp, Global dynamics of symmetric and asymmetric wakes, J. Fluid Mech. 331 (1997) 231–260.
  * Hardy [1971] R. L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 76 (8)(1971) 1905–1915.
  * Hecht [2012] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (3-4)(2012) 251–266.
  * Ho and Huerre [1984] C. -M. Ho, P. Huerre, Perturbed free shear layers, Annu. Rev. Fluid Mech. 16 (1)(1984) 365–422.
  * Huerre and Monkewitz [1990] P. Huerre, P. A. Monkewitz, Local and global instabilities in spatially developing flows, Annu. Rev. Fluid Mech. 22 (1)(1990) 473–537.
  * Hwang and Cossu [2010] Y. Hwang, C. Cossu, Amplification of coherent streaks in the turbulent Couette flow: an input–output analysis at low Reynolds number, J. Fluid Mech. 643 (2010) 333–348.
  * Van Ingen [1956] J. L. Van Ingen, A suggested semi-empirical method for the calculation of the boundary layer transition region, Rep. UTHL-74. University of Technology, Dept. Aero. Engng.
  * Iske [2003a] A. Iske, On the approximation order and numerical stability of local Lagrange interpolation by polyharmonic splines, in: W. Haussmann, K. Jetter, M. Reimer, J. Stöckler (Eds.), Modern developments in multivariate approximation, in: International Series of Numerical Mathematics, Birkh’́auser, Basel. (2003) 153–165.
  * Jackson [1987] C. P. Jackson, A finite-element study of the onset of vortex shedding in flow past variously shaped bodies, J. Fluid Mech. 182 (1987) 23–45.
  * Jacobi and McKeon [2011] I. Jacobi, B. J. McKeon, Dynamic roughness perturbation of a turbulent boundary layer, J. Fluid Mech. 688 (2011) 258–296.
  * Jeun _et al._ [2016] J. Jeun, J. W. Nichols, M. R. Jovanović, Input-output analysis of high-speed axisymmetric isothermal jet noise, Phys. Fluids. 28 (4)(2016) 047101.
  * Jiang and Cheng [2017] H. Jiang, L. Cheng, Strouhal–Reynolds number relationship for flow past a circular cylinder, J. Fluid Mech. 832 (2017) 170–188.
  * Jin _et al._ [2021] B. Jin, S. Symon, S. J. Illingworth, Energy transfer mechanisms and resolvent analysis in the cylinder wake, Phys. Rev. Fluids. 6 (2)(2021) 024702.
  * Jordan and Colonius [2013] P. Jordan, T. Colonius, Wave packets and turbulent jet noise, Annu. Rev. Fluid Mech. 45 (2013) 173–195.
  * Jordinson [1990a] R. Jordinson, The flat plate boundary layer. Part 1. Numerical integration of the Orr—Sommerfeld equation, J. Fluid Mech. 43 (4)(1970) 801–811.
  * Jovanović and Bamieh [2005] M. R. Jovanović, B. Bamieh, Componentwise energy amplification in channel flows, J. Fluid Mech. 534 (2005) 145–183.
  * Kansa [1990b] E. J. Kansa, Multiquadrics—A scattered data approximation scheme with applications to computational fluid-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations, Comput. Math. Appl. 19 (8-9)(1990) 147–161.
  * Larsson _et al._ [2013] E. Larsson, E. Lehto, A. Heryudono, B. Fornberg, Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions, SIAM J. Sci. Comput. 35 (4)(2013) 2096–2119.
  * Le Borne and Leinen [2023] S. Le Borne, W. Leinen, Guidelines for RBF-FD Discretization: Numerical Experiments on the Interplay of a Multitude of Parameter Choices, J. Sci. Comput. 95 (1)(2023) 8.
  * Lele [1992] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1)(1992) 16–42.
  * Lesshafft _et al._ [2019] L. Lesshafft, O. Semeraro, V. Jaunet, A. V. G. Cavalieri, P. Jordan, Resolvent-based modeling of coherent wave packets in a turbulent jet, Phys. Rev. Fluid 4 (6)(2019) 063901.
  * Luchini [2000] P. Luchini, Reynolds-number-independent instability of the boundary layer over a flat surface: optimal perturbations, J. Fluid Mech. 404 (2000) 289–309.
  * Luchini and Bottaro [2014] P. Luchini, A. Bottaro, Adjoint equations in stability analysis, Annu. Rev. Fluid Mech. 46 (2014) 493–517.
  * Mack [1976] M. L. Mack, A numerical study of the temporal eigenvalue spectrum of the Blasius boundary layer, J. Fluid Mech. 73 (3)(1976) 497–520.
  * Mack [1984] M. L. Mack, Boundary-layer linear stability theory, Technical Report No. AGARD-R-709, 1984.
  * Mantič-Lugo _et al._ [2014] V. Mantič-Lugo, C. Arratia, F. Gallaire, Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake, Phys. Rev. Lett. 113 (8)(2014) 084501.
  * Mantič-Lugo _et al._ [2014] V. Mantič-Lugo, C. Arratia, F. Gallaire, A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake, Phys. Fluids. 27 (7)(2015) 074103.
  * Marquet _et al._ [2008] O. Marquet, D. Sipp, L. Jacquin, Sensitivity analysis and passive control of cylinder flow, J. Fluid Mech. 615 (2008) 221–252.
  * Marquet _et al._ [2009] O. Marquet, M. Lombardi, J.-M.. Chomaz, D. Sipp, L. Jacquin, Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities, J. Fluid Mech. 622 (2009) 1–21.
  * Martini _et al._ [2021] E. Martini, D. Rodríguez, A. Towne, A. V. G. Cavalieri, Efficient computation of global resolvent modes, J. Fluid Mech. 919 (2021) A3.
  * Mattingly and Criminale [2010] G. E. Mattingly, W. O. Criminale, The stability of an incompressible two-dimensional wake, J. Fluid Mech. 51 (2)(1972) 233–272.
  * Mattsson and Nordström [2010] K. Mattsson, J. Nordström, Summation by parts operators for finite difference approximations of second derivatives, J. Comput. Phys. 199 (2)(2004) 503–540.
  * McKeon and Sharma [2010] B. J. McKeon, A. S. Sharma, A critical-layer framework for turbulent pipe flow, J. Fluid Mech. 658 (2010) 336–382.
  * Meliga _et al._ [2016] P. Meliga, E. Boujo, F. Gallaire, A self-consistent formulation for the sensitivity analysis of finite-amplitude vortex shedding in the cylinder wake, J. Fluid Mech. 800 (2016) 327–357.
  * Mirzaei [2021] D. Mirzaei, The direct radial basis function partition of unity (D-RBF-PU) method for solving PDEs, SIAM J. Sci. Comput. 43 (1)(2021) A54–A83.
  * Moarref and Jovanović [2012] R. Moarref, M. R. Jovanović, Model-based design of transverse wall oscillations for turbulent drag reduction, J. Fluid Mech. 707 (2012) 205–240.
  * Moarref _et al._ [2014] R. Moarref, M. R. Jovanović, J. A. Tropp, A. S. Sharma, B. J. McKeon, A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization, Phys. Fluids 26 (5)(2014) 051701.
  * Mohseni and Colonius [2000] K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, J. Comput. Phys. 157 (2)(2000) 787–795.
  * Monkewitz [1988] P. A. Monkewitz, Subharmonic resonance, pairing and shredding in the mixing layer, J. Fluid Mech. 188 (1988) 223–252.
  * Monkewitz _et al._ [1993] P. A. Monkewitz, P. Huerre, J.-M. Chomaz, Global linear stability analysis of weakly non-parallel shear flows, J. Fluid Mech. 251 (1993) 1–20.
  * Monokrousos _et al._ [2010] A. Monokrousos, E. Åkervik, L. Brandt, D. S. Henningson, Global three-dimensional optimal disturbances in the Blasius boundary-layer flow using time-steppers, J. Fluid Mech. 650 (2010) 181–214.
  * Moarref and Jovanović [2010] R. Moarref, M. R. Jovanović, Controlling the onset of turbulence by streamwise travelling waves. Part 1. Receptivity analysis, J. Fluid Mech. 663 (2010) 70–99.
  * Noack and Eckelmann [1994] B. R. Noack, H. Eckelmann, A global stability analysis of the steady and periodic cylinder wake, J. Fluid Mech. 270 (1994) 297–330.
  * Oberleithner _et al._ [2014] K. Oberleithner, L. Rukes, J. Soria, Mean flow stability analysis of oscillating jet experiments, J. Fluid Mech. 757 (2014) 1–32.
  * Ohmichi and Yamada [2021] Y. Ohmichi, K. Yamada, Matrix-free triglobal adjoint stability analysis of compressible Navier–stokes equations, J. Comput. Phys. 437 (2021) 110332.
  * Oruç [2022] Ö. Oruç, A strong-form local meshless approach based on radial basis function-finite difference (RBF-FD) method for solving multi-dimensional coupled damped Schrödinger system appearing in Bose–Einstein condensates, Commun. Nonlinear Sci. Numer. Simul. 104 (2022) 106042.
  * Peplinski _et al._ [2015] A. Peplinski, P. Schlatter, D. S. Henningson, Global stability and optimal perturbation for a jet in cross-flow, Eur. J. Mech. B/Fluids. 949 (2015) 438–447.
  * Persson [2005] P.O. Persson, Ph.D. thesis, Massachusetts Institute of Technology, 2005.
  * Pfister _et al._ [2022] J.-L. Pfister, N. Fabbiane, O. Marquet, Global stability and resolvent analyses of laminar boundary-layer flow interacting with viscoelastic patches, J. Fluid Mech. 937 (2022) A1.
  * Pickering _et al._ [2020] E. Pickering, G. Rigas, P. A. S. Nogueira, A. V. G. Cavalieri, O. T. Schmidt, T. Colonius, Lift-up, Kelvin–Helmholtz and Orr mechanisms in turbulent jets, J. Fluid Mech. 896 (2020) A2.
  * Pickering _et al._ [2021] E. Pickering, G. Rigas, O. T. Schmidt, D. Sipp, T. Colonius, Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets, J. Fluid Mech. 917 (2021) A29.
  * Pier [1992] B. Pier, On the frequency selection of finite-amplitude vortex shedding in the cylinder wake, J. Fluid Mech. 458 (2002) 407–417.
  * Ponta and Hassan [2004] F. L. Ponta, H. Aref, Strouhal-Reynolds number relationship for vortex streets, Phys. Rev. Lett. 93 (8)(2004) 084501.
  * Provansal _et al._ [1996] M. Provansal, C. Mathis, L. Boyer, Bénard-von Kármán instability: transient and forced regimes, J. Fluid Mech. 182 (1987) 1–22.
  * Qadri and Schmid [2017] U. A. Qadri, P. J. Schmid, Effect of nonlinearities on the frequency response of a round jet, Phys. Rev. Fluids 2 (4) (2017) 043902.
  * Reddy and Henningson [1993] S. C. Reddy, D. S. Henningson, Energy growth in viscous channel flows, J. Fluid Mech. 252 (1993) 209–238.
  * Reddy _et al._ [1993] S. C. Reddy, P. J. Schmid, D. S. Henningson, Pseudospectra of the Orr–Sommerfeld operator, SIAM J. Appl. Math. 53 (1)(1993) 15–47.
  * Reed _et al._ [1996] H. L. Reed, W. S. Saric, D. Arnal, Linear stability theory applied to boundary layers, Annu. Rev. Fluid Mech. 28 (1)(1996) 389–428.
  * Regan and Mahesh [2017] M . A. Regan, K. Mahesh, Global linear stability analysis of jets in cross-flow, J. Fluid Mech. 828 (2017) 812–836.
  * Ribeiro _et al._ [2020] J. H. M. Ribeiro, C.-A. Yeh, K. Taira, Randomized resolvent analysis, Phys. Rev. Fluids 5 (3) (2020) 033902.
  * Ribeiro _et al._ [2022] J. H. M. Ribeiro, C.-A. Yeh, K. Zhang, K. Taira, Wing sweep effects on laminar separated flows, J. Fluid Mech. 950 (2022) A23.
  * Rigas _et al._ [2021] G. Rigas, D. Sipp, T. Colonius, Nonlinear input/output analysis: application to boundary layer transition, J. Fluid Mech. 911 (2021) A15.
  * Rosenberg and McKeon [2019] K. Rosenberg, B. J. McKeon, Computing exact coherent states in channels starting from the laminar profile: A resolvent-based approach, Phys. Rev. E 100 (2)(2019) 021101.
  * Santos _et al._ [2018] L.G.C. Santos, N. Manzanares-Filho, G.J. Menon, E. Abreu, Comparing RBF-FD approximations based on stabilized Gaussians and on polyharmonic splines with polynomials, Internat. J. Numer. Methods Engrg. 115 (4)(2018) 462–500.
  * Schmid [2007] P. J. Schmid, Nonmodal stability theory, Annu. Rev. Fluid Mech. 39 (2007) 129–162.
  * Schmid and Brandt [2002] P. J. Schmid, L. Brandt, Analysis of fluid systems: Stability, receptivity, sensitivity, Appl. Mech. Rev. 66 (2)(2014)
  * Schmid and Henningson [2001] P. J. Schmid, D. S. Henningson, Stability and transition in shear flows, Springer.
  * Schmidt _et al._ [2017] O. T. Schmidt, A. Towne, T. Colonius, A. V. Cavalieri, P. Jordan, G. A. Brès, Wavepackets and trapped acoustic modes in a turbulent jet: coherent structure eduction and global stability, J. Fluid Mech. 825 (2017) 1153–1181.
  * Schmidt _et al._ [2018] O. T. Schmidt, and A. Towne, G. Rigas, T. Colonius, G. A. Brès Spectral analysis of jet turbulence, J. Fluid Mech. 855 (2018) 953–982.
  * Semeraro _et al._ [2016] O. Semeraro, L. Lesshafft, V. Jaunet, P. Jordan Modeling of coherent structures in a turbulent jet as global linear instability wavepackets: theory and experiment, Int. J. Heat Fluid Flow. 62 (2016) 24–32.
  * Shahane _et al._ [2021] S. Shahane, A. Radhakrishnan, S.P. Vanka, A high-order accurate meshless method for solution of incompressible fluid flow problems, J. Comput. Phys. 445 (2021) 110623.
  * Shahane and Vanka [2022] S. Shahane, S.P. Vanka, Consistency and Convergence of a High Order Accurate Meshless Method for Solution of Incompressible Fluid Flows, arXiv preprint arXiv:2202.02828.
  * Shahane and Vanka [2023] S. Shahane, S.P. Vanka, A semi-implicit meshless method for incompressible flows in complex geometries, J. Comput. Phys. 472 (2023) 111715.
  * Shankar [2017] V. Shankar, The overlapped radial basis function-finite difference (RBF-FD) method: A generalization of RBF-FD, J. Comput. Phys. 342 (2017) 211–228.
  * Shankar and Fogelson [2018] V. Shankar, A. L. Fogelson, Hyperviscosity-based stabilization for radial basis function-finite difference (RBF-FD) discretizations of advection–diffusion equations, J. Comput. Phys. 372 (2018) 616–639.
  * Sharma and McKeon [2013] A. S. Sharma, B. J. McKeon, On coherent structure in wall turbulence, J. Fluid Mech. 728 (2013) 196–238.
  * Sidharth _et al._ [2016] G. S. Sidharth, A. Dwivedi, G. V. Candler, J. W. Nichols, Onset of three-dimensionality in supersonic flow over a slender double wedge, Phys. Rev. Fluids 3 (9)(2018) 093901.
  * Sipp [2012] D. Sipp, Open-loop control of cavity oscillations with harmonic forcings, J. Fluid Mech. 708 (2012) 439–468.
  * Sipp and Marquet [2007] D. Sipp, A. Lebedev, Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows, J. Fluid Mech. 593 (2007) 333-358.
  * Sipp and Marquet [2013] D. Sipp, O. Marquet, Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer, Theor. Comput. Fluid Dyn. 27 (5)(2013) 617–635.
  * Sipp _et al._ [2010] D. Sipp, O. Marquet, P. Meliga, A. Barbagallo, Dynamics and control of global instabilities in open-flows: a linearized approach, Appl. Mech. Rev. 63 (3) (2010).
  * Skene _et al._ [2022] C. S. Skene, C.-A. Yeh, P. J. Schmid, K. Taira, Sparsifying the resolvent forcing mode via gradient-based optimisation, J. Fluid Mech. 944 (2022) A52.
  * Smith and Gamberoni [1956] A. M. O. Smith, N. Gamberoni, Transition, pressure gradient and stability theory, Rep. ES 26388. Douglas Aircraft Co., El Segundo, CA.
  * Smith [1979] F. T. Smith, On the non-parallel flow stability of the Blasius boundary layer, Proc. R. Soc. London Ser. A 366 (1724)(1979) 91–109.
  * Sreenivasan _et al._ [1987] K. R. Sreenivasan, P. J. Strykowski, D. J. Olinger, Hopf bifurcation, Landau equation, and vortex shedding behind circular cylinders, Forum on unsteady flow separation. 1 (1987) 1–13.
  * Suzuki and Colonius [2006] T. Suzuki, T. Colonius, Instability waves in a subsonic round jet detected using a near-field phased microphone array, J. Fluid Mech. 565 (2006) 197–226.
  * Symon _et al._ [2018] S. Symon, K. Rosenberg, S. T. M. Dawson, B. J. McKeon, Non-normality and classification of amplification mechanisms in stability and resolvent analysis, Phys. Rev. Fluids. 3 (5)(2018) 053902.
  * Theofilis [2011] V. Theofilis, Global linear instability, Annu. Rev. Fluid Mech. 43 (2011) 319–352.
  * Tissot _et al._ [2017] G. Tissot, F. C. Lajús Jr., A. V. G. Cavalieri, P. Jordan, Wave packets and Orr mechanism in turbulent jets, Phys. Rev. Fluids 2 (9)(2017) 093901.
  * Tissot _et al._ [2017] G. Tissot, M. Zhang, F. C. Lajús Jr., A. V. G. Cavalieri, P. Jordan, Sensitivity of wavepackets in jets to nonlinear effects: the role of the critical layer, J. Fluid Mech. 811 (2017) 95–137.
  * Tolstykh [2000] A.I. Tolstykh, On using RBF-based differencing formulas for unstructured and mixed structured-unstructured grid calculations, Proc. 16th IMACS World Congress. 228 (2000) 4606–4624.
  * Tominec _et al._ [2022] I. Tominec, E. Larsson, A. Heryudono, A least squares radial basis function finite difference method with improved stability properties, SIAM J. Sci. Comput. 43 (2)(2021) A1441–A1471.
  * Towne _et al._ [2017] A. Towne, A. Cavalieri, P. Jordan , T. Colonius, O. T. Schmidt, V. Jaunet, G. Brès, Acoustic resonance in the potential core of subsonic jets, J. Fluid Mech. 825 (2017) 1113–1152.
  * Towne _et al._ [2022] A. Towne, G. Rigas, O. Kamal, E. Pickering, T. Colonius, Efficient global resolvent analysis via the one-way Navier–Stokes equations, J. Fluid Mech. 948 (2022) A9.
  * Towne _et al._ [2018] A. Towne, O. T. Schmidt, T. Colonius, Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis, J. Fluid Mech. 847 (2018) 821–867.
  * Trefethen _et al._ [2018] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, T. A. Driscoll, Hydrodynamic stability without eigenvalues, Science. 261 (5121)(1993) 578–584.
  * Tuckerman and Barkley [2000] L. Tuckerman, D. Barkley, Numerical methods for bifurcation problems and large-scale dynamical systems, Numerical methods for bifurcation problems and large-scale dynamical systems. 119 (2000) 453-466.
  * Tuckerman and Marcus [1985] L. Tuckerman, P. Marcus, Formation of Taylor vortices in spherical Couette flow, Ninth Int’l Conf. on Numerical Methods in Fluid Dynamics. (1985) 552–556.
  * Unnikrishnan _et al._ [2022] A. Unnikrishnan, S. Shahane, V. Narayan, S. P. Vanka, Shear-driven flow in an elliptical enclosure generated by an inner rotating circular cylinder, Phys. Fluids. 34 (1)(2022) 013607.
  * Vón Kármán [1911] T. Vón Kármán, Über den Mechanismus des Widerstandes, den ein bewegter Körper in einer Flüssigkeit erfährt, Göttingen Nachr. Math. Phys. Kl. 1911 (1911) 509–517.
  * Weisbrot and Wygnanski [1988] I. Weisbrot, I. Wygnanski, On coherent structures in a highly excited mixing layer, J. Fluid Mech. 195 (1988) 137–159.
  * Williamson [1996] C. H. K. Williamson, Vortex dynamics in the cylinder wake, Annu. Rev. Fluid Mech. 28 (1)(1996) 477–539.
  * Woodley and Peake [1997] B. M. Woodley, N. Peake, Global linear stability analysis of thin aerofoil wakes, J. Fluid Mech. 339 (1997) 239–260.
  * Wright and Fornberg [2006] G.B. Wright, B. Fornberg, Scattered node compact finite difference-type formulas generated from radial basis functions, J. Comput. Phys. 212 (1)(2006) 99–123.
  * Yeh and Taira [2019] C.-A. Yeh, K. Taira, Resolvent-analysis-based design of airfoil separation control, J. Fluid Mech. 867 (2019) 572–610.
  * Zebib [1987] A. Zebib, Stability of viscous flow past a circular cylinder, J. Engng Math. 21 (2)(1987) 155–165.