{'fraction_non_alphanumeric': 0.03253345639132441, 'fraction_numerical': 0.010529848554767797, 'mean_word_length': 4.8575992136626125, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling mode interactions in boundary layer flows via the Parabolized Floquet Equations

In this paper, we develop a model based on successive linearization to study interactions between different modes in boundary layer flows. Our method consists of two steps. First, we augment the Blasius boundary layer profile with a disturbance field resulting from the linear Parabolized Stability Equations (PSE) to obtain the modified base flow; and, second, we draw on Floquet decomposition to capture the effect of mode interactions on the spatial evolution of flow fluctuations via a sequence of linear progressions. The resulting Parabolized Floquet Equations (PFE) can be conveniently advanced downstream to examine the interaction between different modes in slowly varying shear flows. We apply our framework to two canonical settings of transition in boundary layers; the H-type transition scenario that is initiated by exponential instabilities, and streamwise elongated laminar streaks that are triggered by the lift-up mechanism. We demonstrate that the PFE capture the growth of various harmonics and provide excellent agreement with the results obtained in direct numerical simulations and in experiments.

In this paper, we develop a linear model to study interactions between different modes in slowlygrowing boundary layer flows. Our method consists of two steps. First, we augment the Blasius boundary layer profile with a disturbance field resulting from the linear Parabolized Stability Equations (PSE) to obtain the modified base flow; and, second, we draw on Floquet decomposition to capture the effect of mode interactions on the spatial evolution of flow fluctuations via a linear progression. The resulting Parabolized Floquet Equations (PFE) can be conveniently advanced downstream to examine the interaction between different modes in slowly-growing spatially evolving flows. We apply our framework to two canonical settings of transition in boundary layers; the H-type transition scenario that is initiated by exponential instabilities, and streamwise elongated laminar streaks that are triggered by the lift-up mechanism. We demonstrate that the PFE capture the growth of various harmonics and provide excellent agreement with the results obtained in direct numerical simulations and in experiments.

I. INTRODUCTION
A thorough understanding of the mechanisms driving laminar-turbulent transition in boundary layer flows is crucial for the prediction of the point of transition and for the design of air and water vehicles. In the past thirty years, remarkable progress has been made on simulating the physics of transitional flows using models with various levels of fidelity. In spite of this, the multi-layer nature of transition and the inherent complexity of the Navier-Stokes (NS) equations have hindered the development of practical control strategies for delaying transition in boundary layer flows [1]. Direct Numerical Simulations (DNS) have opened the way to accurate investigations of the underlying physics of transitional flows [2][3][4]. However, due to their high complexity and large number of degrees of freedom, nonlinear dynamical models that are based on the NS equations are not suitable for analysis, optimization, and control. On the other hand, nontrivial challenges, including lack of robustness, may arise in the model-based control of reduced-order models that are obtained using data-driven techniques [5].
Linearization of the NS equations around the mean-velocity profile results in models that are well-suited for analysis and control synthesis using tools from modern robust control [1]. In particular, stochastically forced linearized NS equations have been used to capture structural and statistical features of transitional [6][7][8][9] and turbulent [10][11][12][13] channel flows. In these models, stochastic forcing may be utilized to model the impact of exogenous excitation sources or to capture the effect of nonlinear terms in the NS equations. Moreover, in conjunction with the parallelflow assumption, the linearized NS equations are convenient for modal and non-modal stability analysis of spatially evolving flows [14,15]. However, this approach does not account for the effect of the spatially evolving base flow on the stability of the boundary layer. Global stability analysis addresses this issue by accounting for the spatially varying nature of the base flow and discretizing all inhomogeneous spatial directions [16][17][18]. Although accurate, this approach leads to problem sizes that may be prohibitively large for flow control and optimization.
In the boundary layer flow, primary disturbances are instigated via receptivity processes that involve internal or external perturbations [19], e.g., acoustic noise, free-stream turbulence, and surface roughness. Depending on the amplitude of these disturbances, the transition process may take various paths to breakdown [20,21]. In particular, primary disturbances can be amplified through modal instability mechanisms or they may experience non-modal amplification, e.g., via transient growth, the lift-up [22,23] and Orr mechanisms [24,25]. Both pathways can intensify disturbances beyond the critical threshold, trigger secondary instabilities, and induce a strong energy transfer from the mean flow into secondary modes [26]. The H-type [26][27][28] and K-type [28,29] transition scenarios are typical cases that are triggered by secondary instability mechanisms. Such mechanisms have also been shown to play an important role in the breakdown of laminar streaks at the later stages of transition [30][31][32][33][34]. All of these are initiated after the significant growth of the primary disturbances which intensify the role of nonlinear interactions. The modulation of the base flow by the primary perturbations precludes the usual normal-mode assumptions made in the derivation of the Orr-Sommerfeld equation. Instead, the physics of these secondary growth mechanisms have been commonly studied using Floquet analysis [26,30,35] and the parabolized stability equations (PSE) [36][37][38].
The PSE were introduced to account for non-parallel and nonlinear effects and thereby overcome challenges associated with analyses based on a parallel-flow assumption. In particular, the PSE were developed as a means to refine predictions of parallel flow analysis in slowly varying flows [37][38][39], e.g., in the laminar boundary layer. The PSE have also been adapted to account for the dynamics of three-dimensional flows that depend strongly on two spatial directions [40][41][42], and more recently, they have been used to model the amplification of disturbances in DNS and wall-modeled large-eddy simulation of transitional boundary layers [43]. In spite of these successes, the nonlinear nature of this framework has hindered their utility in systematic optimal flow control design. In general, the linear PSE provide reasonable predictions for the evolution of individual primary modes such as Tollmien-Schlichting (TS) waves [37]. Moreover, the predictive capability of the linear PSE has been further refined by modeling the effect of nonlinear terms as a stochastic source of excitation [44]. However, secondary growth mechanisms that lead to laminar-turbulent transition of the boundary layer flow originate from interactions between different modes and these interactions cannot be explicitly accounted for using such techniques.
In the transitional boundary layer, primary instability mechanisms can cause disturbances to grow to finite amplitudes and saturate at steady or quasi-steady states. Floquet stability analysis identifies secondary instability modes as the eigen-modes of the NS equations linearized around the modified base flow profile that contains spatially periodic primary velocity fluctuations. In the corresponding eigenvalue problem, the operators inherit a periodic structure from the underlying periodicity of the base flow and, as a result, capture primary-secondary mode interactions. Such representations that account for mode interactions also appear in the analysis of distributed systems with spatially or temporally periodic coefficients [45,46] as well as in the model-based design of sensor-free periodic strategies for controlling transitional and turbulent wall-bounded shear flows [12,[47][48][49].
In this paper, we take a step toward developing low-complexity models that account for harmonic interactions via a linear progression. To capture the dominant mode interactions while taking into account non-parallel effects, we introduce a computational framework that combines concepts from Floquet analysis and the linear PSE. The resulting equations are advanced downstream via a marching procedure. Our framework thus inherits the ability to account for mode interactions from Floquet theory while maintaining the low-complexity of the linear PSE. As a result, the proposed approach not only captures the essential physics of transitional boundary layer flows, but also opens the door to model-based control design.
Our presentation is organized as follows. In Sec. II, we describe the linearized NS equations and the linear PSE. In Sec. III, we derive the proposed Parabolized Floquet Equations (PFE) and explain the key features of our modeling framework. In Sec. IV, we examine the growth of subharmonic modes in a typical H-type transition scenario and, in Sec. V, we employ our framework to study the formation of streaks in the boundary layer flow. We conclude the paper with remarks and outline of future research directions in Sec. VI.

II. BACKGROUND
We first present the equations that govern the dynamics of flow fluctuations in incompressible flows of Newtonian fluids and then provide details on our proposed model for the downstream marching of spatially growing fluctuations in the boundary layer flow.
In a flat-plate boundary layer, with geometry shown in Fig. 1, the dynamics of infinitesimal fluctuations around a two-dimensional base flowū = where v = [ u v w ] T is the vector of velocity fluctuations, p denotes pressure fluctuations, u, v, and w are the streamwise (x), wall-normal (y), and spanwise (z) components of the fluctuating velocity field, and Re 0 is the Reynolds number at the inflow location x 0 . The Reynolds number is defined as Re = U ∞ δ/ν, where δ = νx/U ∞ is the Blasius length scale at the streamwise location x, U ∞ is the free-stream velocity, and ν is the kinematic viscosity. Spatial coordinates are non-dimensionalized by the Blasius length scale δ 0 at the inflow location x 0 , velocities by U ∞ , time by δ 0 /U ∞ , and pressure by ρU 2 ∞ , where ρ is the fluid density. It is customary to use the parallel-flow approximation to study the local stability of boundary layer flows to small amplitude disturbances [14]. This approximation, in conjunction with Floquet theory, has also been used to investigate secondary instabilities that inflict transition [14,26]. However, the parallel-flow approximation excludes the effect of the evolution of the base flow on the amplification of disturbances. This issue can be addressed via global stability analysis which accounts for the spatially varying nature of the base flow by discretizing all inhomogeneous directions. Nevertheless, global analysis of spatially-evolving flows may be prohibitively expensive for analysis, optimization, and control purposes.
The PSE provide a computationally attractive framework for the spatial evolution of perturbations in non-parallel and weakly nonlinear scenarios [37][38][39]. They are obtained by removing terms of O(1/Re 2 ) and higher from the NS equations and are significantly more efficient than conventional flow simulations based on the governing equations. In weakly non-parallel flows, e.g., in the pre-transitional boundary layer, flow fluctuations can be separated into slowly and rapidly varying components via the following decomposition for the fluctuation field q = [ u v w p ] T in (1). For a specific spanwise wavenumber and temporal frequency pair (β, ω), we consider q(x, y, z, t) =q(x, y) χ(x, z, t) + complex conjugate, whereq(x, y) and χ(x, z, t) are the shape and phase functions, and α(x) is the streamwise varying generalization of the wavenumber [39]. This decomposition separates slowly (q(x, y)) and rapidly (χ(x, z, t)) varying scales in the streamwise direction. The ansatz in Eq. (2) provides a representation of oscillatory instability waves such as TS waves.
The ambiguity arising from the streamwise variations of bothq and α is resolved by imposing the condition Ωyq * q x dy = 0, whereq * denotes the complex conjugate transpose of the vectorq. In practice, this condition is enforced through the iterative adjustment of the streamwise wavenumber [37,Sec. 3.2.5]. Following the slow-fast decomposition highlighted in (2), the linearized NS equations are parabolized under the assumption that the streamwise variation ofq and α are sufficiently small to neglectq xx , α xx , α xqx , α x /Re 0 , and their higher order derivates with respect to x, resulting in the removal of the dominant ellipticity in the NS equations. The linear PSE take the form where expressions for the operator-valued matrices L and M can be found in [37]. We next propose a two-step modeling procedure to study the dominant mode interactions in weakly-nonlinear mechanisms that arise in spatially evolving flows.

III. PARABOLIZED FLOQUET EQUATIONS
In the transitional boundary layer flow, primary instabilities can cause disturbances to grow to finite amplitudes and get saturated by nonlinearity. Secondary stability analysis examines the asymptotic growth of the resulting modulated state and is based on the linearized NS equations around the modified base flow Here, u 0 denotes the original base flow and u pr represents the primary disturbance field. Sinceū is typically spatially or temporally periodic, Floquet analysis is invoked to identify the spatial structure of exponentially growing fluctuations aroundū. However, such analysis relies on a parallel flow assumption and it does not explicitly account for the spatially growing nature of the base flow. To account for the interactions of fluctuations with spatially growing modified base flowū in a computationally efficient manner, we introduce a framework which draws on Floquet theory to enhance the linear PSE. Our approach allows us to capture the dominant mode interactions in the fluctuating velocity field while accounting for non-parallel effects in the base flow. Starting from a spatially or temporally periodic initial condition the linear PSE can be marched downstream to obtain the primary disturbance field. For example, such an initial condition can be obtained using stability analysis of the two-dimensional Orr-Sommerfeld equation or transient growth analysis of streamwise constant linearized equations (under the locally-parallel base flow assumption). When the periodic solutions to the linear PSE computation are superposed to the Blasius boundary layer profile, the modulated base flow (4) takes the following form Here, T represents the Blasius boundary layer profile, φ 0 = 1, u m and φ m for m = 0 are the shape and phase functions corresponding to various harmonics that constitute flow structures of the primary disturbance field (such as TS waves or streaks), and u * m = u −m . Note that each harmonic u m of the modified base flowū inherits a similar slow-fast structure from PSE (cf. Eq. (2)) in which the phase function φ m is spanwise or streamwise/temporally periodic. For example, when TS waves are superposed to the Blasius boundary layer profile the phase functions φ m are streamwise and temporally periodic; see Sec. IV for details. The evolution of fluctuations around the modulated base flow (5) can be studied using the following expansion which, similar to PSE, involves a decomposition of disturbances into slowly (q n ) and rapidly (χ n ) varying components. Note that we follow classical Floquet decomposition [26,50] in assuming that the phase functions χ n represent various harmonics of the same fundamental frequency/wavenumber as φ m in Eq. (5). As a result of this assumption the evolution of each harmonic mode in q can contribute to the evolution of its neighboring harmonics via the periodicity of the modulated base flow (5). For spanwise periodic modulations to the base flow, a concrete example of the form of the fluctuation field (6) is discussed in Remark 1.
Remark 1 When spanwise-periodic streaks with a fundamental wavenumber β are superposed to the Blasius boundary layer profile, the modified base flow takes the form and the spatial evolution of fluctuations that account for fundamental harmonics (in z) around this modulated base flow profile can be studied using the Fourier expansion Here, α n (x) is the purely imaginary streamwise wavenumber of various harmonics, which can evolve in the streamwise direction similarly to linear PSE. Note that if α n (x) is identical for all harmonics, we recover the nondispersive wavepacket assumed in Floquet stability analysis; see Sec. V for additional details.
Under the assumptions of linear PSE, the dynamics of fluctuations represented by (6) can be studied using the Parabolized Floquet Equations (PFE) The state in (7),q contains all harmonics of q in the periodic direction, i.e., and the operators L F and M F inherit the following bi-infinite structure from the periodicity of the phase functions φ m in the modified base flow (5), Note that the operator L i,j captures the influence of the jth harmonicq j on the dynamics of the ith harmonicq i . In practice, the bi-infinite structures of the state vector and operators in Eq. (7) are truncated to account for the spatial evolution of a finite number of essential modes.
In what follows, we employ the PFE to examine two canonical problems: • the H-type transition scenario (Sec. IV); and • the formation of streamwise elongated streaks in laminar boundary layer flow (Sec. V).
In the wall-normal direction, homogenous Dirichlet boundary conditions are imposed, where L y denotes the height of the computational domain. We discretize differential operators L F and M F using a pseudospectral scheme with N y Chebyshev collocation points in the wall-normal direction [51] and employ an implicit Euler method to march the PFE (7) in the streamwise direction with constant step-size ∆x. The diagonal lines represent the base flows that enter as coefficients into the linear PSE and PFE, respectively. (b) When the secondary disturbances are of the same frequencies and spatial wavenumbers as the primary disturbances, the dominant harmonics resulting from each PFE iteration (e.g.,q0,q±1, andq±2 in the formation of streaks) can be used to update the modulation to the base flow, iterate the PFE, and compute an equilibrium configuration.

A. Two-step modeling procedure
To model the effect of mode interactions in weakly nonlinear regimes we consider the following two-step procedure: 1. The linear PSE are used to march the primary harmonic and obtain the corresponding velocity profile u pr at each streamwise location.
2. The PFE are used to march all harmonicsq and obtain the spatial evolution of velocity fluctuations around the modified base flowū = u 0 + u pr .
The PFE are thus used to study the effect of dominant harmonic interactions on the growth of disturbances in the streamwise direction. The block diagram in Fig. 2(a) illustrates our modeling procedure. When the secondary disturbances contain the same temporal frequency and spatial wavenumbers as the primary disturbances, the dominant harmonics resulting from the PFE (e.g.,q 0 ,q ±1 , andq ±2 in the formation of streaks) can be subsequently used to update the modulation to the base flow, iterate the PFE computation, and thereby provide an equilibrium configuration. The block diagram in Fig. 2(b) illustrates how our framework can be employed to correct primary disturbances that subsequently modulate the base flow. This should be compared and contrasted to the conventional Floquet analysis in which the frequencies and wavenumbers of the identified secondary instabilities are different from those in the primary disturbances. Our computational experiments in Sec. V demonstrate that the flow state of the feedback interconnection in Fig. 2(b) converges after a certain number of iterations. While the equilibrium configuration in Fig. 2(b) implies that our framework is inherently nonlinear, each step in our iterative procedure is linear. We note that a similar iterative method was successfully utilized for model-based design of spanwise wall oscillations in turbulent channel flows [12].

B. Comparison with nonlinear PSE
In contrast to the nonlinear PSE, which treat the interaction between various modes as a nonlinear input forcing, the PFE introduced in (7) account for a subset of dominant interactions between the primary and secondary modes while maintaining the linear progression of the governing equations. The implementation of the PFE is thus less complex than that of the nonlinear PSE as the explicit evaluation of the nonlinear terms commonly used in the PSE and the associated transformations between physical and Fourier space are avoided. In particular, at each downstream location, the PFE eliminate the inner iterations that are required for converging over nonlinear terms and can become costly as the amplitudes of harmonics grow.
In comparison to the nonlinear PSE, the PFE thus reduce both the implementational and computational complexity. While the PSE involve inner iterations to compute nonlinear terms at each step along the longitudinal domain, the PFE are advanced by inverting a high-dimensional sparse matrix. The computational efficiency of both PSE and PFE methods can be significantly improved using parallel computing packages. However, improving the computational efficiency of such methods in a way that would lead to a fair comparison is out of the scope of the current work.

IV. H-TYPE TRANSITION
We next apply our approach to model an H-type transition scenario in a zero-pressure-gradient boundary layer [26]. This route to transition begins with the exponential growth of two-dimensional TS waves which, upon reaching a critical amplitude, become unstable to secondary disturbances. The modulation of the Blasius profile by the TS waves induces the amplification of otherwise stable oblique modes which have half the frequency of the TS waves. While the linear PSE can be used to characterize the spatial evolution of modes arising from primary instabilities, secondary instabilities that trigger the growth of subharmonic modes call for an expansion in the harmonics of the modulated base flow. Such a growth mechanism cannot be identified via the normal-mode ansatz employed in the linear PSE but it can be captured by marching the PFE, a model resulting from a combination of linear PSE with Floquet decomposition.

A. Setup
At the initial position x 0 and for a real-valued temporal frequency ω, the fundamental mode in the H-type scenario is identified as the eigenvector corresponding to the most unstable complex eigenvalue α from the discrete spectrum of the standard two-dimensional Orr-Sommerfeld equation; see [14, Sec. 7.1.2] for additional details. The linear PSE can be used to march this fundamental mode, which is in the form of a two-dimensional TS wave, and obtain a reasonable prediction of its spatial growth [39]. We use the resulting solution to a primary linear PSE computation to augment the Blasius boundary layer profile u 0 in the base flow for the PFE as Equivalently, the base flowū can be written in the following compact form where u 0 is the Blasius profile. In Eq. (9), α r is the real part of the prescribed streamwise wavenumber of the fundamental mode α(x) resulting from linear PSE, c = ω/α r denotes the phase speed of the fundamental and subharmonic modes in the fixed (laboratory) frame, and [ U T (x, y) V T (x, y) 0 ] T denotes the TS wave whose local amplitude and shape are obtained from the linear PSE computation. Note that the exponential growth resulting from the imaginary part of the wavenumber, α i , is absorbed into the amplitudes of U T (x, y) and V T (x, y).
To study the evolution of subharmonic modes that are triggered via secondary instability mechanisms, we follow the Floquet decomposition which was originally conducted in the moving frame [14, Sec. 8.2] using the following Fourier expansion, but in the fixed (laboratory) frame whereq * n =q −n−1 , and γ n and σ n are the spatial and temporal detuning factors corresponding to the nth subharmonic. The imaginary and real parts of γ n (σ n ) denote the spatial (temporal) growth rate and the detuning in the wavenumber (frequency), respectively. In practice, the detuning factors of the wavenumber and frequency are negligible, i.e., γ n and σ n can be assumed to be purely imaginary. A further assumption of nondispersive wavepackets in accordance with classical Floquet analysis [26] brings the ansatz for the fluctuation field to the following form where Gaster's transformation [52] has been used to replace σ with −cγ and all subharmonic modes are assumed to share a uniform detuning parameter γ. In making these approximations, we have followed [26] in assuming that all Fourier components have the same phase speed, which is consistent with experimental studies [27]. While α r in the PFE computation is prescribed by the solution of linear PSE for the primary disturbance field, we update the spatial growth rate −γ via a similar scheme to the one used for the streamwise wavenumber update in PSE [37,Sec. 3.2.5].
The PFE account for the interaction between different subharmonics by leveraging the slow-fast decomposition inherited from the solution of the linear PSE. By substituting the ansatz (11) and the modulated base flow (9) into the linearized NS equations, we arrive at the PFE which take the form of Eq. (7). For this case study, the operators L F and M F in the PFE (7) are provided in Appendix A.
The procedure for obtaining the results presented in the next subsection can be summarized as follows: 1. Solve the spatial eigenvalue problem corresponding to the Orr-Sommerfeld equations to obtain the initial complex wavenumber α(x 0 ) and the shape of initial TS wave u pr (x 0 ).

Use linear PSE to march the TS wave downstream and obtain α(x) and u pr (x).
3. Augment the Blasius boundary layer profile with the TS wave u pr (x) to obtain the modulated base flowū.
4. Obtain the initial uniform growth rate γ and shape functionq as the most unstable eigenvalue and eigenvectors from standard Floquet analysis [26] at x 0 .
5. Use PFE to compute the spatial evolution of all subharmonic modes around the modulated base flow.

B. Growth of subharmonic secondary instabilities
We next examine the interaction of TS waves with subharmonic secondary instabilities in the classical H-type transition scenario. This problem was initially studied in [50] and has been further explored using both experiments [27] and numerical simulations [28,36].
Following Refs. [26,27], the primary disturbance field is generated by marching a TS wave with a root-mean-square (rms) amplitude of 4.8×10 −3 and frequency ω = 0.496 from Re 0 = 424 to Re = 700 using linear PSE. We subsequently initialize the PFE with the most unstable eigen-mode from the classical Floquet analysis [26] to study the growth of subharmonic secondary instabilities triggered by the TS wave. The initial rms amplitude of the subharmonic mode is 1.46 × 10 −5 and its frequency and spanwise wavenumber are ω = 0.248 and β = 0.132, respectively. The initial spatial growth rate −γ is obtained by applying Gaster's transformation to the temporal growth rate. We consider a truncation of the bi-infinite state q with 2N modes, i.e., n = −N · · · N − 1, and a computational domain with L x × L y = 1100 × 40. Our computations demonstrate that N y = 40, ∆x = 30, and N = 2 (i.e., 4 subharmonic modes), provide sufficient accuracy in capturing the physics of H-type transition. Figure 3(a) shows the rms amplitude of individual modes resulting from experiments [27], DNS [28], nonlinear PSE [36], along with the present PFE computations. Here, the modes are denoted by (l, k), where l stands for the temporal frequency of the harmonic and subharmonic modes as a multiple of the subharmonic mode frequency ω = 0.248 and k represents the spanwise wavenumber as a multiple of the fundamental wavenumber β = 0.132. The amplitude of the (1, 1) mode from the PFE is in excellent agreement with other results. While the amplitude of the (3, 1) mode is somewhat under-predicted, the general trend in the growth of this mode is captured by the PFE. We note that in the absence of interactions between harmonics, the linear PSE results in inaccurate predictions for the amplitude of subharmonic modes; see thin solid lines in Fig. 3(a). Figures 3(b-d) show the amplitude of the streamwise velocity profile of the modes considered in this study normalized by the results of nonlinear PSE. For all three modes, the profiles resulting from PFE are in good qualitative agreement with the result of nonlinear PSE.
Remark 2 As mentioned in Sec. III B, even if numerical efficiency is not a decisive factor, the PFE are conceptually easier to implement than the nonlinear PSE and they do not encounter numerical issues that arise in computing the nonlinear terms. We also note that in the H-type transition scenario, the growth of nonlinear terms can impede the convergence of the nonlinear PSE. Specifically, after secondary modes grow to the same order of magnitude as primary modes, the nonlinear terms become dominant and, at some point, update schemes for α(x) fail to converge. There have been previous efforts to suppress the feedback from secondary to primary modes and to maintain the march of nonlinear PSE through the transitional region; see for example [37,Sec. 3.4.3]. The framework proposed in this paper allows for the formal investigation of such effects by limiting the interactions within the PFE framework to a subset of dominant harmonics of the base flow.

V. STREAMWISE ELONGATED LAMINAR STREAKS
Bypass transition often originates from non-modal growth mechanisms that can lead to streamwise elongated streaks; see for example [53]. The streaks can attain substantial amplitudes (15-20% of the free-stream velocity) and make the flow susceptible to the amplification of high frequency secondary instabilities [54,55]. Secondary instability analysis of saturated streaks has been previously used to analyze the breakdown stage in the transition process [30,35]. However, nonlinear effects that influence the formation of streaks become prominent in earlier stages of transition and before the breakdown of streaks. In this section, we utilize the PFE to capture the interactions between various modes in the amplification of the streaks. We focus on the interaction between different spanwise harmonics and study their contribution to the mean flow distortion (MFD), which in turn affects the energy balance among various harmonics that form streaks. We show that the linear PSE fail to predict such a phenomenon and demonstrate how the PFE provide the means to capture the correct trend in the MFD as well as the resulting velocity distribution.

A. Setup
We trigger the formation of streaks by imposing an initial condition computed via the PSE-based optimization approach introduced in [56]. This optimal initial condition yields the highest amplification of perturbation kinetic energy and it is obtained from the singular value decomposition of a pseudo-propagator which advances arbitrary superpositions of the most unstable eigenfunctions in the non-parallel base flow. The initial perturbation field describes a set of counter-rotating streamwise vortices which give rise to the streaks by means of the lift-up mechanism. We compute the spatial evolution of the perturbation field via the linear PSE and use this solution to augment the Blasius boundary layer base flow profile u 0 for the subsequent PFE computations as which can be written in the following compact form In the linear PSE computations, the real part of the complex wavenumber α is set to zero in accordance with the nature of streamwise elongated streaks and its imaginary part is initialized with a small number (e.g., 10 −10 ). Moreover, the exponential growth resulting from the imaginary part of α is absorbed into the amplitude of U S,1 (x, y) in Eq. (12). The velocity field of streamwise elongated streaks is dominated by the growth of the streamwise component while wall-normal and spanwise components experience viscous decay. As a consequence, we disregard the normal and spanwise components of the solution to linear PSE, and only use the streamwise component U S,1 (x, y) in (12), which is also in agreement with the structure and amplitude of the initial condition. We represent the state in the PFE using the Fourier expansion whereq 0 is the MFD, higher-order harmonics in the spanwise direction represent various streaks of wavelength 2π/(nβ), and α = iα i is the uniform streamwise wavenumber over all harmonics. Similar to the procedure in Sec. IV A, we derive the PFE in the form of Eq. (7) by substituting the ansatz (14) and modulated base flow (12) into the linearized NS equations and rearranging the governing equations for each harmonic. For this case study, the operators L F and M F in PFE (7) are provided in Appendix B. The procedure for obtaining the results presented in the next subsection can be summarized as follows: 1. Compute the initial fluctuation field for maximum streamwise growth using the methodology presented in [56].
2. March the initial fluctuation field downstream for the linear evolution of the optimal streak using linear PSE.
3. Augment the Blasius boundary layer profile with the solution to linear PSE to obtain the modulated base flow u according to Eq. (13).
4. Use the same initial condition as the linear PSE in step 2 and its complex conjugate to initializeq ±1 in Eq. (14) and initialize other harmonics with zero. Also, set the initial growth rate α(x 0 ) to a small imaginary number.
5. Use the PFE to compute the spatial evolution of all harmonic modes around the modulated base flow.

B. Nonlinear evolution of optimal streaks
Although the initial condition imposed at the first downstream location only contains a single spanwise wavenumber, the appreciable amplitudes of the developing streaks lead to modal interactions that introduce additional harmonics and modulate the mean flow. To investigate these harmonic interactions, we consider truncations of the bi-infinite stateq in the PFE (7) to 2N + 1 harmonics in z, i.e., n = −N, · · · , N . We set N = 3 and consider a computational domain with L x × L y = 2000 × 60, N y = 80 collocation points in the wall-normal direction, and a step-size of ∆x = 15. Our computational experiments demonstrate that similar predictions can be achieved with appreciably coarser computational grids, e.g., N y = 40 and ∆x = 30.
The temporal frequency, streamwise and fundamental spanwise wavenumbers are set to ω = 0, α = −10 −10 i, and β = 0.4065, respectively. Note that the small imaginary-valued wavenumber α corresponds to infinitely long structures in the streamwise direction that saturate after a particular streamwise location. Moreover, α = 0 maintains a wellconditioned downstream progression for the PFE computations. We initialize the PFE computation at Re 0 = 467 with zero initial conditions for allq n with n = ±1. The fundamental harmonicq ±1 is initialized with the same initial condition as the primary linear PSE computations and with an rms amplitude of 6.4 × 10 −4 . Since this case study considers the evolution of disturbances with a slowly varying streamwise wavenumber α, we set α x = 0 for both the primary linear PSE and the subsequent PFE computations. To verify the predictions of our framework, we also conduct direct numerical simulations of the nonlinear NS equations (with the same initial conditions) using a second-order finite volume code with 2049 × 257 × 257 grid points in the streamwise, wall-normal, and spanwise dimensions, respectively.
As illustrated in Fig. 4, all harmonics undergo an initial algebraic growth followed by saturation. The solution to the linear PSE accurately predicts the evolution of the fundamental spanwise harmonic; cf. Eq. (14). The PFE accurately predict the growth of the dominant harmonics, and especially the MFD. While a discrepancy is observed for the third harmonic, its contribution to the overall structure of the streaks is negligible. The reasonable prediction of growth trends and generation of the MFD component is a direct consequence of accounting for interactions between different harmonics within our framework because, apart fromq ±1 , all other harmonics were initialized with zero. Figure 5 shows the cross-plane spatial structure of the streaks comprised of all harmonics in the spanwise direction at x = 2400. Comparison of Figs. 5(a) and 5(b) indicates a significant discrepancy between the shape of the structures in the cross-plane if the interaction between modes is not taken into account. Since the first (fundamental) harmonic has much larger amplitude than the second and third harmonics, the velocity distribution resulting from the linear PSE is dominated by the structure of the first harmonic. Furthermore, in the absence of interactions between harmonics, the linear PSE would not be able to generate the MFD and would thus result in inaccurate predictions for the amplitude of higher-order harmonics; see dashed lines in Fig. 4. In Fig. 4, we have used a scaled version of the initial condition for the first harmonic to initialize the linear PSE computations for higher-order harmonics. Figure 5(c) demonstrates excellent agreement of the results from PFE and DNS. To emphasize the contribution of the MFD on the final velocity distribution, Fig. 5(d) shows the cross-plane contours of the streamwise velocity component without the MFD. While the influence of the MFD on the growth of the first, second, and third harmonics has been retained, this figure demonstrates its influence on the shape of the streamwise elongated structures.
Since the evolution of streaks is highly influenced by nonlinear interactions, it is worth examining if the results would change with an increase in the streak amplitude. To test the robustness of our framework we consider a different initial condition which has twice the amplitude as the previous case. We use the same computational configuration as before and initialize all harmonics apart fromq ±1 with zero. As shown in Fig. 6(b), the linear PSE provide a poor prediction for the amplitude of the fundamental harmonic. We use the MFDq 0 , fundamental harmonicsq ±1 , and second harmonicsq ±2 from each run of PFE to update the base flow modulation and rerun the PFE (cf. Fig. 2(b)). The base flow in subsequent iterations is thus given by U (x, y) = U B (x, y) + U S,0 (x, y) + U S,1 (x, y) e iβz + U * S,1 (x, y) e −iβz + U S,2 (x, y) e 2iβz + U * S,2 (x, y) e −2iβz , V (x, y) = V B (x, y) + V S,0 (x, y), and takes the compact formū Here, U S,0 , U S,1 , and U S,2 represent the MFD, first-and second-order harmonics corresponding to the solution of the previous PFE run, respectively, and u 0 contains both the Blasius profile and the MFD. We note that higher-order harmonics are omitted from the base flow modulation in Eq. (16) as they do not significantly influence the profiles that result from the iterative PFE procedure. Furthermore, similar to Eq. (12), special care is taken in modulating the base flow, i.e., higher-order harmonics (|m| ≥ 1) are excluded from the wall-normal and spanwise components of the base flow modulation in agreement with the structure and amplitude of the initial condition. Figure 6(a) demonstrates how iterating the PFE can improve our prediction of the predominantly nonlinear streak evolution. The rms curves for various harmonics converge after 7 iterations of the PFE feedback loop illustrated in  . Note that in the previous case of moderate-amplitude streaks, subsequent iterations were not necessary and accurate results were obtained after one run of the PFE over the streamwise domain. The high number of iterations required for convergence is indicative of the significant role nonlinear terms play in the more challenging case of high-amplitude streaks. In Fig. 6(b), we compare the result from the final iteration (solid lines in Fig. 6(a)) with the result of DNS and nonlinear PSE. We see that the PFE capture the initial algebraic growth, inhibition of growth, and general trend in the saturation of amplitudes. Nonlinear interactions generate an appreciable MFD that alters the mean flow profile and hampers the growth of the principal harmonic in comparison to the single mode computation of linear PSE (cf. thin blue line in Fig. 6(a)). Previous studies have reported the stabilizing effect of nonlinearity on the evolution of unsteady streaks [57] and the boundary layer response to perturbations [58,59]. In Fig. 6(b), while a discrepancy is observed in the prediction of the MFD, dominant trends in the shape and amplitude of velocity profiles are consistently captured in the streamwise domain, and the final amplitudes are in close agreement with the result of nonlinear PSE. However, both nonlinear PSE and PFE seem to under-predict the growth of the fundamental harmonic and MFD. Figures 8 and 9 show the streamwise velocity component of the MFD and first harmonic at x = 1700 and x = 2400, which correspond to the largest error in matching the MFD and the end of the longitudinal domain. Finally, as the cross-plane contour plots of Fig. 7 demonstrate, the PFE provide good predictions for the spatial structure of the streaks that are comprised of various harmonics. In the present case study, nonlinear interactions play a crucial role in the growth of high amplitude streaks. The PFE capture the nonlinear interactions by allowing spanwise modulations to the base state and extending the state variableq over multiple harmonics in the spanwise direction. Subsequent iterations of the PFE feedback loop (Fig. 2(b)) refine our predictions of nonlinear interactions on a sweep-by-sweep basis, i.e., by treating the base flow as a streamwise varying parameter in each individual PFE run, and only updating it for the next run. This is in contrast to nonlinear PSE in which nonlinear interactions are captured by explicitly converging over the corresponding nonlinear terms at each step of the streamwise progression. Regardless of how nonlinear interactions are captured, our results demonstrate the difficulty in accurately capturing the correct growth of these optimal streaks (cf. Fig. 6(b)). While the approximation used by the PFE framework may be seen as a limitation, the encouraging performance of the PFE warrants future study into improving the predictive capability of models that capture harmonic interactions through iterative refinement of the base state and not by explicitly computation of nonlinear terms.

VI. CONCLUDING REMARKS
We have combined ideas from Floquet decomposition and the linear PSE to develop the parabolized Floquet equations which can be used to march primary and secondary instability modes while accounting for dominant mode interactions. Our modeling framework involves two steps: (i) the linear PSE are used to march the primary disturbances in the streamwise direction; (ii) the PFE are used to march velocity fluctuations around the modulated base flow profile while capturing weakly nonlinear effects and the interaction of modes. The developed framework can  account for secondary instabilities as fluctuations around a modulated base flow that includes primary modes generated in step (i). The PFE involve a computationally efficient linear march of various harmonics and can be used as a tool to decipher the role of individual harmonics in the spatial evolution of the fluctuation field. Furthermore, subsequent iterations of the PFE, in which the base flow modulation results from the previous PFE computation, can provide a corrective sequence that improve the quality of prediction. To demonstrate the utility of the proposed modeling framework, we have examined the secondary instability analysis of the H-type transition scenario and the evolution of streamwise streaks. Our computational experiments demonstrate good agreement with DNS and nonlinear PSE. We note that the overall performance of the proposed method relies on the reasonable prediction of the evolution of primary disturbances using linear PSE. In cases where the linear PSE give a poor prediction, an additional source of white or colored stochastic excitation can be used to replicate the effect of nonlinearities and improve the outcome of linear PSE; see [44,Sec. IV]. For this purpose, the spatio-temporal spectrum of stochastic excitation sources can be identified using the recently developed theoretical framework outlined in [13,60]. This methodology can also be used to improve the accuracy of results when nonlinear interactions are of critical importance in the evolution of multimodal dynamics. Implementation of such ideas to further improve the predictive capability of the proposed method is a topic for future research.
In the PFE framework, nonlinear interactions are captured via the interplay between the state and periodic base flow. The equilibrium configuration in Fig. 2(b) provides the means to better approximate nonlinear interactions through iterative refinement of the base state. While this configuration implies that the PFE framework is inherently nonlinear, each run of the PFE is linear and is thus well-suited for feedback control design using the tools from linear systems theory. More specifically, the modes that are marched using PFE modulate the base state as a streamwise varying parameter in subsequent PFE runs and thus should not be assumed as variables that violate the premise of linearity. Based on the fluctuation field generated at each sweep of PFE, an optimal control strategy can be synthesized to perturb the dynamics of subsequent PFE runs; see schematic in Fig. 10. While both the control strategy and dynamics are simultaneously updated, convergence of the fluctuation fieldq would insure that the final control design is indeed optimal. Analyzing the performance of this design strategy and providing theoretical justification for convergence calls for additional in-depth examination. For the case study considered in Sec. IV, the operators L n,m and M n,m in L F and M F from Eq. (7) are of the form: Γ n+ = i ((n − 0.5)α r + γ) U * T + V * T ∂ y , Ω n = U B − 2 i Re ((n + 0.5)α r + γ) . (A1g)

Appendix B: Operators for streamwise streaks
For the case study considered in Sec. V, the operators L n,m and M n,m in L F and M F from Eq. (7) are of the form: where In subsequent PFE runs, the second harmonic and MFD also augment the base flow (cf. Eq. (15)) and U B and V B in Eqs. (B1a) and (B1c) denote the combination of the Blasius profile and the MFD. Note that consistent with nonlinear PSE, the boundary conditions on the MFD require special treatment; see [39]. In addition, the contribution of the pressure gradient to the streamwise velocity momentum is removed, i.e., M 0,0 (1, 4) = 0. Finally, interactions with the second harmonic generated from previous PFE runs are facilitated by the off-diagonal operators: