# Efficient harmonic resolvent analysis via time-stepping

Ali Farghadan Department of Mechanical Engineering, University of Michigan,
Ann Arbor, MI, USA  Junoh Jung Department of Mechanical Engineering,
University of Michigan, Ann Arbor, MI, USA  Rutvij Bhagwat Department of
Mechanical Engineering, University of Michigan, Ann Arbor, MI, USA  Aaron
Towne Email address for correspondence<EMAIL_ADDRESS>Department of
Mechanical Engineering, University of Michigan, Ann Arbor, MI, USA

###### Abstract

We present an extension of the RSVD-$\Delta t$ algorithm initially developed
for resolvent analysis of statistically stationary flows to handle harmonic
resolvent analysis of time-periodic flows. The harmonic resolvent operator, as
proposed by Padovan et al. (2020), characterizes the linearized dynamics of
time-periodic flows in the frequency domain, and its singular value
decomposition reveals forcing and response modes with optimal energetic gain.
However, computing harmonic resolvent modes poses challenges due to $(i)$ the
coupling of all $N_{\omega}$ retained frequencies into a single harmonic
resolvent operator and $(ii)$ the singularity or near-singularity of the
operator, making harmonic resolvent analysis considerably more computationally
expensive than a standard resolvent analysis. To overcome these challenges,
the RSVD-$\Delta t$ algorithm leverages time stepping of the underlying time-
periodic linearized Navier-Stokes operator, which is $N_{\omega}$ times
smaller than the harmonic resolvent operator, to compute the action of the
harmonic resolvent operator. We develop strategies to minimize the algorithm’s
CPU and memory consumption, and our results demonstrate that these costs scale
linearly with the problem dimension. We validate the RSVD-$\Delta t$ algorithm
by computing modes for a periodically varying Ginzburg-Landau equation and
demonstrate its performance using the flow over an airfoil.

##  1 Introduction

Model reduction plays a critical role in the study of fluid mechanics due to
the complex and high-dimensional nature of fluid flows, especially turbulent
ones. In particular, modal decompositions, both data-driven and equation-
based, have proven to be effective at identifying low-dimensional sets of
modes associated with interpretable coherent structures that significantly
influence quantities of engineering interest such as kinetic energy, heat and
mass transfer, and noise emissions (Taira et al., 2017; Towne et al., 2018).
In particular, resolvent analysis has been extensively employed to comprehend,
model, and control statistically stationary turbulent flows such as wall-
bounded turbulence (Dawson & McKeon, 2019), turbulent boundary layer (Sipp &
Marquet, 2013), and jet flows (Pickering et al., 2021). However, real-world
fluid systems often exhibit non-stationary behavior, including periodicity,
rendering resolvent analysis less effective. Examples of flows exhibiting
periodic behavior include internal waves in a stratified fluid (Troy & Koseff,
2005), vortex shedding in the wake of a bluff body (Giannenas et al., 2022),
and pulsatile blood flow (Farghadan & Arzani, 2019). Recently, an extension of
resolvent analysis, known as harmonic resolvent analysis, has been developed
by Padovan et al. (2020). This extension characterizes the perturbation
forcing and responses around a periodic base flow, providing a more effective
tool for analyzing, modeling, and controlling periodically time-varying flows.
The harmonic resolvent operator can be considered a special case of the
harmonic transfer function introduced by Wereley (1990), wherein the operator
exhibits a dominant frequency.

Investigating the harmonic resolvent operator is one way to identify the
optimal spatio-temporal perturbations and their corresponding responses in
fluid systems. Unlike the resolvent operator, the harmonic resolvent operator
takes into consideration the periodicity of the base flow by representing it
in the form of a Fourier series. This allows for characterization of the
dynamic behavior of the system, particularly with respect to its spectral
components and it provides a means to understand and explore the interactions
between input and output modes through the linearized dynamics. Harmonic
resolvent analysis determines the triadic interactions between the input
frequency $\omega_{1}$, the base flow frequency $\omega_{1}-\omega_{2}$, and
the response frequency $\omega_{2}$. The number of Fourier modes of the base
flow is determined by its power spectral density (PSD), which in turn
determines the number of triads interacting with the forcing. Typically,
periodic flows with a dominant periodicity exhibit a multimodal nature with
their most energetic modes being harmonics of the base frequency. When the
base flow is statistically stationary, linearization around the temporal
average is sufficient, and traditional resolvent analysis arises as a special
case.

Recent studies have expanded on this research area to model and control the
dynamics of turbulent flows with periodic motions. Lin et al. (2023) applied
the harmonic resolvent framework for linear time-periodic systems with more
than one dominant frequency, offering a tool to study the flow control of a
plunging cylinder. In a distinct approach, Franceschini et al. (2022) proposed
a novel method aimed at capturing the phase dependency of small-scale
turbulent phenomena in flows exhibiting a periodic limit cycle. This method
builds upon classical spectral proper orthogonal decomposition (SPOD) and
resolvent analyses but introduces a quasi-steady approximation that
effectively separates the high-frequency turbulent component from the slow
periodic oscillation. Moreover, Heidt & Colonius (2023) have pioneered the
theory and algorithm for cyclostationary SPOD (CS-SPOD). This innovative
approach extends the conventional SPOD algorithm to effectively capture the
characteristics of flows characterized by periodic motions. CS-SPOD is closely
tied to harmonic resolvent analysis, akin to the relationship between SPOD and
resolvent analysis (Towne et al., 2018).

From a computational standpoint, computing harmonic resolvent modes presents
similar challenges to resolvent analysis, including the need to invert sparse
matrices in Fourier space. This challenge is further exacerbated in harmonic
resolvent analysis due to expanded system size, where all relevant frequencies
are merged into a unified operator, in contrast to resolvent analysis, where
each frequency has its own independent system. The second shared challenge
involves performing the singular value decomposition (SVD). Specifically,
turbulent flows characterized by periodic motions require a high-resolution
mesh, presenting both challenges as formidable obstacles. A modified
randomized SVD (RSVD) can be used to simultaneously overcome both
computational challenges (Moarref et al., 2013; Ribeiro et al., 2020; Padovan
et al., 2020). Nevertheless, the modified RSVD, referred to as “RSVD-LU” in
this study, becomes impractical for tackling large-scale turbulent flows since
it requires the lower-upper (LU) decomposition of the harmonic resolvent
operator. Recently, RSVD-$\Delta t$, an improvement to RSVD-LU, addresses
bottlenecks in resolvent analysis (Farghadan et al., 2023). The approach of
time-stepping, initially developed by Monokrousos et al. (2010) and further
optimized by Martini et al. (2021), allows for computing the action of an
operator on a vector without solving a linear system in Fourier space, making
it a key component of RSVD-$\Delta t$.

The main objective of this study is to expand the applicability of the
RSVD-$\Delta t$ algorithm, initially crafted for resolvent analysis, to
facilitate the computation of harmonic resolvent modes. A streaming approach
is adopted for time integration, reducing memory usage while maintaining
computational efficiency. Moreover, an efficient transient removal strategy is
introduced, cutting the temporal length of the integration and thereby
minimizing the CPU cost of RSVD-$\Delta t$ for achieving a desired accuracy.
The RSVD-$\Delta t$ algorithm for periodic flows maintains linear scalability
in both CPU time and memory consumption. This scalability enables the
efficient computation of harmonic resolvent modes for large-scale turbulent
flows that were previously out of reach.

The structure of our paper unfolds as follows. We begin by explaining
resolvent analysis in $\S$2, followed by the formulation of harmonic resolvent
analysis in $\S$3, where the computation of harmonic resolvent modes using
RSVD-LU is briefly described. The RSVD-$\Delta t$ algorithm is introduced in
$\S$4. A comparison of the computational complexity of algorithms and
suggested improvements appear in $\S$5. Further details on performance and
error analysis are presented in $\S$6. In $\S$7, we conduct two test cases to
demonstrate the accuracy and cost efficiency of RSVD-$\Delta t$ for two
periodic systems, and we conclude with final remarks in $\S$8.

##  2 Resolvent analysis

Before considering the time-periodic case, we briefly review standard
resolvent analysis for stationary flows. The Navier-Stokes equations can be
linearized around a base flow via Reynolds decomposition of the state,
$\bm{q}=\bar{\bm{q}}+\bm{q}^{\prime}$, where $\bar{\bm{q}}$ represents the
temporal average of the flow, and $\bm{q}^{\prime}$ denotes the time-varying
fluctuations. The dynamics of perturbations are then described by the linear
time-invariant system

$\frac{d\bm{q}^{\prime}}{dt}=\bm{A}(\bar{\bm{q}})\bm{q}^{\prime}+\bm{f}^{\prime}(\bar{\bm{q}},\bm{q}^{\prime}),$
(2.1)

where $\bm{A}\in\mathbb{C}^{N\times N}$ is the linearized Navier-Stokes (LNS)
operator, $N$ is the state dimension of the discretized system, and
$\bm{f}^{\prime}$ represents an external forcing as well as nonlinear terms
that are treated as an additional forcing, as suggested by McKeon & Sharma
(2010). It should be noted that $\bm{A}$ is time-independent because the base
flow is time-invariant. By taking the Fourier transform

$\mathcal{F}(\cdot)=\hat{(\cdot)}(\omega)=\int_{-\infty}^{+\infty}(\cdot)e^{-i\omega
t}\,dt$ (2.2)

in time, we obtain

$\hat{\bm{q}}(\omega)=\bm{R}(\omega)\hat{\bm{f}}(\omega),$ (2.3)

where $\omega$ represents the frequency and $\hat{(\cdot)}$ denotes the
frequency-domain representation of the time-domain vector. The resolvent
operator

$\bm{R}=(\text{i}\omega\bm{I}-\bm{A})^{-1}\in\mathbb{C}^{N\times N}$ (2.4)

is a transfer function that maps the input forcing to the output response.
Here, $\text{i}=\sqrt{-1}$ and $\bm{I}$ is the identity matrix of appropriate
dimension.

To analyze the perturbation characteristics of a system, the transfer function
can be examined across a range of frequencies. Typically, the most amplified
responses are identified by computing the left singular vectors
$\textbf{U}_{R}$ of the resolvent operator,

$\bm{R}=\textbf{U}_{R}\bm{\Sigma}_{R}\textbf{V}_{R}^{*},$ (2.5)

where $(\cdot)^{*}$ denotes complex conjugate transpose. The optimal forcing
inputs that lead to the highest amplifications comprise the right singular
vectors $\textbf{V}_{R}$, and the degree of amplifications are determined by
the singular values $\bm{\Sigma}_{R}$.

The SVD in 2.5 implies that disturbance amplitudes are measured in a Euclidean
norm. More generally,
$||\bm{x}||^{2}_{f}=\langle\bm{x},\bm{x}\rangle_{f}=\bm{x}^{*}\bm{W}_{f}\bm{x}$
computes the $f$-norm of any vector $\bm{x}$, and input and output norms can
be different, _i.e._ , $||\cdot||_{q}=||\cdot||_{f}$ is not necessary.
Additionally, one can define an input matrix $\bm{B}$ to restrict the forcing
and an output matrix $\bm{C}$ to extract the output of interest from the
state. For the sake of notational brevity, we assumed the weight, input, and
output matrices to be the identity, _i.e._ , $\bm{W}=\bm{B}=\bm{C}=\bm{I}$,
throughout this paper. However, our algorithm, with minor adjustments, can
accommodate non-identity weight, input, and output matrices by computing the
modes of the modified resolvent operator

