Linear Resolvent Analysis (RA)

• Linear Resolvent Analysis is an input–output framework that maps harmonic forcing to state responses via a resolvent operator in linearized systems.
• It uses weighted singular-value decomposition to extract optimal forcing and response modes, enabling reduced-order models for flow control and coherent structure analysis.
• Recent advancements include randomized, matrix-free, and data-driven approaches that enhance computational efficiency and extend RA to non-stationary and turbulent flows.

Linear resolvent analysis (RA) is an input–output framework for linearized dynamical systems in which harmonic forcing is mapped to harmonic state or output response through a resolvent operator. For a linear time-invariant system

$\frac{dx}{dt}=A\,x+f,$

harmonic forcing $f(t)=\hat f e^{i\omega t}+c.c.$ yields $(i\omega I-A)\hat x=\hat f$, so the resolvent operator is

$R(\omega)=(i\omega I-A)^{-1},$

with the more general input–output transfer $H(\omega)=C(i\omega I-A)^{-1}$ when measured outputs $y=Cx$ are considered. In fluid mechanics, RA is typically constructed from the Navier–Stokes equations linearized about a base flow or turbulent mean flow, with the nonlinear terms interpreted as an intrinsic forcing. Its central purpose is to identify the most amplified forcings, the most receptive states, and the associated gains at each frequency or frequency–wavenumber pair (Herrmann et al., 2020, Rolandi et al., 2024).

1. Core operator formulation

In the standard formulation, the flow field is decomposed into a mean or base state and fluctuations. After linearization, the fluctuation equations take the form

$\partial_t q' = L[\bar q]\,q' + u',$

where $L[\bar q]$ is the linearized operator and $u'$ collects external actuation and, in the mean-flow formulation, nonlinear interaction terms. Fourier transformation in time, and often in homogeneous spatial directions, converts the problem into a frequency-domain forced linear system. For incompressible and compressible flows alike, the resulting input–output relation is written as

$\hat q = H(\omega)\,\hat u,$

or, when spanwise periodicity is present,

$f(t)=\hat f e^{i\omega t}+c.c.$0

In this form, RA characterizes how disturbances at a prescribed frequency $f(t)=\hat f e^{i\omega t}+c.c.$1 and, where relevant, wavenumber $f(t)=\hat f e^{i\omega t}+c.c.$2, are selectively amplified by the linearized dynamics (Yeh et al., 2018, Lopez-Doriga et al., 2022).

The operator viewpoint is especially important in turbulent flows. Rather than treating turbulence solely as a broadband stochastic phenomenon, RA isolates coherent amplification mechanisms embedded in the linearized equations. In wall-bounded turbulence, pipe flow, airfoil wakes, duct flows, reacting flows, and stenotic transitional flows, the nonlinear term is treated as a forcing on a linear amplifier. This interpretation is not merely formal: it underlies applications to coherent-structure modeling, sensor and actuator placement, and flow control, and it is one of the reasons resolvent analysis has remained closely connected to both stability theory and control-theoretic transfer-function analysis (Rolandi et al., 2024, Villié et al., 14 Jun 2026).

2. Weighted singular-value decomposition and modal meaning

RA becomes operational through a singular-value decomposition of the resolvent operator in an energy norm. In the generic weighted formulation, one defines an inner product

$f(t)=\hat f e^{i\omega t}+c.c.$3

and forms the weighted resolvent

$f(t)=\hat f e^{i\omega t}+c.c.$4

Its singular-value decomposition is

$f(t)=\hat f e^{i\omega t}+c.c.$5

The right singular vectors $f(t)=\hat f e^{i\omega t}+c.c.$6 are the optimal forcing modes, the left singular vectors $f(t)=\hat f e^{i\omega t}+c.c.$7 are the response modes, and the singular values $f(t)=\hat f e^{i\omega t}+c.c.$8 are gains. This yields the low-rank approximation

$f(t)=\hat f e^{i\omega t}+c.c.$9

When $(i\omega I-A)\hat x=\hat f$0, the rank-one approximation is especially informative because a single forcing–response pair dominates the amplification (Herrmann et al., 2020, Rolandi et al., 2024).

The choice of norm is physically consequential. In compressible flows, the Chu energy norm is a standard choice; in other settings, kinetic-energy or application-specific norms are used. In the bluff-body stabilized flame with conjugate heat transfer, for example, the response norm was defined through the squared reaction-rate field, while the forcing norm was defined from the body force in the fluid momentum equations. In that setting, the leading gain $(i\omega I-A)\hat x=\hat f$1 was obtained as the largest eigenvalue of a Hermitian weighted problem, and the corresponding optimal forcing and response modes identified the flow structures that maximize heat-release response under harmonic forcing (Chen et al., 2024).

The singular structure is also the basis for physical interpretation. In turbulent pipe flow, good agreement between leading SPOD and resolvent modes is observed in a large region of parameter space, particularly where the gain separation $(i\omega I-A)\hat x=\hat f$2 is large and the lift-up mechanism is active. In that region, the optimal forcing takes the form of quasi-streamwise vortices and the response takes the form of streaks, making RA a reduced-order model for near-wall coherent structures. In airfoil separation control, the dominant response mode distinguishes actuation that remains confined to the suction-side shear layer from actuation that extends into the wake, a distinction tied directly to whether the induced mixing promotes reattachment or undesirable wake excitation (Abreu et al., 2020, Yeh et al., 2018).

3. Equation-based, inverse-free, randomized, and matrix-free computation

The classical route to RA requires direct access to the linearized operator and, in large problems, repeated action of the resolvent and its adjoint on vectors. “An Invitation to Resolvent Analysis” emphasizes two implementation paths: matrix-forming assembly of the Jacobian and matrix-free action of the Jacobian and its adjoint. In either case, the resolvent is not formed by direct inversion in realistic problems; instead, one solves linear systems of the form

$(i\omega I-A)\hat x=\hat f$3

and then computes leading singular triplets with Krylov or randomized algorithms (Rolandi et al., 2024).

Several variants reduce this cost. Marques Ribeiro, Yeh, and Taira introduced randomized resolvent analysis, in which a random test matrix $(i\omega I-A)\hat x=\hat f$4 is used to sketch the action of the resolvent on a low-dimensional subspace. On a turbulent post-stall NACA 0012 airfoil with operator size $(i\omega I-A)\hat x=\hat f$5, MATLAB’s svds for $(i\omega I-A)\hat x=\hat f$6 took $(i\omega I-A)\hat x=\hat f$7 s and $(i\omega I-A)\hat x=\hat f$8 GB, whereas the randomized algorithm with $(i\omega I-A)\hat x=\hat f$9 took $R(\omega)=(i\omega I-A)^{-1},$0 s and $R(\omega)=(i\omega I-A)^{-1},$1 GB. With a physics-informed sketch $R(\omega)=(i\omega I-A)^{-1},$2, the sketch size could be reduced from $R(\omega)=(i\omega I-A)^{-1},$3 to $R(\omega)=(i\omega I-A)^{-1},$4 with no loss in accuracy in that example (Ribeiro et al., 2019).

Barthel, Gomez, and McKeon proposed a variational formulation in which resolvent response modes are defined as stationary points of

$R(\omega)=(i\omega I-A)^{-1},$5

subject to $R(\omega)=(i\omega I-A)^{-1},$6. This converts the problem into a generalized eigenvalue problem for

$R(\omega)=(i\omega I-A)^{-1},$7

avoids explicit matrix inversion, and permits approximation in any chosen basis. In the examples reported there, the method achieved model reduction of up to two orders of magnitude, with $R(\omega)=(i\omega I-A)^{-1},$8 reduction in compute time and $R(\omega)=(i\omega I-A)^{-1},$9 memory savings for a streamwise-developing turbulent boundary layer (Barthel et al., 2021).

A different matrix-free strategy was introduced by Martini et al. through time-domain integration of the direct and adjoint linearized equations. Their transient-response and steady-state periodic-forcing methods compute all frequencies of interest simultaneously and provide an order-of-magnitude speedup when compared to previous matrix-free time stepping methods. The steady-state periodic-forcing method also recovers suboptimal modes for a discrete set of frequencies, rather than only the optimal mode (Martini et al., 2020).

4. Data-driven and time-localized generalizations

A major development is equation-free RA. In “Data-driven resolvent analysis,” the operator is replaced by an eigendecomposition learned from transient data via dynamic mode decomposition (DMD). The procedure collects unforced trajectories initialized with distinct perturbations, computes exact DMD to approximate eigenpairs of the underlying operator, projects the resolvent onto the learned eigenbasis, performs a weighted SVD in the reduced space, and lifts the forcing and response modes back to physical space. In reduced coordinates, the dynamics are

$H(\omega)=C(i\omega I-A)^{-1}$0

and the reduced weighted SVD is performed at cost $H(\omega)=C(i\omega I-A)^{-1}$1. The method is equation-free and adjoint-free, but it assumes a linearized time-invariant and asymptotically stable flow, requires sufficiently rich transient data, and depends on accurate DMD eigenvalues and biorthogonality of direct and adjoint modes (Herrmann et al., 2020).

Classical RA is tied to a Fourier decomposition in time, and therefore to statistically stationary or time-periodic settings. Two recent extensions relax that restriction. Wavelet-based resolvent analysis replaces the Fourier transform with a continuous wavelet transform, defines a time–frequency-localized operator $H(\omega)=C(i\omega I-A)^{-1}$2, and computes an SVD at each scale–time pair $H(\omega)=C(i\omega I-A)^{-1}$3. For statistically stationary channel flow at $H(\omega)=C(i\omega I-A)^{-1}$4, the wavelet and Fourier formulations agree to better than $H(\omega)=C(i\omega I-A)^{-1}$5, while in oscillating and abruptly forced shear flows the wavelet modes track the evolving mean profile and isolate transient growth mechanisms such as the Orr mechanism (Ballouz et al., 2024).

Sparse space-time resolvent analysis incorporates the temporal dimension through a discrete time-differentiation operator and promotes localization with an $H(\omega)=C(i\omega I-A)^{-1}$6-norm penalization. The resulting nonlinear eigenproblem is solved by an inverse power method. In statistically stationary turbulent channel flow, the sparse variant isolates a single vortex–streak unit cell rather than a periodic train; in time-varying systems, including a turbulent Stokes boundary layer and a channel flow subjected to a sudden lateral pressure gradient, it identifies temporally localized amplification events not directly evident from standard resolvent analysis (Lopez-Doriga et al., 2024).

5. Applications in aerodynamics, wall turbulence, reacting flow, and biofluid mechanics

RA has been used directly for flow-control design. In the NACA 0012 study at chord-based Reynolds number $H(\omega)=C(i\omega I-A)^{-1}$7 and Mach number $H(\omega)=C(i\omega I-A)^{-1}$8, resolvent gains peaked at chord-based Strouhal numbers $H(\omega)=C(i\omega I-A)^{-1}$9–$y=Cx$0 and small $y=Cx$1, indicating effective actuation parameters for open-loop separation control. Large-eddy simulations of approximately $y=Cx$2 controlled cases showed drag reduction up to $y=Cx$3 and lift increase up to $y=Cx$4, and the modal metric

$y=Cx$5

correlated strongly with the measured control performance (Yeh et al., 2018).

In wall-bounded turbulence, RA has been applied to smooth walls, riblets, ducts, pipes, and stratified channels. For riblets, a linear penalization term $y=Cx$6 extends the governing equations, producing a block-coupled resolvent operator that predicts both deterioration in performance with increasing riblet size and the emergence of spanwise rollers resembling Kelvin–Helmholtz vortices. The limited riblet shape optimization reported there yielded $y=Cx$7. In rectangular ducts, RA showed that finite-span geometry alters both gains and mode shapes and that secondary flow can either enhance or suppress amplification depending on the regime. In stratified channel flow, leading resolvent modes reproduced the balance of energy budget terms and coherent-structure trends across friction Richardson numbers $y=Cx$8 using only mean profiles and a limited range of representative scales (Chavarin et al., 2018, Lopez-Doriga et al., 2022, Ahmed et al., 2021).

Reacting and thermally coupled flows have also become a significant application domain. In a bluff-body stabilized flame with conjugate heat transfer, RA of the monolithically coupled fluid–solid Jacobian revealed that heat fluctuations are maximized when heat transfer between the fluid and solid is minimized, with anchored-flame states exhibiting strong low-$y=Cx$9 amplification and flashback states showing overall low amplification. In a turbulent hydrogen-air slot flame at Reynolds number $\partial_t q' = L[\bar q]\,q' + u',$0, Karlovitz number $\partial_t q' = L[\bar q]\,q' + u',$1, and equivalence ratio $\partial_t q' = L[\bar q]\,q' + u',$2, both shifted SPOD and RA identified Kelvin–Helmholtz wave packets over $\partial_t q' = L[\bar q]\,q' + u',$3–$\partial_t q' = L[\bar q]\,q' + u',$4 Hz; the velocity fluctuation mode shapes agreed well, while progress-variable and heat-release agreement improved when a generalized active-flame closure calibrated with high-fidelity data replaced the classical linearized RANS–EBU closure (Chen et al., 2024, Talasikar et al., 29 Mar 2026).

In stenotic transitional flow at $\partial_t q' = L[\bar q]\,q' + u',$5, RA about the LES mean identified a low-frequency sinuous stationary mode and an intermediate-frequency amplification region in which the most amplified fluctuations are axisymmetric. At intermediate frequencies, the optimal response mode showed both high gain separation and strong alignment with the leading SPOD mode. The same low-rank structure was then used to reconstruct turbulent kinetic energy and turbulent wall shear stress from the optimal response mode, especially in the immediate post-stenotic zone where axisymmetric fluctuations dominate the wall-shear signal (Villié et al., 14 Jun 2026).

6. Limits, misconceptions, and active directions

RA is often described as a linear theory, but its practical use in turbulence depends on a precise division of labor between linear amplification and nonlinear forcing. The cylinder-wake study at $\partial_t q' = L[\bar q]\,q' + u',$6 makes the limitation explicit: although resolvent analysis achieves a suitable energy balance for each harmonic, it does not correctly model energy transfer between temporal frequencies. DNS showed that nonlinearity neither produces nor consumes energy overall, yet it redistributes energy through both a forward cascade from low frequencies to high frequencies and a considerable inverse cascade from high frequencies to low frequencies. Higher harmonics therefore draw almost all their energy from nonlinear transfer rather than direct linear amplification of the mean (Jin et al., 2020).

This limitation is closely related to a common misconception: RA does not, by itself, supply a closure for nonlinear transfer. The same cylinder-wake analysis points to remedies that require additional information, including extended resolvent operators that include selected harmonic base flows or additional mean–harmonic couplings, harmonic resolvent formulations based on a time-periodic base flow, low-order self-consistent models coupling the mean and fundamental harmonic, and data-driven or physics-informed eddy-viscosity models calibrated to reproduce DNS triadic transfer (Jin et al., 2020).

A second misconception is that RA is confined to stable, equation-based, statistically stationary settings. Classical formulations do have such assumptions, and unstable operators require discounting or temporal filtering. In the airfoil-control formulation, a discounting factor $\partial_t q' = L[\bar q]\,q' + u',$7 with $\partial_t q' = L[\bar q]\,q' + u',$8 larger than the maximal growth rate shifts the spectrum into the stable half-plane and yields a bounded gain over the finite-time horizon $\partial_t q' = L[\bar q]\,q' + u',$9. In the data-driven setting, asymptotic stability is likewise required so that DMD captures decaying transients. Yet the broader literature now includes randomized, inverse-free, equation-free, wavelet-based, and sparse space-time formulations, each relaxing a different bottleneck—operator inversion, explicit adjoints, governing-equation access, or temporal stationarity (Yeh et al., 2018, Herrmann et al., 2020, Ballouz et al., 2024, Lopez-Doriga et al., 2024).

A plausible implication is that contemporary RA is best understood not as a single algorithm, but as a family of input–output formulations centered on the same spectral object—the resolvent operator—and adapted to distinct constraints on data, scale, non-stationarity, and physics. Across these variants, the enduring invariant is the decomposition of a linearized dynamics into forcing modes, response modes, and gains, with the quality of the resulting model governed by norm choice, mean-flow fidelity, gain separation, and the adequacy of the assumed forcing representation (Rolandi et al., 2024, Barthel et al., 2021).