$\tilde{\bm{R}}=\bm{W}^{1/2}_{q}\bm{C}\bm{R}\bm{B}\bm{W}^{-1/2}_{f}.$ (2.6)

For more details, please refer to Farghadan et al. (2023).

##  3 Harmonic resolvent analysis

In this section, we provide an overview of harmonic resolvent analysis. After
reviewing its basic formulation, we discuss the singularity issue that can
arise in some special flows, the truncation of the harmonic resolvent
operator, and two variations of harmonic resolvent analysis, namely cross-
frequency amplification and subharmonic resolvent. Lastly, a brief overview is
provided for the state-of-the-art algorithm used to compute the harmonic
resolvent modes and gains.

###  3.1 Formulation

Analogous to resolvent analysis, by inserting the generalized Reynolds
decomposition $\bm{q}(t)=\bar{\bm{q}}(t)+\bm{q}^{\prime}(t)$ into the Navier-
Stokes equations, the linearization around a periodic base flow
$\bar{\bm{q}}(t)=\bar{\bm{q}}(t+T)$ leads to the linear time-periodic system

$\frac{d\bm{q}^{\prime}}{dt}=\bm{A}_{p}\bm{q}^{\prime}+\bm{f},$ (3.1)

where $\bm{A}_{p}(t)=\bm{A}_{p}(t+T)\in\mathbb{C}^{N\times N}$ is the periodic
LNS operator, $\omega_{f}$ is the fundamental frequency, and
$T=2\pi/\omega_{f}$ is defined as the fundamental cycle length. As both
$\bm{A}_{p}$ and $\bm{q}^{\prime}$ are time-periodic, the Fourier transform of
(3.1) incorporates interactions between all pairs of frequencies. For
simplicity, we drop the prime notation for fluctuations from here on out.

Expanding $\bm{A}_{p}(t)$ and $\bm{q}(t)$ in terms of the Fourier series

$\begin{gathered}\bm{A}_{p}(t)=\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t},\\\
\bm{q}(t)=\sum_{j=-\infty}^{\infty}\hat{\bm{q}}_{j}e^{\text{i}j\omega_{f}t},\end{gathered}$
(3.2)

where $(\cdot)_{j}$ denotes $j^{th}$ harmonic of $\omega_{f}$, and
substituting these expansions into (3.1), we obtain

$\left[\bm{T}\hat{\bm{q}}\right]_{k}=\text{i}k\omega_{f}\hat{\bm{q}}_{k}-\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{k-j}\hat{\bm{q}}_{j}=\hat{\bm{f}}_{k}.$
(3.3)

Here, $\bm{T}$ is an infinite dimensional block matrix of the form

$\bm{T}=\begin{bmatrix}\ddots&\vdots&\vdots&\vdots&\vdots&\udots\\\
\dots&\bm{R}_{-1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&\hat{\bm{A}}_{-3}&\dots\\\
\dots&\hat{\bm{A}}_{1}&\bm{R}_{0}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&\dots\\\
\dots&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{1}^{-1}&\hat{\bm{A}}_{-1}&\dots\\\
\dots&\hat{\bm{A}}_{3}&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{2}^{-1}&\dots\\\
\udots&\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},$ (3.4)

where the diagonal entries are in fact the inverse of resolvent operators at
various frequencies,

$\bm{R}_{j}^{-1}=\text{i}j\omega_{f}\bm{I}-\hat{\bm{A}}_{p,0},\forall
j\in\mathbb{Z}.$ (3.5)

The harmonic resolvent operator

$\bm{H}=\bm{T}^{-1}$ (3.6)

transfers the harmonic inputs to the outputs for periodic base flows,

$\hat{\bm{q}}=\bm{H}\hat{\bm{f}}.$ (3.7)

The SVD of the harmonic resolvent operator

$\bm{H}=\textbf{U}_{H}\bm{\Sigma}_{H}\textbf{V}_{H}^{*}$ (3.8)

unveils the most amplified responses $\textbf{U}_{H}$ that correspond to
optimal forcing $\textbf{V}_{H}$, each accompanied by associated amplification
magnitudes $\bm{\Sigma}_{H}$. The matrix $\bm{H}$ encompasses the base flow,
the fundamental frequency, and higher harmonics. Whereas individual
frequencies can be analyzed separately in a standard resolvent analysis, all
frequencies are coupled within $\bm{H}$ in harmonic resolvent analysis and
must be considered simultaneously. Furthermore, the modes of the modified
harmonic resolvent operator, defined as

$\tilde{\bm{H}}=\bm{W}^{1/2}_{q}\bm{C}\bm{H}\bm{B}\bm{W}^{-1/2}_{f},$ (3.9)

can be computed using our algorithm with the same minor adjustments as
detailed in Farghadan et al. (2023). In a special case where the base flow is
time-independent, _i.e._ , $\hat{\bm{A}}_{p,j}=0$ for $j\neq 0$, $\bm{H}$
reduces to

$\bm{H}=\begin{bmatrix}\ddots&\vdots&\vdots&\vdots&\udots\\\
\dots&\bm{R}_{-1}&0&0&\dots\\\ \dots&0&\bm{R}_{0}&0&\dots\\\
\dots&0&0&\bm{R}_{1}&\dots\\\
\udots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},$ (3.10)

so each frequency is decoupled. Therefore, resolvent analysis can be perceived
as a special case of harmonic resolvent analysis.

###  3.2 Singular or near-singular nature of $\bm{T}$

The operator $\bm{T}$ tends to be ill-conditioned for many periodic systems
(Padovan et al., 2020). A system is considered ill-conditioned when the ratio
of the largest and smallest singular values is large
$(\kappa=\sigma_{\text{max}}/\sigma_{\text{min}}\sim O(10^{4})$ or higher).
The smallest singular value of $T$ can be shown to reach machine precision
zero if the periodic base flow satisfies the Navier-Stokes equations (Padovan
et al., 2020; Leclercq & Sipp, 2023), _i.e._ , when

$\frac{d\bm{\bar{q}}}{dt}=\bm{\mathcal{N}}(\bm{\bar{q}}),$ (3.11)

where $\bm{\mathcal{N}}$ is the nonlinear Navier-Stokes operator. Taking a
time derivative of both sides,

$\frac{d\left(\frac{d\bm{\bar{q}}}{dt}\right)}{dt}=\frac{d\left(\bm{\mathcal{N}}(\bm{\bar{q}})\right)}{dt}=\frac{\partial\bm{\mathcal{N}}}{\partial\bm{\bar{q}}}\frac{d\bm{\bar{q}}}{dt}=\bm{A}_{p}\frac{d\bm{\bar{q}}}{dt}.$
(3.12)

Taking a Fourier transform of (3.12) and recalling that $\bm{T}$ is obtained
as shown in (3.3), $\widehat{d\bm{\bar{q}}/dt}$ is a non-trivial solution of

$\bm{T}\bm{q}=0,$ (3.13)

residing within the null space of $\bm{T}$, thereby proving the singularity of
$\bm{T}$. In systems where the base flow approximately satisfies (3.11),
$\bm{T}$ becomes nearly singular. Consequently, irrespective of the specific
periodic system under consideration, $\bm{T}$ might be poorly conditioned,
presenting even more challenges when dealing with large systems.

Padovan et al. (2020) demonstrated that defining the harmonic resolvent
operator in a manner that eliminates the phase shift imparted by the Fourier
coefficients of $d\bm{\bar{q}}/dt$ is effective, without adversely affecting
the other dominant amplification mechanisms. By constraining $\bm{T}$ to a
subspace $\bm{U}^{\perp}$ that is orthogonal to the direction of the phase
shift $\widehat{d\bm{\bar{q}}/dt}$, the range of the harmonic resolvent
operator is limited to $\bm{U}^{\perp}$. The key concept here is to eliminate
the singular vectors associated with the smallest singular value of $\bm{T}$
without affecting the other singular values and vectors of $\bm{T}$. Given
that the phase shift resides in the null space of $\bm{T}$, the desired right
singular vector
$\bm{v}=\widehat{\frac{d\bm{\bar{q}}}{dt}}/||\widehat{\frac{d\bm{\bar{q}}}{dt}}||$
is readily available. We can also solve

$\bm{T}^{*}\bm{\tilde{u}}=\bm{v}$ (3.14)

and compute the corresponding left singular vector as
$\bm{u}=\bm{\tilde{u}}/||\bm{\tilde{u}}||$. Upon obtaining $\bm{v}$ and
$\bm{u}$, the process becomes straightforward: projecting out these modes from
the forcing and response terms effectively restricts $\bm{T}$ to
$\bm{U}^{\perp}$. These steps are elucidated in detail in Padovan et al.
(2020).

To compute the action of the harmonic resolvent operator on a given vector (or
matrix), our goal is to determine the vector (or matrix) $\bm{q}$ that
satisfies $\bm{q}=\bm{H}\bm{f}$. In brief, to compute the action of a singular
$\bm{H}$ on a vector $\tilde{\bm{f}}$, the component along $\bm{u}$ must be
projected out as

$\tilde{\bm{f}}_{in}=\tilde{\bm{f}}-\bm{u}(\bm{u}^{*}\tilde{\bm{f}}),$ (3.15)

and this will provide a response

$\bm{T}\bm{\tilde{q}}=\tilde{\bm{f}}_{in}\rightarrow\bm{\tilde{q}}=\bm{H}\tilde{\bm{f}}_{in},$
(3.16)

which is orthogonal to the direction of the phase shift. Nevertheless, to
mitigate potential round-off errors, it is advisable to refine the output by
projecting out the component along $\bm{v}$ as

$\tilde{\bm{q}}_{out}=\bm{\tilde{q}}-\bm{v}(\bm{v}^{*}\bm{\tilde{q}}).$ (3.17)

Similarly, when computing the action of singular $\bm{H}^{*}$ on a vector, the
initial step involves projecting out the component along $\bm{v}$. After the
response is obtained, to mitigate round-off errors, the component along
$\bm{u}$ is projected out.

###  3.3 Truncating the harmonic resolvent operator

Computing harmonic resolvent modes is contingent upon $\bm{T}$ being of finite
size. The size and block sparsity pattern of $\bm{T}$ are determined by two
key variables: the number of frequencies within the base flow, $N_{b}$, and
within the response, $N_{\omega}$. Here, we elucidate the impact of each
variable on the operator.

The rows of $\bm{T}$ consist of $N_{b}$-stencil blocks, dictated by the
presence of non-zero elements in $\hat{\bm{A}}_{p,j}$. In the context of
periodic flows, the time-dependent base flow can be expressed using a Fourier
series expansion

$\bar{\bm{q}}(t)=\sum_{j\omega_{f}\in\Omega_{\bar{q}}}\hat{\bar{\bm{q}}}_{j}e^{\text{i}j\omega_{f}t}.$
(3.18)

The set $\Omega_{\bar{q}}$ encompasses the relevant frequencies associated
with the dominant flow structures. Typically, $\Omega_{\bar{q}}$ is a subset
of $j\omega_{f}$ with $j\in\mathbb{Z}$ and encapsulates the majority of energy
within the periodic flows. By retaining the base periodicity and a limited
number of higher harmonics, it becomes possible to simplify the base flow
representation, resulting in $N_{b}$ frequencies within $\Omega_{\bar{q}}$ and
$N_{b}$ non-zero elements within each row of $\bm{T}$.

On the other hand, we aim to resolve perturbations occurring at temporal
frequencies $\omega\in\Omega_{q}$, where $\Omega_{q}$ represents a subset of
$j\omega_{f}$ with $j\in\mathbb{Z}$, although the specific choice of harmonics
for $\Omega_{q}$ may vary, as further elucidated in $\S$3.5. Typically,
$\Omega_{\bar{q}}$ is a subset of $\Omega_{q}$, _i.e._ , the perturbation
frequency content spans a wider range of harmonics (Padovan et al., 2020). The
number of frequencies $N_{\omega}$ expanding both $\bm{f}$ and $\bm{q}$
determines the dimensions of the block matrices, yielding a size of
$(N_{\omega}N)\times(N_{\omega}N)$ for $\bm{T}$.

We provide a simple example to visually observe the harmonic resolvent
operator. Suppose we use $\Omega_{\bar{q}}=\\{-2,-1,0,1,2\\}\omega_{f}$ to
expand the LNS operator and $\Omega_{q}=\\{-4,-3,\dots,3,4\\}\omega_{f}$ to
expand the perturbations. In this case, the stencil length would be $N_{b}=5$,
and the size of the block matrix would be $(9N)\times(9N)$. The harmonic
resolvent operator is then

$\bm{H}=\begin{bmatrix}\bm{R}_{-4}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0&0&0\\\
\hat{\bm{A}}_{1}&\bm{R}_{-3}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0&0\\\
\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{-2}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0&0\\\
0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{-1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0&0\\\
0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{0}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0&0\\\
0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{1}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}&0\\\
0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{2}^{-1}&\hat{\bm{A}}_{-1}&\hat{\bm{A}}_{-2}\\\
0&0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{3}^{-1}&\hat{\bm{A}}_{-1}\\\
0&0&0&0&0&0&\hat{\bm{A}}_{2}&\hat{\bm{A}}_{1}&\bm{R}_{4}^{-1}\end{bmatrix}^{-1}.$
(3.19)

The forcing and response vectors can be represented as

$\hat{\bm{f}}=\begin{bmatrix}\hat{\bm{f}}_{-4}\\\ \hat{\bm{f}}_{-3}\\\
\hat{\bm{f}}_{-2}\\\ \hat{\bm{f}}_{-1}\\\ \hat{\bm{f}}_{0}\\\
\hat{\bm{f}}_{1}\\\ \hat{\bm{f}}_{2}\\\ \hat{\bm{f}}_{3}\\\
\hat{\bm{f}}_{4}\\\
\end{bmatrix},\hat{\bm{q}}=\begin{bmatrix}\hat{\bm{q}}_{-4}\\\
\hat{\bm{q}}_{-3}\\\ \hat{\bm{q}}_{-2}\\\ \hat{\bm{q}}_{-1}\\\
\hat{\bm{q}}_{0}\\\ \hat{\bm{q}}_{1}\\\ \hat{\bm{q}}_{2}\\\
\hat{\bm{q}}_{3}\\\ \hat{\bm{q}}_{4}\\\ \end{bmatrix}.$ (3.20)

The accuracy of the obtained resolvent system in Fourier space depends on how
well the chosen frequencies capture the relevant flow information, and it is
affected by the truncation limit imposed on the frequency range. Although it
is feasible to reduce the infinite-dimensional harmonic operator to a finite
dimension based on perturbation frequencies, the choice of frequency range is
consequential. To demonstrate the influence of truncation on response modes,
consider a periodic flow with a dominant first harmonic, such that the base
flow can be expressed as

$\bar{\bm{q}}(t)=\hat{\bar{\bm{q}}}_{0}+\hat{\bar{\bm{q}}}_{1}e^{\text{i}\omega
t}+\hat{\bar{\bm{q}}}_{-1}e^{-\text{i}\omega t}.$ (3.21)

Assume that the LNS equations are subjected to a constant forcing
$\bm{f}(t)=\hat{\bm{f}}_{0}$. The LNS operator exhibits frequency content at 0
and $\pm\omega$, indicating that a constant forcing will directly induce a
response at these frequencies. However, higher harmonics are also indirectly
stimulated, and the triggering of responses at higher harmonics follows a
cascade pattern, with norms gradually decaying towards zero, _i.e._ ,
$\lim_{j\to\infty}||\hat{\bm{q}}_{j}||=0$. Hence, the truncation must be set
at a sufficiently high level to preserve the response norms.

In practice, by defining $\Omega_{q}$ for both input perturbation and the
response, we compute accurate harmonic resolvent modes as long as the norms of
higher frequency modes become less relevant within the spectrum. In other
words, if we set both forcing and response frequencies sufficiently high and,
upon computing harmonic resolvent modes, find that the norms of input and
output modes with higher frequency content are relatively small, it indicates
convergence, and including higher frequencies has negligible impact on the
results. In a case, for instance, where response modes with higher frequency
content remain relatively important, extending the frequency range of the
output modes is necessary without changing the input frequency content. In our
algorithm, the input and output frequency ranges are adjustable as needed.

###  3.4 Cross-frequency harmonic resolvent analysis

Harmonic resolvent analysis captures the interaction between various
frequencies within $\Omega_{\bar{q}}$ and $\Omega_{q}$, and the typical
objective is to identify the optimal forcing mode, potentially spanning
various frequencies, that produces the optimal response across a range of
frequencies. Alternatively, one can compute the unit-norm forcing at one
frequency $\omega_{1}\in\Omega_{q}$ that triggers the most amplified response
at another frequency $\omega_{2}\in\Omega_{q}$. This is achieved by defining

$\bm{H}_{\omega_{2},\omega_{1}}=\bm{C}\bm{H}\bm{B}$ (3.22)

and computing the modes and gains following the same procedure as before.
Here, matrices $\bm{B}$ and $\bm{C}$ are constructed to extract the
$\omega_{1}$ and $\omega_{2}$ frequencies from the forcing and response modes,
respectively. This represents a special case of the modified harmonic
resolvent operator in (2.6).

###  3.5 Subharmonic resolvent analysis

The frequency content of the forcing inputs is typically composed of harmonics
of the base flow frequency. However, Padovan & Rowley (2022) have shown that
the harmonic resolvent system may exhibit sensitivity to subharmonic inputs.
For instance, in the case where the fundamental base flow frequency is denoted
by $\omega_{f}$, the linearized system may be sensitive to forcing inputs with
$\gamma<\omega_{f}$, such as $\omega_{f}/2$. More generally, if the frequency
content of a periodic LNS operator is limited to the frequency range
$\Omega_{\bar{q}}$, it may be desirable to compute input-output modes with
frequencies outside of this range, _i.e._ , , $\gamma\notin\Omega_{\bar{q}}$.
In such cases, the subharmonic resolvent modes become relevant.

To identify unique subharmonic inputs, we can take advantage of the linear
nature of the harmonic resolvent operator, which allows interaction with the
fundamental frequency $\omega_{f}$ and its harmonics. It can be demonstrated
that the interval $(-\omega_{f}/2,\omega_{f}/2]$ encompasses all sets of
subharmonic frequencies where the input-output modes are unique. As discussed
in $\S$3, a forcing with a frequency $\gamma\in(-\omega_{f}/2,\omega_{f}/2]$
can only trigger responses at frequencies $\gamma+j\omega_{f}$, where $j$ is
an integer. Therefore, the input and output modes exist within the frequency
range $\Omega_{\gamma}=\gamma\oplus\Omega_{\bar{q}}$, where $\oplus$ denotes
element-wise addition. Harmonics of $\gamma$ trigger a unique set of outputs
as long as $j\gamma\leq|\omega_{f}|$. Due to the linear relationship, we can
show that all harmonics are isolated sets that need to be studied separately.
Note that $-\omega_{f}/2$ is excluded since
$\Omega_{-\omega_{f}/2}=\Omega_{\omega_{f}/2}$, resulting in redundancy. For
$\gamma\in\Omega_{\gamma}$, one can derive the subharmonic resolvent system as
(Padovan & Rowley, 2022)

$\hat{\bm{q}}=(\text{i}\gamma\bm{I}-\bm{T})^{-1}\hat{\bm{f}}.$ (3.23)

To further elucidate the point, consider an example where the subharmonic
frequency of interest is $\omega_{f}/5$ along with its harmonics. Note that
the valid harmonics are the ones that fulfill the condition
$j\omega_{f}/5\in(-\omega_{f}/2,\omega_{f}/2]$, which includes the set
$\\{-2/5,-1/5,0,1/5,2/5\\}$. While
$\Omega_{0}=0\oplus\\{-N_{b}/2,\ldots,N_{b}/2\\}$, where $N_{b}$ is the number
of active frequencies within the base flow, is identical to harmonic inputs,
the sets $\Omega_{-2/5}$, $\Omega_{-1/5}$, $\Omega_{1/5}$, and $\Omega_{2/5}$
represent distinct input-output modes. In general, for any two distinct
subharmonic frequencies $\gamma_{1}$ and $\gamma_{2}$ within the interval
$(-\omega_{f}/2,\omega_{f}/2]$, the corresponding subharmonic modes in
$\Omega_{\gamma_{1}}$ are decoupled from those in $\Omega_{\gamma_{2}}$ as
there exists no common frequency $\gamma_{1}-\gamma_{2}\in\Omega_{\bar{q}}$.
While our algorithm described in $\S$4 is described for harmonic resolvent
analysis, it can be easily extended to compute subharmonic resolvent modes, as
described in Appendix A.

###  3.6 Computing harmonic resolvent modes using RSVD-LU

The RSVD algorithm (Halko et al., 2011) is a randomized linear algebra
technique developed to efficiently identify the singular vectors with the
highest gains in a given matrix. In the context of harmonic resolvent
analysis, the matrix of interest is the harmonic resolvent operator. This
operator, similar to the resolvent operator, relies on the inverse of a sparse
matrix, as shown in (3.6), making the computation of modes challenging.
Modifying the original RSVD algorithm is necessary to address these
challenges, enabling harmonic resolvent analysis. The modifications for
resolvent analysis have been extensively documented in the literature by
Ribeiro et al. (2020) and Farghadan et al. (2023), while the outline for the
harmonic resolvent operator can be found in the appendix of Padovan et al.
(2020).

In brief, the computation involves computing the action of $\bm{H}$ and
$\bm{H}^{*}$ for sketching the range and image of the harmonic resolvent
operator, respectively, and approximating the leading harmonic modes.
Computing actions of $\bm{H}$ and $\bm{H}^{*}$, however, requires solving
linear systems

$\begin{gathered}\bm{T}\hat{\bm{q}}=\hat{\bm{f}},\\\
\bm{T}^{*}\hat{\bm{q}}=\hat{\bm{f}},\end{gathered}$ (3.24)

respectively, in Fourier space. These linear systems are solved via LU-
decomposition of $\bm{T}$, but due to its large size and ill-conditioned
nature, computing its LU decomposition can be a formidable computational
obstacle, particularly for flows with three inhomogeneous directions. We
propose an alternative way in the following section.

##  4 Computing harmonic resolvent modes using time stepping

In this section, we show how time stepping can be used to efficiently compute
harmonic resolvent modes. The same concept was first used by Monokrousos et
al. (2010) in the context of the resolvent analysis and further optimized by
Martini et al. (2021) and Farghadan et al. (2023). Here, we introduce the
extension of our algorithm to compute the harmonic resolvent modes and gains.

###  4.1 Computing the action of $\bm{H}$ using time stepping

Figure 1: Flowchart depicting the action of $\bm{H}$ on a forcing comprised of
$N_{\omega}$ frequencies within the RSVD-LU (top route) and the RSVD-$\Delta
t$ (bottom route) algorithms. Both routes produce the same result, but the
bottom route is computationally advantageous for large systems.

Computing the action of $\bm{H}$ in the Fourier space poses a bottleneck
within the RSVD-LU algorithm. To overcome this limitation, we propose to use
time-stepping as an effective surrogate approach. Given that
$\Omega_{\bar{q}}$ represents the frequency content of the base flow around
which $\bm{T}$ is constructed, we define the time-domain linearized operator
as

$\bm{A}_{p}(t)=\sum_{j\omega_{f}\in\Omega_{\bar{q}}}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}.$
(4.1)

The action of $\bm{H}$ on $\hat{\bm{f}}$ in Fourier space is expressed in
(3.7), where the frequency content of the forcing lies within $\Omega_{q}$.

The steady-state solution of (3.1) is given by

$\hat{\bm{q}}_{s}=\bm{H}\hat{\bm{f}},$ (4.2)

when subjected to the forcing represented as

$\bm{f}=\sum_{j\omega_{f}\in\Omega_{q}}\hat{\bm{f}}_{j}e^{\text{i}j\omega_{f}t}.$
(4.3)

The steady-state response can be obtained through time-stepping before taking
a Fourier transform. The time-domain forcing is constructed to include
frequencies identical to those present in $\hat{\bm{f}}$. Since the time-
domain and frequency-domain frequencies are the same, both approaches yield
identical results, _i.e._ , $\hat{\bm{q}}_{s}=\hat{\bm{q}}$. In practice, the
time-stepping method is inherently discrete, and the frequencies are retained
up to the maximum limit determined by the chosen time step (below the Nyquist
frequency (Nyquist, 1928)). By employing a suitably small time step, we can
capture frequencies up to the desired limit.

The equivalence between computing the action of $\bm{H}$ in both the RSVD-LU
and RSVD-$\Delta t$ algorithms is depicted in figure 1. In the upper route,
within the RSVD-LU algorithm, the LNS equations are first transformed into
Fourier space before solving a coupled linear system in the frequency domain
for a given forcing. In the lower route, within the RSVD-$\Delta t$ algorithm,
the LNS equations are initially integrated in the time domain before being
transformed into Fourier space. Both routes produce identical outputs so long
as numerical artifacts are minimized.

###  4.2 RSVD-$\Delta t$ for harmonic resolvent analysis

By recognizing the equivalence between time stepping and solving linear
systems in the frequency domain, we can effectively address the limitations of
the RSVD-LU algorithm, leading to the development of RSVD-$\Delta t$. The
algorithmic steps for extending RSVD-$\Delta t$ to harmonic systems are
presented in Algorithm 1.

Algorithm 1 RSVD-$\Delta t$ for harmonic resolvent analysis

1:Input parameters: $\bm{A}_{p},k,q,\Omega_{q},\Omega_{\bar{q}},$TSS,
$dt,T_{t}$

2:$\hat{\bm{\varTheta}}\leftarrow\;$randn$(NN_{\omega},k)$ $\triangleright$
Create random test matrices

3:$\hat{\bm{Y}}\leftarrow\;$DirectAction$(\bm{A}_{p},\hat{\bm{\varTheta}},$TSS,
$dt,T_{t})$ $\triangleright$ Sample the range of $\bm{H}$

4:if $q>0$ then $\triangleright$ Optional power iteration

5: $\hat{\bm{Y}}\leftarrow\;$PI$(\bm{A}_{p},\hat{\bm{Y}},q,$TSS, $dt,T_{t})$
$\triangleright$ Power iteration with time-stepping

6:$\hat{\bm{Q}}_{\varOmega}\leftarrow\;$qr$(\hat{\bm{Y}}_{\varOmega})$
$\triangleright$ Build the orthonormal subspace $\hat{\bm{Q}}_{\varOmega}$

7:$\hat{\bm{S}}\leftarrow\;$ AdjointAction$(\bm{A}_{p}^{*},\hat{\bm{Q}},$TSS,
$dt,T_{t})$ $\triangleright$ Sample the image of $\bm{H}$

8:$(\tilde{\bm{U}}_{H},\bm{\varSigma}_{H},\bm{V}_{H})\leftarrow\;$svd$(\hat{\bm{S}}_{H})$
$\triangleright$ Obtain $\bm{\varSigma}_{H},\bm{V}_{H}$

9:$\bm{U}_{H}\leftarrow\;\hat{\bm{Q}}_{\varOmega}\tilde{\bm{U}}_{H}$
$\triangleright$ Recover $\bm{U}_{H}$

10:Output parameters: $\bm{U}_{H},\bm{\varSigma}_{H},\bm{V}_{H}$

Algorithm 1: $k,q,\Omega_{q},\Omega_{\bar{q}}$ are common parameters with
RSVD-LU. $(\cdot)_{\varOmega}$ indicates all frequencies are merged into a
single column, and TSS is an abbreviation for time-stepping schemes (_e.g._ ,
backward Euler). DirectAction and AdjointAction are functions that solve the
direct and adjoint LNS equations, respectively, with a predefined forcing.

The algorithm proposed in this study is based on the RSVD-$\Delta t$ algorithm
introduced by Farghadan et al. (2023), and we will provide a concise overview
of the key steps. Initially, a random matrix $\hat{\bm{\varTheta}}$ is
generated (line 2) to compute the sketch of $\bm{H}$ through time stepping
(line 3). DirectAction is a function that computes the action of $\bm{H}$ onto
a given forcing using time-stepping and transforms the steady-state responses
to Fourier space. Next, an optional power iteration is performed (lines 4 and
5) via $q$ successive applications of DirectAction and AdjointAction,
enhancing the accuracy of harmonic resolvent modes. QR decomposition is
performed to obtain $\hat{\bm{Q}}_{\varOmega}$ (line 6) before constructing
$\bm{S}$ via time stepping (line 7). AdjointAction is the function that
computes the action of $\bm{H}^{*}$ onto a given forcing using time-stepping
and transforms the steady-state responses to Fourier space. An SVD (line 8) is
conducted to obtain the optimal forcing modes
$\bm{V}_{H}\in\mathbb{C}^{NN_{\omega}\times k}$ and gains
$\bm{\varSigma}_{H}\in\mathbb{R}^{k\times k}$ of $\bm{H}$. Lastly, the optimal
response modes $\bm{U}_{H}\in\mathbb{C}^{NN_{\omega}\times k}$ are recovered
in line 9.

Compared to the RSVD-$\Delta t$ algorithm for resolvent analysis, two notable
differences are present. In the extension of RSVD-$\Delta t$, the response
modes at various frequencies are merged prior to both the QR decomposition and
the SVD steps. This deviation from the original algorithm is motivated by the
nature of the problem that requires interactions between frequencies, unlike
in resolvent analysis where the QR decomposition and SVD are performed
separately for each individual frequency of interest. The second difference is
the generation of LNS operators over time, whereas in resolvent analysis, the
LNS operator remains constant throughout the integration.

If $\bm{T}$ approaches machine precision singularity, it is necessary to
project out the forcing along the $\bm{u}$ direction before the DirectAction
and along $\bm{v}$ before the AdjointAction, as described in $\S$3.2.
Moreover, to mitigate numerical artifacts, it is advisable to project out the
response along $\bm{v}$ and $\bm{u}$ after the DirectAction and AdjointAction,
respectively.

##  5 Computational complexity

A brief comparison is conducted between memory requirements and CPU cost of
the RSVD-LU and RSVD-$\Delta t$ algorithms. We also propose strategies to
minimize the memory consumption. For more detailed information on the CPU and
memory scaling of the LU decomposition of the resolvent operator, we refer
readers to Farghadan et al. (2023). To offer empirical insights and
confirmation of the theoretical scalings discussed below, we present a two-
dimensional test case in $\S$7.

###  5.1 Memory benefits of RSVD-$\Delta t$

In the case of resolvent analysis, the memory scaling of the LU decomposition
is empirically $O(N^{1.2})$ and $O(N^{1.6})$ for two- and three-dimensional
systems, respectively (Towne et al., 2022; Farghadan et al., 2023). The memory
scaling of LU decomposition is expected to be worse for the harmonic resolvent
operator, as the lower- and upper-triangular matrices are denser. When
$\bm{T}$ is singular or nearly singular, the computational task becomes even
more challenging, as discussed in Padovan et al. (2020).

On the other hand, the memory usage of RSVD-$\Delta t$ for harmonic resolvent
analysis is contingent on the sizes of various matrices, specifically the
sparse LNS operators $\hat{\bm{A}}_{p}\in\mathbb{C}^{N\times N\times N_{b}}$,
and the dense forcing matrix $\hat{\bm{F}}\in\mathbb{C}^{NN_{\omega}\times k}$
and response matrix $\hat{\bm{Q}}\in\mathbb{C}^{NN_{\omega}\times k}$, all of
which are stored in Fourier space. The scaling for resolvent analysis is
$O(NN_{\omega})$, and the sole extra space is for the storage of LNS
operators.

One approach to generating LNS operators for all time points involves
retaining $N_{b}$ coefficients in Fourier space and subsequently constructing
$\bm{A}_{\varOmega}=\\{\bm{A}_{p,1},\bm{A}_{p,2},\bm{A}_{p,3},\dots,\bm{A}_{p,N_{t}}\\}$
once for all time steps within a period. However, this approach can be memory-
intensive for large systems, especially when $N_{t}\sim O(10^{3}-10^{5})$. An
alternative strategy to reduce memory consumption is to create LNS operators
on the fly, as shown in figure 2. As a result, no time-domain LNS operator is
permanently stored in memory. Each LNS operator is created on-the-fly using
the discrete Fourier transform (DFT) of $N_{b}$ Fourier coefficients,

$\bm{A}_{p}(t)=\sum_{j=-N_{b}/2}^{N_{b}/2}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}.$
(5.1)

Typically, $\bm{A}_{p}$ is very sparse, with $\texttt{nnz}(\bm{A}_{p})\sim
O(N)$, where $\texttt{nnz}(\cdot)$ represents the count of non-zero elements,
resulting in an overall memory consumption of $O(NN_{b})$. The methods
employed for generating LNS operators, both in the streaming and memory-
intensive approaches, mirror those used for creating forcing terms in
resolvent analysis, as elaborated in Farghadan et al. (2023). Additionally,
the CPU cost of this procedure scales as $O(\texttt{nnz}(\bm{A}_{p})N_{b})$ or
$O(NN_{b})$.

Figure 2: Schematic of the action of $\bm{H}$ with streaming DFT/iDFT to
transform between the Fourier and time domains.

Streaming DFT is indeed effective in reducing memory consumption;
nevertheless, additional memory savings can be specifically achieved for real-
valued LNS operators. In many cases, $\bm{A}_{p}$ is real-valued, _i.e._ ,
confined to $\mathbb{R}^{N\times N}$, yielding $\hat{\bm{A}}_{p,j}$ =
$\bar{\hat{\bm{A}}}_{p,-j}$, where $\bar{(\cdot)}$ denotes the complex
conjugate. Hence, $\hat{\bm{A}}_{p,j}$ with $j\geq 0$ suffices to construct
$\bm{A}_{p}$ as

$\bm{A}_{p}(t)=\hat{\bm{A}}_{p,0}+2\mathcal{R}e\left(\sum_{j=1}^{N_{b}/2}\hat{\bm{A}}_{p,j}e^{\text{i}j\omega_{f}t}\right),$
(5.2)

resulting in halving the number of Fourier coefficients of the LNS operator to
$\lfloor N_{b}/2\rfloor+1$. Here, $\mathcal{R}e$ represents the real part.

Another opportunity to reduce memory usage for real-valued LNS matrices arises
from the symmetry between positive and negative frequencies. Rewriting (3.3)
for $-k\omega_{f}$ yields

$\text{i}(-k\omega_{f})\hat{\bm{q}}_{-k}-\sum_{j=-\infty}^{\infty}\hat{\bm{A}}_{-k+j}\hat{\bm{q}}_{-j}=\hat{\bm{f}}_{-k},$
(5.3)

and the complex conjugate version of (3.3) becomes

$-\text{i}k\omega_{f}\bar{\hat{\bm{q}}}_{k}-\sum_{j=-\infty}^{\infty}\bar{\hat{\bm{A}}}_{k-j}\bar{\hat{\bm{q}}}_{j}=\bar{\hat{\bm{f}}}_{k}.$
(5.4)

Both equations (5.3) and (5.4) yield the same output for a given
$\bar{\hat{\bm{f}}}_{k}=\hat{\bm{f}}_{-k}$ since
$\bar{\hat{\bm{A}}}_{k-j}=\hat{\bm{A}}_{-k+j}$, inducing the harmonic
resolvent response and forcing modes containing
$\bm{U}_{H,j}=\bar{\bm{U}}_{H,-j}$ and $\bm{V}_{H,j}=\bar{\bm{V}}_{H,-j}$ for
$j\neq 0$, respectively. Hence, we only keep $\lfloor N_{\omega}/2\rfloor+1$
Fourier coefficients of the forcing and response modes, reducing memory
consumption by half when storing these dense matrices.

###  5.2 CPU benefits of RSVD-$\Delta t$

The computational cost of computing harmonic resolvent modes using RSVD-LU is
predominantly determined by the LU decomposition of $\bm{T}$. In the case of
resolvent analysis, the CPU scaling of the LU decomposition is $O(N^{1.5})$
and $O(N^{2})$ for two- and three-dimensional systems, respectively (Towne et
al., 2022; Farghadan et al., 2023). When, in general, block matrix $\bm{T}$
contains $\hat{\bm{A}}_{p,i}$ terms, and yields a denser matrix than block
diagonal with a complex sparsity pattern, the scaling is expected to be
$O((N_{\omega}N)^{1.5})$ or $O((N_{\omega}N)^{2})$ or worse.

The cost of using RSVD-$\Delta t$ is directly tied to the cost of computing
the actions of $\bm{H}$ and $\bm{H}^{*}$. In resolvent analysis, the total CPU
cost is proportional to the state dimension $N$. This cost can be broken down
into three main components: $(i)$ time-stepping, _e.g._ , Adams-Bashforth
methods, which scales as $O(N)$ when $\bm{A}$ is sparse, $(ii)$ creating
forcing from Fourier space to the time domain, which scales as
$O(NN_{t}N_{\omega})$, $(iii)$ taking the response from the time domain to
Fourier space, which scales as $O(NN_{\omega}^{2})$. Among these components,
the time-stepping part is the most significant contributor to the overall
cost. The CPU cost of harmonic resolvent analysis can be assessed similarly to
resolvent analysis, with one exception. For periodic flows, the LNS operator
varies over time, whereas it remains constant for statistically stationary
flows. We showed in $\S$5.1 that the cost of creating LNS operators over time
scales as $O(NN_{b})$ for sparse matrices.

In summary, the CPU and memory costs of RSVD-$\Delta t$ exhibit linear scaling
with respect to the state dimension $N$, regardless of whether the base flow
is steady or periodic. Moreover, the dominant terms in the CPU cost are
independent of the number of retained frequencies $N_{\omega}$. Both of these
properties make the RSVD-$\Delta t$ algorithm considerably more scalable than
RSVD-LU for harmonic resolvent analysis.

##  6 Minimizing the CPU cost of RSVD-$\Delta t$

Similar to time stepping in the context of resolvent analysis, periodic
systems may also experience a slow decay of transient response, which is
undesirable, since we require the steady-state response in isolation to
compute the action of the harmonic resolvent operator using time stepping.
This phenomenon might be linked to instances when $\bm{T}$ is near
singularity. Regardless, in cases where the LNS equations are absolutely
stable but exhibit slowly decaying modes, we propose a strategy that can
effectively reduce the CPU cost of time stepping.

###  6.1 Stability analysis: Floquet theorem and transient response

The steady-state response of the LNS equations, subject to forcing as
described in equation (3.1), is of interest for our analysis. In the case of a
linear time-invariant system where $\bm{A}$ is independent of time, the decay
rate is controlled by the least-damped eigenvalue of $\bm{A}$. This decay rate
can be estimated by integrating a homogeneous system over a sufficiently long
period. However, for a linear time-periodic system, the decay rate is
determined by the least-damped Floquet exponent.

Assuming $\bm{\Phi}(t)$ is the fundamental solution of (3.1), the monodromy
matrix is constructed as

$\bm{M}=\bm{\Phi}(0)^{-1}\bm{\Phi}(T),$ (6.1)

where $T$ is the period of the flow. When $\bm{\Phi}(0)=\bm{I},\bm{\Phi}(t)$
is referred to as the principal fundamental matrix, and the monodromy matrix
simplifies to $\bm{\Phi}(T)$. In other words, we evaluate the principal
fundamental matrix at the end of the first period. The eigenvalues of
$\bm{M}$, known as Floquet multipliers $\mu_{j}$, allow us to determine the
Floquet exponents $\lambda_{j}=\log(\mu_{j})/T$. The least-damped Floquet
exponent has the smallest real part and it determines the decay rate of the
transient response. If any mode is located in the unstable half-plane, _i.e._
, $\mathcal{R}e(\lambda_{j})>0$, the system is globally unstable.

It is worth noting that the transient length remains independent of the
steady-state period. This characteristic is particularly important as it
ensures that the length of integration does not rely on $T$, thereby avoiding
potential deterioration in the performance of RSVD-$\Delta t$ for flows with
low periodicity. In practice, both finding the principal fundamental matrix
and determining the Floquet exponents can be computationally expensive,
especially for large-scale systems. As with autonomous cases, it is possible
to run a homogeneous simulation to analyze the long-term behavior of the
transient response for periodic flows. This equivalence has been illustrated
on a test case in $\S$7.

If $\bm{T}$ is singular but all other modes are stable, then the Floquet
exponent with smallest real part will be zero with
$\widehat{d\bm{\bar{q}}/dt}$ as the corresponding Floquet mode. This can be
seen from (3.12), where the phase shift direction is a non-trivial solution of
the homogenous system. Leclercq & Sipp (2023) have also shown that
$d\bm{\bar{q}}/dt$ is associated with a zero Floquet exponent under the
assumption that the base flow satisfies the Navier-Stokes equations as in
(3.11).

###  6.2 Efficient transient removal

Our strategy leverages the periodic nature of the steady-state part and the
linear evolution of the transient part of the solution to directly compute and
eliminate the undesired transient portion. For a given pair of solutions
$\bm{q}_{1}$ and $\bm{q}_{2}$ at times $t_{1}$ and $t_{2}=t_{1}+T$,
respectively, they can be expressed as a sum of their steady-state and
transient components as

$\displaystyle\bm{q}_{1}$ $\displaystyle=\bm{q}_{s,1}+\bm{q}_{t,1},$ (6.2)
$\displaystyle\bm{q}_{2}$ $\displaystyle=\bm{q}_{s,2}+\bm{q}_{t,2},$

where $\bm{q}_{t,1}$ and $\bm{q}_{t,2}$ represent the transient part which
decays to zero as $t\to\infty$ and $\bm{q}_{s,1}$ and $\bm{q}_{s,2}$ denote
the steady-state part which evolves periodically. The time distance between
two snapshots is one period, hence

$\bm{q}_{s,2}=\bm{q}_{s,1}.$ (6.3)

Also, the evolution of the transient part can be expressed as

$\bm{q}_{t,2}=\bm{\Phi}(T)\bm{q}_{t,1},$ (6.4)

where $\bm{\Phi}$ is the principal fundamental matrix of (3.1). Simplifying
(6.2), (6.3), and (6.4) leads to

$(\bm{\Phi}(T)-\bm{I})\bm{q}_{t,1}=\bm{b},$ (6.5)

where $\bm{b}=\bm{q}_{2}-\bm{q}_{1}$ is known.

By obtaining $\bm{q}_{t,1}$ from solving (6.5), the steady-state solution can
be recovered as $\bm{q}_{s,1}=\bm{q}_{1}-\bm{q}_{t,1}$. The central challenge
in removing the transient is the computational cost associated with solving
the linear system in (6.5). Instead, we employ Petrov-Galerkin (or Galerkin)
projection to obtain an approximate solution to (6.5) in a more cost-effective
manner.

A low-dimensional representation of the transient part can be expressed as

$\bm{q}_{t,1}=\bm{\phi}\bm{\beta}_{1},$ (6.6)

where $\bm{\phi}\in\mathbb{C}^{N\times r}$ is an orthonormal low-dimensional
test basis with $r\ll N$, and $\bm{\beta}_{1}\in\mathbb{C}^{r}$ represents the
coefficients describing the transient response in this basis. Substituting
(6.6) into (6.5), the linear system

$(\bm{\Phi}(T)-\bm{I})\bm{\phi}\bm{\beta}_{1}=\bm{b}$ (6.7)

is overdetermined. We use Petrov-Galerkin projection with a low-dimensional
trial basis $\bm{\psi}\in\mathbb{C}^{N\times r}$ to close (6.7), yielding

$\tilde{\bm{M}}\bm{\beta}_{1}=\bm{\psi}^{*}\bm{b},$ (6.8)

where

$\tilde{\bm{M}}=\bm{\psi}^{*}(\bm{\Phi}(T)-\bm{I})\bm{\phi}$ (6.9)

is the map between coefficients. Solving (6.8) for $\bm{\beta}_{1}$ is
inexpensive due to its reduced dimension, and from (6.6) the transient
estimate is obtained as

$\bm{q}_{t,1}=\bm{\phi}(\bm{\psi}^{*}\bm{\phi}-\tilde{\bm{M}})^{-1}\bm{\psi}^{*}\bm{b}.$
(6.10)

The reduced operator $\tilde{\bm{M}}$ is obtained by integrating the columns
of $\bm{\phi}$ for one period, giving $\bm{\Phi}(T)\bm{\phi}$, and then
projecting $(\bm{\Phi}(T)-\bm{I})\bm{\phi}$ against the columns of
$\bm{\psi}$. This integration happens once and its CPU cost is equivalent to
integrating the LNS equations for $r$ periods. The analogous process occurs
for the adjoint equations.

How should the test and trial bases be selected? We have found a test basis
spanning the least-damped modes to be effective. Rather than determining these
modes through an exhaustive Floquet analysis, we employ a homogenous
simulation as illustrated in the top row of the flowchart in figure 3. By
initiating a sufficiently long simulation from a normalized random initial
condition, the system asymptotically converges to the least-damped modes. The
longer the integration length, the lower-dimensional the space becomes. As
demonstrated in a test case in $\S$7, the dimension of the test bases can, in
some cases, be as small as a single column. Our test cases employ Galerkin
projection with identical test and trial bases.

Figure 3: Flowchart depicting the action of $\bm{H}$ using time stepping (a)
with efficient transient removal strategy and (b) with natural transient
decay.

The obtained test basis effectively captures the transient response at a
specific phase $\theta$. Utilizing $\bm{\phi}$ and performing Galerkin
projection, we obtain $\tilde{\bm{M}}$. The time-stepping is then carried out
for a sufficiently long duration to allow the initial transients to vanish
before applying (6.8) to determine the transient at phase $\theta$ as shown in
the middle row of the flowchart in figure 3. Once the steady-state response at
the same phase is updated as $\bm{q}_{s,1}=\bm{q}_{1}-\bm{q}_{t,1}$, it is
used as a new initial condition, which is synchronized with the forcing and
will not excite a transient response (within the span of the bases used to
construct $\tilde{\bm{M}}$). Then, integration for one period is sufficient,
and $N_{\omega}$ steady-state snapshots are collected as desired. The
resulting steady-state snapshots are identical to those obtained by a
prolonged wait in the case of natural decay, as illustrated in the bottom row
of the flowchart in figure 3. This procedure is also applicable to the adjoint
LNS equations. Given that the adjoint equations differ from the LNS equations,
a new basis needs to be constructed for them. Implementing this strategy can
reduce the integration length by a factor of 10 or more, depending on the
desired accuracy.

##  7 Test cases

We assess the effectiveness of RSVD-$\Delta t$ in periodic systems using two
test cases. First, we use a modified Ginzburg-Landau equation to validate the
RSVD-$\Delta t$ algorithm and to illustrate the transient removal strategy.
Second, we consider a periodic flow around an airfoil, similar to a test case
from Padovan et al. (2020), to evaluate the accuracy and cost-savings of the
RSVD-$\Delta t$ algorithm compared to RSVD-LU.

###  7.1 Periodically varying complex Ginzburg-Landau equation

The one-dimensional complex Ginzburg-Landau equations is a widely used model
for understanding and controlling the non-modal growth of instabilities in
transitional and turbulent shear flows (Chomaz et al., 1988; Hunt & Crighton,
1991; Bagheri et al., 2009; Chen & Rowley, 2011; Cavalieri et al., 2019).
Here, a modified system is used as an inexpensive test case to validate our
algorithm. The original complex Ginzburg-Landau system follows the form of
(2.1) with

$\begin{gathered}\mathcal{A}=-\nu\frac{\partial}{\partial
x}+\gamma\frac{\partial^{2}}{\partial x^{2}}+\mu(x),\\\
\mu(x)=(\mu_{0}-c_{\mu}^{2})+\frac{\mu_{2}}{2}x^{2},\\\
\bm{B}=\bm{C}=\bm{I}.\end{gathered}$ (7.1)

Following Bagheri et al. (2009), we set
$c_{\mu}=0.2,\mu_{2}=-0.01,\gamma=1-\text{i},$ and $\nu=2+0.2\text{i}$. This
system is globally stable when $\mu_{0}<\mu_{0,cr}\approx 0.3977$ (Bagheri et
al., 2009). By substituting
$\mu_{0}(t)=\bar{\mu}_{0}+\mu_{p}sin(\omega_{f}t-\pi/2)$ in place of
$\mu_{0}$, we transform the Ginzburg-Landau equations to a periodically
varying system with fundamental frequency $\omega_{f}$, following the form of
(3.1). We set $\mu_{0}=0.395$, $\mu_{p}=0.1$, and $\omega_{f}=0.1$; the mean
value of $\mu(t)$ is equal to $\mu_{0}$, whose value is close to the critical
value $\mu_{0,cr}$ to promote significant growth. The periodic linear
operators $\bm{A}_{p,j}$ are constructed using a central finite-difference
method for $x\in[-50,50]$ with $N=1000$ grid points, and we set
$\bm{W}=\bm{I}$ on account of the uniform grid.

Nine operators are built within the interval $[0,T)$, with
$T=2\pi/\omega_{0}\approx 62.83$, such that $\bm{A}_{p}(t)=\bm{A}_{p}(t+T)$,
before conducting a DFT to obtain $\hat{\bm{A}}_{p}$. Figure 4(a) depicts the
Frobenius norm of $\hat{\bm{A}}_{p}$ up to the fourth harmonic, normalized
relative to the maximum norm. The spectrum reveals that only the mean and the
first harmonic are active, as expected since $\mu_{0}(t)$ oscillate at a
single frequency. Accordingly, $\bm{A}_{p}$ can be reconstructed at any given
time using $\hat{\bm{A}}_{p,0}$, and $\hat{\bm{A}}_{p,\pm 1}$, and the
harmonic resolvent operator is constructed around
$\Omega_{\bar{q}}=\\{-1,0,1\\}\omega_{f},$ _i.e._ , $N_{b}=3$. Notably, the
first harmonic’s norm is several orders of magnitude smaller than that of the
base flow. Nonetheless, it should not be disregarded, as we will elucidate its
significance in subsequent discussions.

Figure 4: Ginzburg-Landau test case: (a) the normalized Frobenius norm of
$\hat{\bm{A}}_{p,\omega}$ up to the fourth harmonic. (b) The norm of the
transient response over time (blue) along with the expected decay rate from
Floquet analysis (red).

An essential initial step involves ensuring the stability of $\bm{A}_{p}$.
Typically, this is achieved by integrating the homogeneous system, _i.e._ ,
integrating (3.1) with zero forcing, starting from a normalized random initial
condition. However, the compact size of this system allows us to validate the
time-stepping approach for computing the least-damped decaying mode against a
Floquet stability analysis. The least-damped Floquet exponent, computed as
$\lambda_{l}=-0.0021-0.0477\text{i}$, is expected to govern the decay rate of
the transient response. Figure 4(b) displays the norm of snapshots of the
homogeneous system response over time, along with a reference line
representing $e^{\lambda_{l}t}$. The natural decay rate aligns with the
reference line, as expected, ensuring that the transient run will ultimately
converge to the least-damped Floquet mode. Additionally, the system is not
singular, as the transient dynamics never reach a naturally stable state, and
none of the Floquet exponents have a zero real value.

Before proceeding to the computation of harmonic resolvent modes and gains
using RSVD-$\Delta t$ and RSVD-LU, the transient removal strategy is
demonstrated, and computing the actions of $\bm{H}$ and $\bm{H}^{*}$ are
validated for a given random forcing containing
$\Omega_{q}\in\\{-10,-9,...,9,10\\}\omega_{f}$, _i.e._ , $\omega\in[-1,1]$
with $\Delta\omega=0.1$. The responses have been computed up to the $20^{th}$
harmonic, _i.e._ , $\omega\in[-2,2]$, using both time stepping and directly in
Fourier space as outlined in (3.7). We used a $4^{th}$ order Runge–Kutta (RK4)
integration scheme and a time step of $dt=0.001$ and $T_{t}=15000$ for this
purpose.

Figure 5(a) depicts the spectrum of response norms. In order to highlight the
substantial difference in output generated by the harmonic resolvent operator
around the base flow ($\Omega_{\bar{q}}=\\{0\\}\omega_{f},N_{b}=1$), the norm
of the responses for the same input forcing is also presented. Figure 5(b)
displays the relative errors. The harmonic resolvent with $N_{b}=1$ produces
inaccurate results across the entire spectrum and has zero response for
$|\omega|>1$ since the modes at different frequencies are entirely decoupled.
Therefore, neglecting $\hat{\bm{A}}_{p,\pm 1}$ is impractical, as mentioned
earlier, and induces significant alterations in the harmonic resolvent modes
and gains. On the other hand, the time-stepping output exhibits almost
machine-precision accuracy within $|\omega|\leq 1$, and its accuracy
diminishes as response norms decrease within $|\omega|>1$. The validation of
$\bm{H}^{*}$ follows a similar procedure but is omitted here for brevity.

Figure 5: Ginzburg-Landau test case: (a) the norms of the Fourier space
response for $N_{b}=3$ (black) and $N_{b}=1$ (blue), along with the norm of
the steady-state response obtained via time-stepping (red). (b) The relative
errors of the norms, with the norms from Fourier space with $N_{b}=3$ (black
line in (a)) serving as the reference.

The transient norm decay rate, as shown in figure 4(b), follows
$e^{-0.0021t}$. This indicates that approximately 40 periods are necessary for
the transient norm to decay to 1%. However, we can significantly reduce this
duration by employing our transient removal strategy. To begin, we precompute
the test basis by conducting a sufficiently long homogeneous simulation, and
we find that a single mode suffices. Next, we integrate the orthogonalized
basis for one period to obtain $\tilde{M}$ and complete the preprocessing
step.

Employing the same forcing as the validation case, we obtain the responses in
time at the identical phase (initiated at $t=0$ or $\theta=0$ with a time
interval of $\Delta t=T$), and compare them with the solutions after the
removal of transients. Figure 6(a) presents the norm of the responses, which
illustrates the effectiveness of our strategy commencing from the second
period. This plot suggests that after two periods using our transient removal
strategy with Galerkin projection, we obtain an accurate steady-state snapshot
which can be used as the synchronized initial condition for the forced
equations. Figure 6(b) demonstrates that the snapshots ultimately align
perfectly with the updated snapshot. Hence, using our strategy, a total of
three periods are required – two periods before transient removal and one
period after the initial condition is obtained – to acquire all $N_{\omega}$
steady-state snapshots, a 10-fold speed-up compared to the natural decay. A
similar observation was made regarding the adjoint equations. In brief, our
strategy can significantly reduce computational costs, potentially by one or
more orders of magnitude, depending on the desired level of accuracy.

Figure 6: Ginzburg-Landau test case: (a) the response norms before (black) and
after (red) transient removal at the same phase, _i.e._ , snapshots with
$\Delta t=T$ interval. (b) Comparison of the normalized response at the second
period in blue, the $50^{th}$ period in black, and the second period after
transient removal in red.

Finally, RSVD-$\Delta t$ employing an RK4 integration scheme with $dt=0.003$
and $T_{t}=2T$, is carried out along with our efficient transient removal
strategy to compute the harmonic resolvent modes for $k=5$ and $q=1$. To
establish a reference, the harmonic resolvent modes are also computed using
RSVD-LU with an equivalent number of test vectors and power iterations. Figure
7(a) illustrates that RSVD-$\Delta t$ computes the five leading gains
identically to RSVD-LU. The gain relative error in figure 7(b) confirms that
relative errors remain below $10^{-9}$ across all modes. The inner products of
response and forcing modes, computed via RSVD-$\Delta t$ and RSVD-LU,
demonstrate parallel directions as the relative errors remain below $10^{-8}$
across all modes. The inner-product error between two unit-norm vectors
$\bm{v}_{1}$ and $\bm{v}_{2}$ is defined as

$e_{ip}=1-\langle\bm{\bm{v}_{1}},\bm{v_{2}}\rangle.$ (7.2)

Figure 7: Ginzburg-Landau test case: (a) the five leading gains using
RSVD-$\Delta t$ (blue) and RSVD-LU (red). (b) The corresponding gain relative
error (purple) and relative inner product errors (7.2) between the response
modes (cyan) and the forcing modes (green) computed via RSVD-$\Delta t$ and
RSVD-LU.

###  7.2 Flow over an airfoil

Our second test case consists of the flow over a NACA0012 airfoil at a low
Reynolds number of $Re=200$ and a steep angle of attack of
$\alpha=20^{\circ}$, which is dominated by periodic motions. This problem
serves as a benchmark for computational cost and accuracy comparisons. A
direct numerical simulation is conducted using the “CharLES” compressible flow
solver. To mimic the incompressible simulations by Padovan et al. (2020), we
set the Mach number to 0.05. The simulation employs a C-shaped mesh with a
total of 62,000 cells. The leading edge of the airfoil is located at the
origin $(x/L_{c},y/L_{c})=(0,0)$, where $L_{c}=1$ is the chord length. The
computational domain spans $x/L_{c}\times y/L_{c}\in[-49,50]\times[-50,50]$. A
constant time step of $\Delta tU_{\infty}/L_{c}=6.88\times 10^{-5}$ is used,
corresponding to a Courant-Friedrichs-Lewy number of 0.91, where $U_{\infty}$
is the inflow streamwise velocity. Integration continues for a sufficient
duration until the flow reaches a periodic state. Figure 8(a) displays the
power spectral density computed from the transverse velocity at
$(x/L_{c},y/L_{c})=(3.0,-0.43)$, where a pronounced vortex shedding occurs
behind the trailing edge. The shedding frequency is $St_{f}=0.114$, where the
Strouhal number is defined as $St=\frac{\omega L_{c}sin(\alpha)}{2\pi
U_{\infty}}$. The consistency between the dominant frequency identified in the
power spectral density and obtained from our linearized code has been
confirmed in our previous work (Jung et al., 2023).

Figure 8: Airfoil test case: (a) the PSD spectrum based on transverse velocity
at $(x/L_{c},y/L_{c})=(3.0,-0.43)$. (b) The normalized Frobenius norm of
$\hat{\bm{A}}_{p,St}$ up to the eighth harmonic.

The linearized operators are constructed at 100 time points within
$T=2\pi/St_{f}\approx 55$. Performing a DFT on the $\bm{A}_{i}$, we generate
100 LNS operators $\hat{\bm{A}}_{j}$ in Fourier space to be used for harmonic
resolvent analysis. The Frobenius norms $||\hat{\bm{A}}_{p,St}||_{F}$,
depicted in figure 8, indicate that $\bm{A}_{p}(t)$ can be accurately
reconstructed using up to the 5th harmonic. Given that the linearized operator
$\bm{A}_{p}(t)$ is real-valued for this problem, we retained only $\lfloor
N_{b}/2\rfloor+1=6$ coefficients, _i.e._ , the zeroth and 5 positive
harmonics. For time stepping, $dt=0.0045$ is chosen to ensure the stability
and accuracy of the RK4 integration scheme. The input and output perturbation
frequency range spans up to the $7^{th}$ harmonic, _i.e._ , $N_{\omega}=15$.
The domain of interest is $x/L_{c}\times y/L_{c}\in[-4,12]\times[-2.5,2.5]$,
identical to that in Padovan et al. (2020). We seek optimal forcing and
response under the Chu energy norm (Chu, 1965), which has been used in several
previous studies (Towne et al., 2018; Schmidt et al., 2018; Heidt & Colonius,
2023).

A homogeneous simulation of the time-periodic system exhibits a slow decay
rate, so we employ the transient removal strategy. Three snapshots constitute
the test basis, obtained after a sufficiently long time interval where most
initial transients are eliminated. A random forcing, encompassing the set
$\Omega_{q}\in\\{-7,-6,...,6,7\\}St_{f}$ frequencies, is applied to the LNS
equations for over 30 periods. The snapshot norms slowly approach towards the
steady-state norm but remain far from converging, as illustrated in figure 9.
Utilizing our strategy, the relative norm error falls below the 1% threshold
after 20-30 periods in most cases with different random forcings. The efficacy
of our strategy is further emphasized in figure 9(b), where the 1% relative
error is achieved before 10 periods for the optimal forcing (from RSVD-$\Delta
t$ output). Ultimately, we set $T_{t}=29T$ and conducted our simulations for a
total duration of 30 periods to compute the actions of $\bm{H}$ and
$\bm{H}^{*}$. Extrapolating the natural decay rate suggests that, in the
absence of our transient-removal strategy, the norms reach the 1% threshold
after around 2000 periods, demonstrating over a 60-fold acceleration in time
stepping and overall algorithm performance with our strategy. This observation
holds true for the adjoint equations as well.

Figure 9: Airfoil test case: (a) the response norms before (black) and after
(red) transient removal, and the true steady-state norm (blue), at the same
phase, _i.e._ , snapshots with $\Delta t=T$ interval, for a random forcing.
(b) The relative error between steady-state response norm and the response
norms before (black) and after (red) transient removal.

The number of test vectors is set to $k=5$, and $q=2$ power iteration proves
sufficient for the convergence of both gains and modes. The suboptimal forcing
and response modes are shown in figure 10. The results from RSVD-$\Delta t$
closely align with the existing data obtained from RSVD-LU by Padovan et al.
(2020), despite variations in the CFD solver, boundary conditions, domain
setup, and energy norm. The gain plot in figure 11(a) exhibits a similar
pattern to Padovan et al. (2020), with the exception that we have retained the
optimal gain associated with the phase shift (the exact values differ due to
differences in the problem setup). We noted that the leading mode and gain,
_i.e._ , the mode associated with the phase shift, converged without power
iterations due to a substantial two-orders-of-magnitude gap between the
optimal and suboptimal gains.

Figure 10: Airfoil test case: real part of the vorticity field computed from
(a, c, e, g) the input mode and (b, d, f, h) the output mode associated with
the first suboptimal gain.

From a computational standpoint, a fair comparison between the RSVD-LU and
RSVD-$\Delta t$ algorithms is feasible when both utilize the same parameters.
However, the use of $N_{b}=11$ proves to be excessively memory-intensive for
the RSVD-LU algorithm, exceeding the available 3.5 TB of memory on our
cluster. The CPU and memory scaling plots in figure 11(b) and (c),
respectively, depict the LU decomposition of $\bm{T}$ with $N_{b}=3$,
maintaining a constant $N$ while varying the number of blocks $N_{\omega}$
from 5 to 15 with $\Delta N_{\omega}=2$. On the same plots, we present the
cost of computing the action of $\bm{H}$ on a vector using time-stepping. The
total duration is set to 30 periods to obtain an accurate solution. We only
vary $N_{\omega}$ while keeping $N_{b}=3$, the overall cost of RSVD-$\Delta t$
remains constant and unaffected by changes in $N_{\omega}$. This implies that
the creation of LNS operators and the time-stepping process incur
significantly higher costs compared to the transformations between Fourier
space and the time domain for forcing and response. The memory of RSVD-$\Delta
t$ to store LNS operators remains independent of $N_{\omega}$, as expected
(see $\S$5.1), while storing $\bm{Q}$ and $\bm{F}$ matrices grows linearly
with number of frequencies $N_{\omega}$.

RSVD-LU is estimated to require 1543 CPU-hours for the LU decomposition of
$\bm{T}$, and an additional 15 CPU-hours for solving the LU-decomposed system
for each test vector, yielding a total cost of $\text{CPU}_{\text{{RSVD-
LU}}}\approx 1543+15\times 5\times 2\times 2=1843$ CPU-hours for $k=5,q=1$. On
the other hand, employing RSVD-$\Delta t$ with the aforementioned time-
stepping parameters and $N_{b}=11$ and $N_{\omega}=15$ incurs a cost of
approximately 7 CPU-hours per period for a test vector, resulting in a total
cost of $\text{CPU}_{\text{{RSVD-$\Delta t$}}}\approx 30\times 7\times 5\times
2\times 2=4200$ CPU-hours. The cost of QR decomposition and the final SVD are
ignored as they are orders of magnitude faster. In terms of memory
consumption, RSVD-LU peaks at $\text{RAM}_{\text{RSVD-LU}}\approx 2036GB$ for
the LU decomposition of $\bm{T}$. RSVD-$\Delta t$ , on the other hand, stores
$\lfloor N_{b}/2\rfloor+1=6$ sparse matrices, each of size 0.4 GB, and three
dense matrices of size $Nk(\lfloor N_{\omega}/2\rfloor+1)=313630\times 5\times
8$, totaling $\text{RAM}_{\text{{RSVD-$\Delta t$}}}\approx 0.4\times
6+0.185\times 3\approx 3GB$. This translates to a memory saving of
approximately three orders of magnitude, even with an unbalanced number of
base frequencies. Overcoming the memory consumption hurdle, which is typically
the limiting factor in practice, RSVD-$\Delta t$ emerges as a viable tool for
analyzing larger-scale flows compared to RSVD-LU.

Figure 11: Airfoil test case: (a) the five leading gains using RSVD-$\Delta t$
. (b) The CPU-hour, and (c) memory usage scaling of RSVD-LU (red) and
RSVD-$\Delta t$ (blue) to compute the action of $\bm{H}$ onto a vector, _i.e._
, $k=1$. The memory usage of RSVD-$\Delta t$ is decomposed into memory
required to store LNS operators (solid) and forcing and response matrices
(dashed).

##  8 Conclusions

This paper presents an extension of the RSVD-$\Delta t$ algorithm that was
originally developed for resolvent analysis of steady base flows to handle
harmonic resolvent analysis of periodic flows. Specifically, we demonstrate
that the time stepping technique employed within RSVD-$\Delta t$ can
effectively replace the actions of $\bm{H}$ and its complex conjugate
$\bm{H}^{*}$ in Fourier space. In terms of CPU cost, the RSVD-$\Delta t$
algorithm exhibits a scaling of $O(N)$ for both resolvent and harmonic
resolvent analyses. This offers a significant advantage, considering that the
LU decomposition of $\bm{T}$ scales as $O((N_{\omega}N)^{c})$ with $c\geq
1.5$. Regarding memory usage, our algorithm only stores relevant matrices in
Fourier space and exhibits $O(NN_{\omega})$ scaling; in contrast, the memory
peak of LU decomposition of $\bm{T}$ empirically scales with
$O((N_{\omega}N)^{1.5})$ or worse. One difference between applying
RSVD-$\Delta t$ to resolvent and harmonic resolvent analyses lies in the
necessity to update the LNS operators during the time-stepping process for the
latter. These operators are efficiently generated on the fly from their
$N_{b}$ Fourier components.

The error sources in our algorithm align with those in resolvent analysis,
extensively investigated in Farghadan et al. (2023). A unique contribution of
this paper is the introduction of a novel transient removal strategy for
harmonic resolvent analysis. While reminiscent of our approach for resolvent
operators, the challenge in harmonic resolvent analysis lies in the
intertwining of all retained frequencies. Our strategy takes advantage of the
differing evolution of the steady-state and transient components of the
response and Petrov-Galrkin or Galerkin projections.

A validation of RSVD-$\Delta t$ against RSVD-LU is conducted using a periodic
Ginzburg-Landau system. Additionally, we verify the governing role of the
Floquet exponent in determining the decay rate of the least-damped mode of the
system. Extending the application of our algorithm, we analyze a two-
dimensional flow passing an airfoil, providing insights into the CPU and
memory complexities. Despite dealing with a mid-sized flow scenario, a
substantial memory gap persists between RSVD-LU and RSVD-$\Delta t$
algorithms. The computed forcing and response modes corresponding to the first
suboptimal gain closely resemble those obtained by Padovan et al. (2020). The
adaptation of the transient removal strategy for these test cases
significantly enhanced the performance of RSVD-$\Delta t$, resulting in
10-fold speed-up for the Ginzburg-Landau problem and 60-fold speed-up for the
airfoil to reach a 1% relative error. The speed-up is even greater for lower
error tolerances.

Moreover, our algorithm exhibits versatility, accommodating non-identity
weight, input, and output matrices. This capability extends to computing modes
for the modified harmonic resolvent operator, and in particular, the analysis
of cross-frequency amplification mechanisms via
$\bm{H}_{\omega_{2},\omega_{1}}$. Our algorithm tackles the computation of
subharmonic resolvent modes if desired. Implementation of the RSVD-$\Delta t$
algorithm within Petsc (Balay et al., 2019) and Slepc (Hernandez et al., 2005)
environments leverage parallel computations for enhanced efficiency. Our time-
stepping approach allows for a matrix-free application. This can be
implemented using any code equipped with linear direct and adjoint
capabilities, bypassing the explicit formation of $\bm{T}$ (de Pando et al.,
2012; Martini et al., 2021). Finally, the time-stepping within our algorithm
can also be used along with other SVD algorithms, _e.g._ , Arnoldi and power
iteration, if desired.

In turbulent flows, where the state dimension $N$ can become very large due to
higher resolution requirements, CPU cost and memory requirements can swiftly
impose constraints on the applicability of RSVD-LU. This limitation applies
even to steady problems and is further exacerbated for harmonic resolvent
analysis due to the inflated frequency-domain operators necessitated by
triadic frequency coupling. Our algorithm circumvents this issue by
strategically operating in the time domain, avoiding the need for limiting
operations like LU decomposition, leading to linear cost scaling. Because of
this key distinction, the RSVD-$\Delta t$ could enable harmonic resolvent
analysis of previous intractable turbulent flows.

## Acknowledgements

This work was funded in part by Air Force Office of Scientific Research
(AFOSR) grant no. FA9550-20-1-0214 and by a Catalyst Grant from the Michigan
Institute for Computational Discovery and Engineering (MICDE). We acknowledge
the University of Michigan’s Great Lakes cluster for providing the essential
computational resources to conduct all the simulations in this study.

##  Appendix A RSVD-$\Delta t$ for the subharmonic resolvent operator

In $\S$3.5, we provided an overview of subharmonic resolvent analysis. In this
appendix, we briefly outline the application of RSVD-$\Delta t$ in computing
subharmonic resolvent modes and gains. Consider the frequency of interest as
$\gamma\in\Omega_{\gamma}$. This set specifies the perturbation frequency,
while the base flow frequency content is confined to $\Omega_{\bar{q}}$.

To compute the actions of $\bm{H}$ and $\bm{H}^{*}$ using time-stepping, we
must adhere to both the base flow frequency, enforcing a duration of
$T=2\pi/\omega_{f}$, and the perturbation frequency, enforcing a duration of
$T_{p}=2\pi/\gamma$, in order to obtain $N_{\omega}$ steady-state solutions to
(3.1). Thus, we need to integrate for $T_{sub}$ such that
$T_{sub}/T\in\mathbb{N}$ and $T_{sub}/T_{p}\in\mathbb{N}$, during which
$N_{\omega}$ equidistant snapshots are saved. For instance, if we consider
$\gamma=\omega_{f}/5$, requiring $T_{p}=2\pi/(\omega_{f}/5)=5T$, integrating
for $T_{sub}=5T$ is sufficient to obtain the steady-state solutions. All the
other steps remain the same as introduced in Algorithm 1.

## References

  * Bagheri et al. (2009) Bagheri, S., Henningson, D. S., Hoepffner, J. & Schmid, P. J. 2009 Input-output analysis and control design applied to a linear model of spatially developing flows. Applied Mechanics Reviews 62 (2).
  * Balay et al. (2019) Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W. & others 2019 PETSc users manual. Argonne National Laboratory.
  * Cavalieri et al. (2019) Cavalieri, A. V. G., Jordan, P. & Lesshafft, L. 2019 Wave-packet models for jet dynamics and sound radiation. Applied Mechanics Reviews 71 (2).
  * Chen & Rowley (2011) Chen, K. K. & Rowley, C. W. 2011 H2 optimal actuator and sensor placement in the linearised complex Ginzburg–Landau system. Journal of Fluid Mechanics 681, 241–260.
  * Chomaz et al. (1988) Chomaz, J. M., Huerre, P. & Redekopp, L. G. 1988 Bifurcations to local and global modes in spatially developing flows. Physical Review Letters 60 (1), 25.
  * Chu (1965) Chu, B.-T. 1965 On the energy transfer to small disturbances in fluid flow (part i). Acta Mechanica 1 (3), 215–234.
  * Dawson & McKeon (2019) Dawson, S. T. M. & McKeon, B. J. 2019 On the shape of resolvent modes in wall-bounded turbulence. Journal of Fluid Mechanics 877, 682–716.
  * Farghadan & Arzani (2019) Farghadan, A. & Arzani, A. 2019 The combined effect of wall shear stress topology and magnitude on cardiovascular mass transport. International Journal of Heat and Mass Transfer 131, 252–260.
  * Farghadan et al. (2023) Farghadan, A., Martini, E. & Towne, A. 2023 Scalable resolvent analysis for three-dimensional flows. arXiv preprint arXiv:2309.04617 .
  * Franceschini et al. (2022) Franceschini, L., Sipp, D., Marquet, O., Moulin, J. & Dandois, J. 2022 Identification and reconstruction of high-frequency fluctuations evolving on a low-frequency periodic limit cycle: application to turbulent cylinder flow. Journal of Fluid Mechanics 942, A28.
  * Giannenas et al. (2022) Giannenas, A. E., Laizet, S. & Rigas, G. 2022 Harmonic forcing of a laminar bluff body wake with rear pitching flaps. Journal of Fluid Mechanics 945, A5.
  * Halko et al. (2011) Halko, N., Martinsson, P. & Tropp, J. A. 2011 Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review 53 (2), 217–288.
  * Heidt & Colonius (2023) Heidt, L. & Colonius, T. 2023 Spectral proper orthogonal decomposition of harmonically forced turbulent flows. arXiv preprint arXiv:2305.05628 .
  * Hernandez et al. (2005) Hernandez, V., Roman, J. E. & Vidal, V. 2005 Slepc: A scalable and flexible toolkit for the solution of eigenvalue problems. ACM Transactions on Mathematical Software (TOMS) 31 (3), 351–362.
  * Hunt & Crighton (1991) Hunt, R. E. & Crighton, D. G. 1991 Instability of flows in spatially developing media. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 435 (1893), 109–128.
  * Jung et al. (2023) Jung, J., Bhagwat, R. & Towne, A. 2023 Resolvent-based estimation of laminar flow around an airfoil. AIAA Paper $\\#$2023-0077 .
  * Leclercq & Sipp (2023) Leclercq, C. & Sipp, D. 2023 Mean resolvent operator of a statistically steady flow. Journal of Fluid Mechanics 968, A13.
  * Lin et al. (2023) Lin, C.-T., Tsai, M.-L. & Tsai, H.-C. 2023 Flow control of a plunging cylinder based on resolvent analysis. Journal of Fluid Mechanics 967, A41.
  * Martini et al. (2021) Martini, E., Rodríguez, D., Towne, A. & Cavalieri, A. V. G. 2021 Efficient computation of global resolvent modes. Journal of Fluid Mechanics 919.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. Journal of Fluid Mechanics 658, 336–382.
  * Moarref et al. (2013) Moarref, R., Sharma, A. S., Tropp, J. A. & McKeon, B. J. 2013 Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels. Journal of Fluid Mechanics 734, 275–316.
  * Monokrousos et al. (2010) Monokrousos, A., Åkervik, E., Brandt, L. & Henningson, D. S. 2010 Global three-dimensional optimal disturbances in the blasius boundary-layer flow using time-steppers. Journal of Fluid Mechanics 650, 181–214.
  * Nyquist (1928) Nyquist, H. 1928 Certain topics in telegraph transmission theory. Transactions of the American Institute of Electrical Engineers 47 (2), 617–644.
  * Padovan et al. (2020) Padovan, A., Otto, S. E. & Rowley, C. W. 2020 Analysis of amplification mechanisms and cross-frequency interactions in nonlinear flows via the harmonic resolvent. Journal of Fluid Mechanics 900.
  * Padovan & Rowley (2022) Padovan, A. & Rowley, C. W. 2022 Analysis of the dynamics of subharmonic flow structures via the harmonic resolvent: Application to vortex pairing in an axisymmetric jet. Physical Review Fluids 7 (7), 073903.
  * de Pando et al. (2012) de Pando, M. F., Sipp, D. & Schmid, P. J. 2012 Efficient evaluation of the direct and adjoint linearized dynamics from compressible flow solvers. Journal of Computational Physics 231 (23), 7739–7755.
  * Pickering et al. (2021) Pickering, E., Towne, A., Jordan, P. & Colonius, T. 2021 Resolvent-based modeling of turbulent jet noise. The Journal of the Acoustical Society of America 150 (4), 2421–2433.
  * Ribeiro et al. (2020) Ribeiro, J. H. M., Yeh, C.-A. & Taira, K. 2020 Randomized resolvent analysis. Physical Review Fluids 5 (3), 033902.
  * Schmidt et al. (2018) Schmidt, O. T., Towne, A., Rigas, G., Colonius, T. & Brès, G. A. 2018 Spectral analysis of jet turbulence. Journal of Fluid Mechanics 855, 953–982.
  * Sipp & Marquet (2013) Sipp, D. & Marquet, O. 2013 Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer. Theoretical and Computational Fluid Dynamics 27 (5), 617–635.
  * Taira et al. (2017) Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W., Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S., Theofilis, V. & Ukeiley, L. S. 2017 Modal analysis of fluid flows: An overview. AIAA Journal 55 (12), 4013–4041.
  * Towne et al. (2022) Towne, A., Rigas, G., Kamal, O., Pickering, E. & Colonius, T. 2022 Efficient global resolvent analysis via the one-way Navier-Stokes equations. Journal of Fluid Mechanics 948, A9.
  * Towne et al. (2018) Towne, A., Schmidt, O. T. & Colonius, T. 2018 Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. Journal of Fluid Mechanics 847, 821–867.
  * Troy & Koseff (2005) Troy, C. D. & Koseff, J. R. 2005 The instability and breaking of long internal waves. Journal of Fluid Mechanics 543, 107–136.
  * Wereley (1990) Wereley, N. M. 1990 Analysis and control of linear periodically time varying systems. PhD thesis, Massachusetts Institute of Technology.