1. Introduction
The boundary layer transition is of great significance in both flow physics and engineering applications. Many associated questions, such as what physical effect is dominant in the transition and how a flow breaks down to turbulence, are not fully understood. It is currently known that the receptivity, transient growth, eigenmode growth, parametric instability and resonance and bypass mechanisms can be involved in the transition process (Morkovin, Reshotko & Herbert Reference Morkovin, Reshotko and Herbert1994). Early knowledge of the boundary layer transition was partly provided by linear eigenmode analysis of Tollmien–Schlichting instability waves. However, in hypersonic states, a family of additional higher modes of acoustic nature were found to be important by linear stability theory (LST) with the parallel-flow assumption (Mack Reference Mack1965, Reference Mack1984). Particularly, the Mack second mode has received much attention due to its high growth rate if the Mach number exceeds approximately 4.
Despite the success of linear stability analysis, our understanding of the mechanisms of the boundary layer transition is still insufficient. Although the second mode appears to be preponderant in the linear instability stage, the upstream receptivity and transient growth as well as the downstream nonlinear stage are less considered jointly. In the flow past a blunt body, the initial amplitude of the first-mode waveband is potentially able to overtake that of the second mode if the receptivities across the bow shock wave (Fedorov et al. Reference Fedorov, Ryzhov, Soudakov and Utyuzhnikov2013), entropy layer (Goparaju & Gaitonde Reference Goparaju and Gaitonde2021) or roughness (Kuester & White Reference Kuester and White2015) are considered. Meanwhile, the integrated N-factor of the first mode can exceed that of the second mode (Masad Reference Masad1993). Direct numerical simulation (DNS) on a Mach 15 flow over a parabola by Zhong (Reference Zhong2001) showed that the maximum receptivity coefficient of the first mode was over twice that of the second mode in a wide frequency range of the forcing. In experimental research, fluctuations in a conventional wind tunnel are usually the most energetic in a relatively low-frequency range (Schneider Reference Schneider2001). A recent comparative study by Hildebrand et al. (Reference Hildebrand, Choudhari, Deegan, Huang and Duan2022) of various hypersonic quiet or noisy wind-tunnel environments showed that the power spectral density of background noise is the largest at a frequency f = O(1) kHz and exponentially decays at higher frequencies. These findings indicate that the first mode may not be weak in the vicinity of the leading edge. Therefore, due solely to the upstream effect, it is already difficult to determine beforehand whether the first or the second mode is dominant in various geometric configurations and free-stream conditions.
Considering the downstream nonlinear instability stage, the complete breakdown mechanisms and relative significance of the first and second modes in hypersonic states remain unclear. The complicated parametric resonance in a real environment is usually decomposed into individual regimes, including fundamental, subharmonic, oblique and combination resonance or breakdown. In a controlled hypersonic transition, the dominant primary wave is usually assumed or found to be the two-dimensional (2-D) second mode. Accordingly, numerous experimental (Chokani Reference Chokani1999; Kennedy et al. Reference Kennedy, Laurence, Smith and Marineau2018a, Reference Kennedy, Jewell, Paredes and Laurence2022) and DNS (Sivasubramanian & Fasel Reference Sivasubramanian and Fasel2015; Hader & Fasel Reference Hader and Fasel2019; Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2020) studies found that second-mode fundamental breakdown, with a frequency of hundreds of kilohertz, was dominant. In these observations, streamwise hot streaks with heat transfer overshoot and evident interactions between the second mode and its higher harmonics were common phenomena. Under the condition of a Boeing/AFOSR Mach 6 quiet tunnel, the second-mode fundamental breakdown was also significant in transitions induced by a wave packet (Sivasubramanian & Fasel Reference Sivasubramanian and Fasel2014) or random forcing (Hader & Fasel Reference Hader and Fasel2018). As a lesser mechanism, second-mode subharmonic resonance was also detected (Bountin, Shiplyuk & Maslov Reference Bountin, Shiplyuk and Maslov2008; Hader & Fasel Reference Hader and Fasel2019; Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2020). The second-mode oblique breakdown was numerically investigated by Franko & Lele (Reference Franko and Lele2013) and Hartman, Hader & Fasel (Reference Hartman, Hader and Fasel2021) but scarcely reported by experiments.
In other existing situations, however, low-frequency bands not belonging to the second-mode instability can be considerably energetic and participate in the nonlinear behaviour of multiple instability waves. These low frequencies could correspond to the first mode, Görtler vortices, wind-tunnel noise or the product of nonlinear interactions between high-frequency side bands. In addition to conventional ones such as the Görtler secondary instability (Ren & Fu Reference Ren and Fu2015; Chen, Huang & Lee Reference Chen, Huang and Lee2019) and first-mode oblique breakdown (Mayer, Wernz & Fasel Reference Mayer, Wernz and Fasel2007; Guo et al. Reference Guo, Shi, Gao, Jiang, Lee and Wen2022b), new breakdown scenarios involving these low-frequency components are possible. For example, a 20–30 kHz component was pronounced and even dominant on a flared cone in the Peking University Mach 6 quiet wind tunnel (Zhu et al. Reference Zhu, Zhang, Chen, Yuan, Wu, Chen, Lee and Elhak2016; Xiong et al. Reference Xiong, Yu, Lin, Zhao and Wu2020). To address this problem, Chen, Zhu & Lee (Reference Chen, Zhu and Lee2017) employed the nonlinear parabolised stability equation (NPSE) to revisit the interaction between the low-frequency and the second mode, and indicated the possibility of combination resonance. However, further support for the corresponding experiments or full three-dimensional (3-D) DNS was not provided. Zhu et al. (Reference Zhu, Zhu, Gu, Lee and Smith2022) also experimentally demonstrated that a nearly adiabatic wall or heated wall enabled the first mode to be comparable to or overtake the second mode in frequency spectra. Under other wind-tunnel conditions, solid experimental evidence of interactions between the second mode and the low-frequency component was presented by Kimmel & Kendall (Reference Kimmel and Kendall1991), Bountin et al. (Reference Bountin, Shiplyuk and Maslov2008), Munoz, Heitmann & Radespiel (Reference Munoz, Heitmann and Radespiel2014) and Craig et al. (Reference Craig, Humble, Hofferth and Saric2019). Generally, however, the second mode was considerably stronger than the low-frequency component in the primary instability stage. The present paper serves to show another possibility that the low-frequency optimal disturbance can be as important as the optimal planar wave, which evolves into the second mode downstream. Furthermore, the present study aims to reveal the breakdown scenario and competition mechanism if the low-frequency component and the high-frequency second mode are equally significant and both result in primary instabilities.
In addition to the long-term discussions and debates on the first and second modes, non-modal instabilities have also attracted attention. Here, the terminologies ‘modal’ and ‘non-modal’ in the present paper denote the long-term characteristic of the eigenvalue problem for normal modes and the short-term characteristic of the initial value problem in local analyses, respectively (Schmid Reference Schmid2007; Kerswell Reference Kerswell2018). With the appearance of numerous relevant studies in the 1990s, the physical image of transition had to be reconsidered. For instance, Trefethen et al. (Reference Trefethen, Trefethen, Reddy and Driscoll1993) claimed that although eigenmodes of the linear Navier–Stokes (N–S) equation decayed monotonically, the non-modality of the operator led to an algebraic transient growth of the disturbances by factors of O(105). Hanifi, Schmid & Henningson (Reference Hanifi, Schmid and Henningson1996) identified the ‘lift-up’ mechanism of non-modal growth in compressible boundary layers. Since then, remarkable accomplishments have been witnessed in linear and nonlinear non-modal theories (Schmid Reference Schmid2007; Kerswell Reference Kerswell2018). One of the effective analytical tools is resolvent analysis, which seeks the maxima of the optimal gain in the frequency or wavenumber space. This approach has been intensively applied to identify the combination of the most responsive forcing and most amplified state in dynamic systems of wall-bounded shear flows (Ahmadi et al. Reference Ahmadi, Valmorbida, Gayme and Papachristodoulou2019; Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and Mckeon2021; Martini et al. Reference Martini, Rodríguez, Towne and Cavalieri2021; Rigas, Sipp & Colonius Reference Rigas, Sipp and Colonius2021). Other recently popular analysis tools for optimal disturbances in spectral space are represented by the direct-adjoint PSE (Paredes et al. Reference Paredes, Choudhari, Li and Chang2016) and the one-way N–S equation (Towne et al. Reference Towne, Rigas, Kamal, Pickering and Colonius2022) approaches. These spatial marching methods require less memory usage than resolvent analysis based on matrix decomposition, which show advantages in cases with large numbers of grid nodes. The ‘one-way’ method further improves the robustness and accuracy for complicated problems compared with the PSE.
Based on the aforementioned novel methods, non-modal instabilities have been recognised to provide a crucial route to the high-speed boundary layer transition. Paredes, Choudhari & Li (Reference Paredes, Choudhari and Li2017, Reference Paredes, Choudhari and Li2019) applied the NPSE to show that the optimal streak with proper spacing stabilised the first and second modes. Thus, interactions between modal and non-modal instabilities are of interest. Non-modal instability was also found to increase with higher Mach number (Tempelmann, Hanifi & Henningson Reference Tempelmann, Hanifi and Henningson2012). In addition, non-modal instability is capable of evolving into the only dominator. Both NPSE and DNS (Paredes, Choudhari & Li Reference Paredes, Choudhari and Li2020) as well as experimental (Kennedy et al. Reference Kennedy, Jewell, Paredes and Laurence2022) studies found that large nose bluntness could suppress the second mode and exhibit instabilities dominated by non-modal disturbances under flight conditions. The present paper is motivated to study the complete dynamics of multiple optimal disturbances in the linear, nonlinear and transitional regimes. The characterisation of early non-modal and subsequent modal instabilities of optimal disturbances is naturally included in a conventional framework of resolvent analysis and the PSE. The 3-D DNS is utilised to study the multiple instabilities in the late stage. The article is organised as follows. A description of flow conditions and numerical and theoretical tools is given in § 2. The computational set-up for the DNS, including the perturbation strategy of external forcings, is detailed in § 3. Results and discussions on the linear and nonlinear development of optimal disturbances are presented in § 4. Concluding remarks are provided in § 5.
2. Flow conditions and methodology
The flow condition of the Ludwieg-type short-duration wind tunnel Tranzit-M by Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013) is considered here, whose hypersonic boundary layer stability and control has been widely studied. The free-stream conditions are given as follows: Mach number M ∞ = 6.0, stagnation temperature T 0 = 354.08 K, static temperature T ∞ = 43.18 K and unit Reynolds number Re 1∞ = 1.05 × 107 m−1. The wall is assumed to be isothermal with a temperature of Tw = 293 K. The corresponding ratio of the wall temperature to the laminar-flow recovery value is Tw/Trec = 0.954.
The flat-plate model for DNS has a length of Lx = 0.5 m with a sharp leading edge. A Cartesian coordinate system is constructed with the origin at the leading edge, the x direction along the direction of the free-stream velocity, the y direction perpendicular to the flat plate and the z direction satisfying the right-hand rule. The size of the effective computational domain for 3-D DNS is Lx × Ly × Lz = 0.5 m × 0.1 m × 0.01
${\rm \pi}$ m, while additional sponge zones are placed near the outflow boundary to minimise the reflection of disturbances (Mani Reference Mani2012). The spanwise length of the computational domain Lz is designed to be exactly four times the spanwise wavelength of the linearly optimal oblique wave for a convenient Fourier analysis, as discussed later.
2.1. Direct numerical simulation and base-flow solver
The 3-D compressible N–S equations in the Cartesian coordinate system are written in the following conservation form:
\begin{equation}\frac{{\partial {\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}}}}{{\partial z}} = \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}},\end{equation}
where Q = (ρ, ρu, ρv, ρw, ρE)T denotes the vector of conservative variables and F, G and H are the vectors of inviscid and viscous fluxes. Detailed expressions can be found in Anderson (Reference Anderson1995). Here, ρ is the density; u, v and w are the velocity components in the x, y and z directions, respectively; and E is the total energy per unit mass. The right-hand side forcing vector 
$\boldsymbol{f}^{\prime} = (\,{f^{\prime}_1},{f^{\prime}_2},{f^{\prime}_3},{f^{\prime}_4},{f^{\prime}_5})$ in (2.1) is introduced to trigger the boundary layer transition in unsteady simulations. Specifically in the present paper, the operator 
$\boldsymbol{\mathcal{B}}$ constrains f′ to the wall-normal profile only at x 0 = 0.04 m, which looks like a Dirac delta function in the streamwise direction. This forcing vector is assumed to possess a harmonic form with a small amplitude, whose details are determined by resolvent analysis later. The assumption of a calorically perfect Newtonian fluid is made with a constant specific heat ratio of γ = 1.4. Sutherland's law is adopted to calculate the dynamic viscosity μ, while the thermal conductivity κ is computed with a constant Prandtl number of 0.72. In addition, the vector Q is assumed to be decomposed into a 2-D steady base flow and a small disturbance as
\begin{equation}{\boldsymbol{Q}}(x,y,z,t) = {{\boldsymbol{Q}}_b}(x,y) + {\boldsymbol{Q}^{\prime}}(x,y,z,t),\end{equation}
where the subscript ‘b’ and the prime denote base-flow and perturbation variables, respectively. The laminar base flow, satisfying the 2-D N–S equation, is the initial condition of the 3-D DNS and for stability analyses. Meanwhile, the steady base flow is confirmed to converge via the same numerical scheme and mesh as those of the 3-D DNS. Below, unless otherwise stated, the local base-flow and perturbed primitive variables are non-dimensionalised by the corresponding free-stream base-flow quantities, except that pressure p is by 
${\rho _\infty }u_\infty ^2$. The reference length scale Lref = 1 mm is used for non-dimensionalisation.
Both the full 3-D DNS and laminar base-flow simulation are performed using an in-house multi-block parallel finite-volume solver called PHAROS (Hao, Wang & Lee Reference Hao, Wang and Lee2016; Hao & Wen Reference Hao and Wen2020; Hao et al. Reference Hao, Cao, Wen and Olivier2021). This solver has successfully resolved the 3-D instability problem of a hypersonic boundary layer over a double cone via DNS, obtaining excellent agreement with experimental results (Hao et al. Reference Hao, Fan, Cao and Wen2022). The inviscid fluxes are calculated using the modified Steger–Warming scheme (MacCormack Reference Maccormack2014), which is extended to a higher order using the seventh-order weighted essentially non-oscillatory scheme (Jiang & Shu Reference Jiang and Shu1996). The viscous fluxes are discretised by the second-order central difference scheme. The three-stage third-order total variation diminishing Runge–Kutta method is employed for time marching. The boundary conditions are given as follows: free-stream conditions are imposed on the far-field boundary, extrapolation is used on the outflow boundary and isothermal, no-slip and no-penetration conditions are employed on the wall boundary. For 3-D DNS, a periodic condition is used on the spanwise boundary.
2.2. Resolvent analysis
Starting from the linearised N–S equations with a specified base flow, resolvent analysis aims at finding the most amplified disturbances in the desired parametric and spatial domains. Throughout the paper, the terminology ‘global resolvent analysis’ refers to the recently popular framework of resolvent analysis (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019), which resolves the disturbance shape on the x–y plane globally instead of at a ‘local’ station. A more detailed description of the algorithm can be found in the work of Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). Substituting (2.2) into (2.1) and then subtracting the base-flow equation yields the following form:
\begin{equation}
\frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} + \frac{{\partial
{\boldsymbol{F}^{\prime}}}}{{\partial x}} + \frac{{\partial
{\boldsymbol{G}^{\prime}}}}{{\partial y}} + \frac{{\partial
{\boldsymbol{H}^{\prime}}}}{{\partial z}} =
{\boldsymbol{N}}^{\prime} +
\boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}},
\end{equation}where N′ is the nonlinear higher-order term. The complete form of (2.3) can be found in Paredes (Reference Paredes2014). Dropping the higher-order term of (2.3) yields the linearised N–S equations as
\begin{equation}\frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}^{\prime}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}^{\prime}}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}^{\prime}}}}{{\partial z}} = \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}}.\end{equation}Equation (2.4) can be reformulated as
\begin{equation}\frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} = \boldsymbol{\mathcal{A}}{\boldsymbol{Q}^{\prime}} + \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}},\end{equation}
where 
$\boldsymbol{\mathcal{A}}$ is the linearised N–S operator related to the Jacobian matrices.
Bi-Fourier transforms can be performed in temporal and spanwise directions, which allows a solution Q′ to be written in the following form:
\begin{equation}{\boldsymbol{Q}^{\prime}}(x,y,z,t) = \hat{{\boldsymbol{Q}}}(x,y)\,\textrm{exp(i}\beta z - \textrm{i}\omega t) + \textrm{c}\textrm{.c}\textrm{.}\end{equation}
Here, 
$\hat{{\boldsymbol{Q}}}$ is the complex eigenfunction, β is the spanwise wavenumber, ω is the real angular frequency and c.c. denotes the complex conjugate. Since Q′ is the linear response of the forcing, f′ can be written in a similar form:
\begin{equation}{\boldsymbol{f}^{\prime}}(x,y,z,t) = \hat{{\boldsymbol{f}}}(x,y)\,\textrm{exp(i}\beta z - \textrm{i}\omega t) + \textrm{c}\textrm{.c}\textrm{.}\end{equation}Substituting (2.6) and (2.7) into (2.5) gives
\begin{equation}\hat{{\boldsymbol{Q}}} = \boldsymbol{\mathcal{RB}}\hat{{\boldsymbol{f}}},\quad \boldsymbol{\mathcal{R}} = {( - \textrm{i}\omega \boldsymbol{\mathcal{I}} - \boldsymbol{\mathcal{A}})^{ - 1}},\end{equation}
which reflects the relationship between the forcing and its linear response (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019). Here, 
$\boldsymbol{\mathcal{I}}$ is the identity operator.
In resolvent analysis, the purpose is to seek the forcing and response that maximise the energy amplification, i.e. the optimal gain σ defined by
\begin{equation}{\sigma ^2}(\beta ,\omega ) = \mathop {\max }\limits_{\hat{{\boldsymbol{f}}}} \frac{{{{\|{\hat{{\boldsymbol{Q}}}} \|}_E}}}{{{{\|{\boldsymbol{\mathcal{B}}\hat{{\boldsymbol{f}}}} \|}_E}}}.\end{equation}Here, Chu's energy (Chu Reference Chu1965) is used to calculate the energy norm as
\begin{equation}{\|{\hat{{\boldsymbol{Q}}}} \|_E} = \iint_\varOmega {{{\boldsymbol{Q}}^ \ast }\boldsymbol{\mathcal{M}}\hat{{\boldsymbol{Q}}}\,\textrm{d}\kern 0.06em x\,\textrm{d}y} ,\end{equation}
where the asterisk denotes the conjugate transpose, 
$\varOmega$ is the domain for the resolvent analysis and 
$\boldsymbol{\mathcal{M}}$ is the weight operator given by Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). Furthermore, dimensionless Chu's energy containing primitive variables serves as an important indicator of the disturbance evolution in the present study. In detail, the form of (2.10) at a local station rather than an overall integral is given by
\begin{gather}{E_{Chu}}(x) = \frac{1}{2}\int_0^\infty {\left[ {\bar{\rho }({{u^{\prime}}^2} + {{v^{\prime}}^2} + {{w^{\prime}}^2}) + \frac{{\bar{T}}}{{\gamma Ma_\infty^2\bar{\rho }}}{{\rho^{\prime}}^2} + \frac{{\bar{\rho }}}{{\gamma (\gamma - 1)Ma_\infty^2\bar{T}}}{{T^{\prime}}^2}} \right]\textrm{d}y} ,\end{gather}where all the quantities are dimensionless, as stated in § 2.1. This definition consists of the fluctuation kinetic energy and a positive definite thermodynamic energy, which is commonly adopted in non-modal stability analyses. In 3-D DNS, the definition of EChu is replaced by the maximum of (2.11) in the spanwise direction.
As shown by Sipp & Marquet (Reference Sipp and Marquet2013), Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019) and Dwivedi et al. (Reference Dwivedi, Sidharth, Nichols, Candler and Jovanovic2019), the optimisation problem (2.9) can be transformed into an eigenvalue problem as
\begin{equation}{\boldsymbol{\mathcal{B}}^\ast }{\boldsymbol{\mathcal{M}}^{ - 1}}{\boldsymbol{\mathcal{R}}^ \ast }\boldsymbol{\mathcal{MRB}}\,\hat{{\boldsymbol{f}}} = {\sigma ^2}\hat{{\boldsymbol{f}}}.\end{equation}
The operator 
$\boldsymbol{\mathcal{A}}$ is discretised utilising the same methods as for (2.1) except that the inviscid fluxes are computed by the modified Steger–Warming scheme near discontinuities detected by a modified Ducros shock sensor (Hendrickson, Kartha & Candler Reference Hendrickson, Kartha and Candler2018) and a central scheme in smooth regions. The settings of the boundary conditions are identical to those in DNS. The resulting discrete eigenvalue problem is solved by ARPACK software for a given β and ω of the regular mode (Sorensen et al. Reference Sorensen, Lehoucq, Yang and Maschhoff1996). Specifically, 
$\boldsymbol{\mathcal{R}}$ and its conjugate transpose are computed using the SuperLU library (Li et al. Reference Li, Demmel, Gilbert, Grigori, Shao and Yamazaki1999). The largest eigenvalue σ corresponds to the optimal gain, and the associated eigenfunction represents the optimal forcing. The optimal response can be obtained via (2.8). Prior grid convergence of resolvent analysis is conducted, and the converged results are further compared with the DNS and PSE results. More details about the resolvent analysis solver and the associated validation cases are described in Hao et al. (Reference Hao, Cao, Guo and Wen2023).
2.3. Parabolised stability equation analysis
Driven by the need to develop an efficient and accurate tool for stability analyses, Herbert & Bertolotti (Reference Herbert and Bertolotti1987), Herbert (Reference Herbert1988) and Bertolotti (Reference Bertolotti1991) proposed the fundamental idea of linear PSE (LPSE) and NPSE. Within the scope of convective instabilities, they assumed that the modal shape function 
$\hat{\boldsymbol{\psi} }$ of the disturbance is streamwise slowly varying, while the rapidly distorted part is absorbed into the exponential wavefunction. In detail, the disturbance is decomposed into the following form with a finite truncation of series:
\begin{equation}{\boldsymbol{\psi} ^{\prime}}(x,y,z,t) = \sum\limits_{m ={-} M}^M {\sum\limits_{n ={-} N}^N {{{\hat{{\boldsymbol{\psi} }}}_{\!mn}}(x,y)\,\textrm{exp}\left( {\textrm{i}\int_{{x_{0,\textit{P}}}}^x {{\alpha_{mn}}\,\textrm{d}\kern 0.06em\tilde{x}} + \textrm{i}n{\beta_0}z - \textrm{i}m{\omega_0}t} \right)} ,}\end{equation}
where the vector ψ = (ρ, u, v, w, T)T and the superscript ‘T’ represents the transpose, 
${\hat{\boldsymbol{\psi}}_{\!mn}}$ and αmn are the shape function and complex streamwise wavenumber of the Fourier mode (m, n), respectively, and β 0 and ω 0 are the specified fundamental spanwise wavenumber and angular frequency, respectively. Moreover, x 0,P is the distance from the inflow station of the PSE to the leading edge, and M and N are the modal numbers kept in the truncated Fourier series. Specially, mode (0, 0) is physically the mean flow distortion (MFD). Mode (0, n) with a non-zero n usually corresponds to a longitudinal vortex mode, which can evolve into a streak downstream with the streamwise vorticity gradually damped.
Substituting (2.13) into (2.3), dropping the forcing vector and neglecting the terms of 
$O(1/Re_0^2)$, where Re 0 = ρ∞u∞l 0/
$\mu_{\infty}$ and l 0 = (
$\mu_{\infty}$x 0,P/ρ∞u∞)1/2, finally gives rise to the NPSE as
\begin{equation}({\boldsymbol{\mathcal{L}}_0} + {\boldsymbol{\mathcal{L}}_1}){\hat{\boldsymbol{\psi}}_{\!mn}} + {\boldsymbol{\mathcal{L}}_2}\frac{{\partial {{\hat{\boldsymbol{\psi}}}_{mn}}}}{{\partial x}} + \frac{{d{\alpha _{mn}}}}{{\partial x}}{\boldsymbol{\mathcal{L}}_3}{\hat{\boldsymbol{\psi}}_{\!mn}} = {\boldsymbol{\mathcal{F}}_{mn}}.\end{equation}
The base-flow-related operators 
${\boldsymbol{\mathcal{L}}_0}$, 
${\boldsymbol{\mathcal{L}}_1}$, 
${\boldsymbol{\mathcal{L}}_2}$ and 
${\boldsymbol{\mathcal{L}}_3}$ arise from the effects of locally parallel flow, non-parallel base flow, non-local shape function and non-local streamwise wavenumber, respectively, and 
${\boldsymbol{\mathcal{F}}_{mn}}$ is the nonlinear term. Dropping 
${\boldsymbol{\mathcal{F}}_{mn}}$ in (2.14) yields the LPSE, and then only maintaining the local operator 
${\mathcal{L}_0}$ will generate an eigenvalue problem corresponding to LST. Detailed expressions of (2.14) are given by Paredes (Reference Paredes2014). The boundary condition of the shape function is given by the common Dirichlet type as
\begin{equation}{\hat{u}_{mn}},{\hat{v}_{mn}},{\hat{w}_{mn}},{\hat{T}_{mn}} = 0,\quad \textrm{at}\ y = 0\ \textrm{and}\ y \to \infty ,\end{equation}except for the wall-normal velocity for the MFD at infinity:
\begin{equation}\textrm{d}{\hat{v}_{00}}/\textrm{d}y = 0,\quad \textrm{at}\ y \to \infty ,\end{equation}to include the growth effect of the boundary layer thickness arising from the MFD. To close the physical problem, an iterative scheme is used for the streamwise wavenumber of non-MFD modes during spatial marching as
\begin{equation}{\alpha _{mn,new}} = {\alpha _{mn,old}} - \textrm{i}\frac{1}{{{E_{mn}}}}\int_0^\infty {\bar{\rho }\left( {\hat{u}_{mn}^ \ast \frac{{\partial {{\hat{u}}_{mn}}}}{{\partial x}} + \hat{v}_{mn}^ \ast \frac{{\partial {{\hat{v}}_{mn}}}}{{\partial x}} + \hat{w}_{mn}^ \ast \frac{{\partial {{\hat{w}}_{mn}}}}{{\partial x}}} \right)\textrm{d}y} ,\end{equation}where
\begin{equation}{E_{mn}} = \int_0^\infty {\bar{\rho }({{|{{{\hat{u}}_{mn}}} |}^2} + {{|{{{\hat{v}}_{mn}}} |}^2} + {{|{{{\hat{w}}_{mn}}} |}^2})\,\textrm{d}y} .\end{equation}The streamwise wavenumber of the MFD is not updated, however, to improve the numerical robustness. This indicates that the growth of the MFD is entirely reflected in the shape function. The iteration procedure (2.17) and algebraic solution of (2.14) continue until the residual norm of αmn is less than 10−7. Subsequently, the procedure is marched to the next streamwise station until the specified end boundary.
Regarding the numerical implementation, our in-house code CHASES is employed, which integrates the LST, LPSE, NPSE and sensitivity analyses. A series of validation cases for LST and LPSE have been carefully compared with both theoretical (Guo et al. Reference Guo, Gao, Jiang and Lee2020, Reference Guo, Gao, Jiang and Lee2021, Reference Guo, Shi, Gao, Jiang, Lee and Wen2022a) and DNS (Guo et al. Reference Guo, Shi, Gao, Jiang, Lee and Wen2022b; Cao et al. Reference Cao, Hao, Guo, Wen and Klioutchnikov2023) results. The validation of the NPSE code is shown in Appendix A. In terms of the numerical scheme for PSE, the second-order backward Euler scheme is used for streamwise marching, while the fourth-order central difference method is applied for the wall-normal difference. The Vigneron technique is adopted to eliminate the numerical instability resulting from residual ellipticity (Vigneron, Tannehill & Rakich Reference Vigneron, Tannehill and Rakich1978). In the LST analysis which also provides the initial condition for the PSE, the global algorithm via the Chebyshev pseudo-spectral method and a local algorithm using the fourth-order compact finite difference are employed for wall-normal discretisation (Malik Reference Malik1990). The LST eigenvalue problem is solved by the Intel math kernel library.
The presented PSE results have been examined by grid convergence studies. The initial conditions of the PSE are mostly optimal responses obtained from resolvent analysis at x 0,P = 0.045 m, downstream of x 0 = 0.04 m where the optimal forcings are localised in resolvent analysis and DNS. Additionally, LST profiles are also used for initialisation in § 4.1 for comparison. Those newly generated Fourier modes due to nonlinear interactions are assumed to be initially zero in § 4.2 or imposed by Fourier-transformed DNS profiles to improve the accuracy in § 4.3.
3. Case description and simulation strategy
Provided the base flow is prepared, the research procedure starts from a global resolvent analysis to identify all local maxima of optimal gain, which correspond to multiple optimal disturbances. Subsequently, the profiles of optimal disturbances or forcings serve to initiate the PSE and 3-D DNS to reveal the linear and nonlinear behaviour of disturbances.
The grid size for the main computational domain of the 3-D DNS is 4001 × 251 × 97, whose x–y mesh is entirely identical to that of the 2-D steady simulation. The orthogonal structured mesh is uniformly distributed in the streamwise and spanwise directions. In the wall-normal direction, the mesh is stretched with a minimum spacing of 2 × 10−5 m on the wall, and the transitional and turbulent boundary layers are ensured to contain 150 to 210 grid cells inside. The dimensionless grid spacings are Δx + ≈ 1.88, 
$\Delta y_{min}^ + \approx 0.30$ and Δz + ≈ 4.92, evaluated by the friction velocity and wall kinematic viscosity in the fully developed turbulent region x = 0.37 m via either 3D-DNS statistics or the skin friction by the van Driest II formula for Cf (Franko & Lele Reference Franko and Lele2013; Guo et al. Reference Guo, Shi, Gao, Jiang, Lee and Wen2022b). In the linear and early nonlinear stages, the mesh contains at least 45 grid cells for each wavelength of the optimal planar wave (also the second mode downstream) in the x direction and 24 cells exactly for each wavelength of the optimal oblique wave (also the first mode downstream) in the z direction. Grid convergence of the amplitude evolution of the optimal disturbances has been confirmed in these instability stages.
Figure 1 displays the contours of the base-flow velocity and density fields. The dashed line, corresponding to x 0 = 0.04 m, represents the location where the forcings are added in the resolvent analysis and 3-D DNS. This type of forcing appears to be a Dirac delta function looking in the streamwise direction and a constantly activated perturbation in the wall-normal, spanwise and temporal dimensions.
Figure 1. Contours of (a) streamwise velocity and (b) density of the laminar base flow non-dimensionalised by the free-stream value. Dashed line: location where the optimal forcings are imposed on the base flow.
In the DNS, the optimal forcings f′ are directly added to the right-hand side of the N–S equations only at x 0 = 0.04 m, which resembles those of the resolvent analysis. The DNS forcings consist of those corresponding to the optimal planar wave 
${\boldsymbol{f}^{\prime\!\!}_p}$, the optimal oblique wave 
${\boldsymbol{f}^{\prime\!\!}_o}$ and an additional background noise term 
${\boldsymbol{f}^{\prime\!\!}_n}$. The resolvent analysis of these optimal disturbances is detailed in § 4.1, and the resulting parametric set-up for the perturbation strategy in DNS is briefly introduced here.
The mathematical form of the forcing in DNS is given by
\begin{gather}{\boldsymbol{f}^{\prime}}({x_0},y,z,t) = {{\boldsymbol{f}^{\prime}\!\!}_p}({x_0},y,z,t) + {{\boldsymbol{f}^{\prime}\!\!}_o}({x_0},y,z,t) + {{\boldsymbol{f}^{\prime}\!\!}_n}({x_0},y,z,t),\end{gather}
\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_p}({x_0},y,z,t) = {\epsilon_1}\,{\hat{{\boldsymbol{f}}}_{\!p}}({x_0},y)\,\textrm{exp(i}{\beta _p}z - \textrm{i}{\omega _p}t) + \textrm{c}\textrm{.c}\textrm{.,}\end{gather}
\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_o}({x_0},y,z,t) = {\epsilon_2}[\,{\hat{{\boldsymbol{f}}}_{\!o}}({x_0},y)\,\textrm{exp(i}{\beta _o}z - \textrm{i}{\omega _o}t) + {\hat{{\boldsymbol{f}}}_{ \!- o}}({x_0},y){\rm exp} ( - \textrm{i}{\beta _o}z - \textrm{i}{\omega _o}t)] + \textrm{c}\textrm{.c}\textrm{.,}\end{gather}
\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_n}({x_0},y,z,t) = {\epsilon_3}[2r(y,z,t) - 1].\end{gather}
The complete description of the parameters in (3.1) is provided later, and a general introduction to the simulation cases is given first. Settings of DNS cases for comparative studies are listed in table 1. Case A1 is considered as the baseline one, and case A2 is introduced to include the background noise effect in real environments, i.e. the term 
${\boldsymbol{f}^{\prime\!\!}_n}$ is non-zero in (3.1). The noise term 
${\boldsymbol{f}^{\prime\!\!}_n}$ has an identical form to the forcing f′ as a right-hand side term of N–S equations. For the baseline case A1 and cases A3–A6, the amplitude 
${\epsilon_3}$ of the term 
${\boldsymbol{f}^{\prime\!\!}_n}$ is set to zero. In comparison, the effect of background noise is involved in case A2, where 
${\epsilon_3} = 0.1$. In (3.1d), 
$r \in [0,1]$ is the white noise signal provided by a random number generator. The chosen value of 
${\epsilon_3}$ makes the norm of 
${\boldsymbol{f}^{\prime\!\!}_n}$ at least one order of magnitude larger than those of 
${\boldsymbol{f}^{\prime\!\!}_p}$ and 
${\boldsymbol{f}^{\prime\!\!}_o}$ simultaneously. This set-up allows adequate selection of additional potentially unstable frequencies and wavenumbers and tends to approach the flow state in a noisy environment. In DNS, the white noise term is added every 50 time steps for case A2 to allow its convection in the meshed domain (Cao et al. Reference Cao, Hao, Guo, Wen and Klioutchnikov2023). This interval corresponds to an angular frequency of 20ωp, where ωp is the angular frequency of the optimal planar wave. We have also simulated another DNS case forced only by the same noise term. It is found that no transition to turbulence occurs upstream of x = 0.5 m. Therefore, the considered white noise intensity should be insufficient to trigger transition in the unstable region of the eigenmodes. The white noise only serves to add initial randomness to the laminar flow region of case A2.
Table 1. Computational set-up for DNS cases.
The DNS cases A1 and A2 for the main discussion in the following contents impose a specific initial amplitude of the forcings. Without loss of generality, additional cases A3–A6 are simulated to examine the effects of the amplitude ratio and absolute amplitude of the forcings. As shown in table 1, cases A3–A5 are used to clarify the effect of the amplitude ratio by adjusting the modulus of the amplitude parameters 
${\epsilon_1}$ and 
${\epsilon_2}$ in (3.1), which are specified later. The responses of the mass-flow-rate fluctuations slightly downstream at x = 0.045 m are as follows: close but unequal in case A3, where 
${(\rho u)^{\prime}_{max,p}} = 2{(\rho u)^{\prime}_{max,o}} = 0.01$; planar wave dominates in case A4, where 
${(\rho u)^{\prime}_{max,p}} = 10{(\rho u)^{\prime}_{max,o}} = 0.01$; and oblique wave dominates in case A5, where 
${(\rho u)^{\prime}_{max,o}}\, = \,10{(\rho u)^{\prime}_{max,p}} = 0.01$. Another case with a weaker perturbation intensity is simulated to show the impact of the absolute amplitude, where 
${(\rho u)^{\prime}_{max,p}} = {(\rho u)^{\prime}_{max,o}} = 0.002$ at x = 0.045 m. Other settings remain identical to those of case A1.
For further details, first, let the fundamental angular frequency ω 0 and spanwise wavenumber β 0 be defined by ω 0Lref/u ∞ = 0.1 and β 0Lref = 0.8, respectively. In (3.1), the following resolvent analysis gives ωp = 10ω 0 and βp = 0 for the optimal planar wave, and ωo = 3ω 0 and βo = β 0 for the optimal oblique wave. Interestingly, ωp = 10ω 0, with dimensional frequency of f ≈ 125.8 kHz, is nearly the second-mode frequency corresponding to the maximum wall pressure fluctuation under the considered experiment condition by Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013). Thus, the optimal frequency should be representative. For the baseline case A1 in table 1, the complex amplitudes are specified by 
${\epsilon_1} = (1.554 \times {10^{ - 3}},\,2.1835 \times {10^{ - 7}})$ for the planar wave and 
${\epsilon_2} = ( - 7.6119 \times {10^{ - 4}}, - 6.4655 \times {10^{ - 5}})$ for the oblique wave. The modulus of the input forcing amplitude 
${\epsilon_1}$ or 
${\epsilon_2}$ determines the output amplitude of optimal disturbances linearly, while the phase of the complex amplitude results in the initial phase of output disturbances. These amplitudes are determined beforehand by resolvent analysis and PSE, which would result in the following modal properties. For the optimal planar wave, the initial amplitude forces the maximum dimensionless mass flow rate fluctuation at x = 0.045 m to be 
${(\rho u)^{\prime}_{max,p}} = 1 \times {10^{ - 2}}$ with a zero phase angle. In other words, the complex amplitude 
${\epsilon_1}$ can be obtained by a linear relation as 
${\epsilon _1} /1 \times {10^{ - 2}} = {\epsilon_{1,try}}/{(\widehat {\rho u})_{max,p,try}}$, where the subscript ‘try’ indicates a test run of resolvent analysis for linear amplitude scaling. Meanwhile, for the optimal oblique wave, the amplitude ratio to that of the planar wave is scaled based on the resolvent analysis with an equivalent 1:1 initial forcing norm (see § 2.2), which yields 
${(\rho u)^{\prime}_{max,o}} \approx 1.01 \times {10^{ - 2}}$ downstream at x = 0.045 m. As mentioned in § 1, it is reasonable to impose an initial amplitude on the first mode no smaller than that of the second mode. Furthermore, for the oblique wave propagated with a negative wave angle, the relationship between its forcing shape 
$\hat{f}_{-o}$ and 
$\hat{f}_{o}$ obeys the symmetry condition described in Park & Park (Reference Park and Park2013). The profiles of the forcing shape and phase angle are shown in figure 2, where the phase angle is in the range of (−180°, 180°]. Part of the non-zero forcing is located outside the boundary layer, which is consistent with the usual features of non-modal disturbances (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019). For the planar wave, as shown in figure 2(c), the energetic forcings are mostly above the relative sonic line for Mack second mode. This observation may indicate a difference in properties between non-modal and modal disturbances.
Figure 2. Profiles of the five components of the optimal forcing that are imposed at x = 40 mm for DNS: (a) modulus and (b) phase angle corresponding to the oblique wave, and (c) modulus and (d) phase angle corresponding to the planar wave. The plotted forcing 
$\hat{f}_{1}$ is non-dimensionalised by ρ ∞u ∞/Lref, 
$\hat{f}_{2}$, 
$\hat{f}_{3}$ and 
${\hat{f}_4}$ by 
${\rho _\infty }u_\infty ^2/{L_{ref}}$, and 
$\hat{f}_{5}$ by 
${\rho _\infty }u_\infty ^3/{L_{ref}}$. Dashed line in (a): location of the boundary layer thickness δ 0.99. Dashed-dotted line in (c): location of the relative sonic line of the second mode (angular frequency equals ωp), which satisfies U = cr − a. Here, U, a and cr represent the base streamwise velocity, base sound speed and phase velocity of the Mack second mode.
4. Results
4.1. Linearly optimal disturbances
Global resolvent analysis is performed over a wide range of spanwise wavenumbers and angular frequencies to determine the most amplified disturbances. The computational domain is extended from x 0 = 0.04 m to xe = 0.2 m, which is almost the pre-transitional region in the 3-D DNS. Figure 3 displays the contour of the optimal gain in the spectral domain. Generally, three types of optimal disturbances are identified in this free-stream condition, namely the optimal planar wave (10ω 0, 0β 0) appearing on the line β = 0, the optimal oblique wave (3ω 0, β 0) and the streak (10−5ω 0, 1.5β 0). Note that further lowering the frequency does not change the optimal gain of the streak. Among the three optimal disturbances, the planar and oblique waves possess the largest optimal gain, which below are called the (Fourier) modes (10, 0) and (3, 1), respectively. Figure 4 compares the modal shape of the streamwise velocity between the LPSE in figure 4(a,c,e) and the resolvent analysis in figure 4(b,d,f). In addition, the effect of adding the streak is compared in figures 4(g) and 4(h) by NPSE. The initial amplitudes of optimal disturbances are based on the resolvent analysis with a 1:1 forcing amplitude ratio, as stated in §§ 2 and 3. The agreement between LPSE and the resolvent analysis is excellent. The incident long-wavelength and upright short-wavelength patterns of the oblique and planar waves are observed. Based on figures 4(g) and 4(h), the streak does not essentially alter the combined disturbance shape owing to the weak gain in figure 3. Therefore, it is sensible to exclude the linearly optimal streak in 3-D DNS and consider the optimal oblique and planar waves only. Furthermore, it is observed that the coexistence of multiple optimal disturbances gives rise to an outer oblique-wave-dominant layer and an inner planar-wave-dominant layer.
Figure 3. Contours of optimal gain in the parameter space of the spanwise wavenumber and the angular frequency. Three plotted circles: local maxima of optimal gain corresponding to three types of optimal disturbances. Here, ω 0Lref/u ∞ = 0.1 and β 0Lref = 0.8.
Figure 4. Streamwise velocity contours of optimal disturbances u′, normalised by each maximum, by LPSE for (a) oblique wave, (c) planar wave and (e) streak, by resolvent analysis for (b) oblique wave, (d) planar wave and (f) streak, and by NPSE with the resolvent-based initial amplitude ratio for (g) oblique wave + planar wave and (h) oblique wave + planar wave + streak. Solid line: location of the base-flow boundary layer thickness δ 0.99. The maximum u′ for normalisation is 0.061 (a,b), 0.039 (c,d), 0.012 (e,f), 0.16 (g) and 0.058 (h).
Figure 5 shows a quantitative evaluation of the evolution of the optimal disturbances. The N-factor based on the Chu's energy is defined by
\begin{equation}{N_s} = 0.5\,\textrm{ln}({E_{Chu}}/{E_{Chu,0}}),\end{equation}
where EChu is defined by (2.11) and EChu ,0 denotes the value at x = 0.04 m. Thus, the N-factor characterises the exponential amplification of Chu's energy, and the streamwise derivative dNs/
${\rm d}\kern 0.06em x$ indicates the corresponding growth rate. A satisfactory agreement is achieved between LPSE and the resolvent analysis. In the meantime, LPSE initialised by LST profiles provides the evolution of the first and second modes denoted by open symbols. In figure 5(b), the vertical dashed lines mark the locations starting from which the difference in dNs/
${\rm d}\kern 0.06em x$ between the closed and open symbols is within 1%. The two locations, at x ≈ 0.136 m for the planar wave and x ≈ 0.164 m for the oblique wave, categorise the regions into upstream non-modal instability and downstream purely modal instability. Compared with the planar wave, the oblique wave possesses a stronger non-modal transient growth and weaker modal growth. This tendency is consistent with existing knowledge that the second mode is mostly dominant in the linear instability stage of hypersonic transitions. However, due to the prominent transient growth, the oblique wave has only a slightly lower N-factor than the planar wave at xe = 0.2 m. This result shows the possibility that the oblique first mode is as important as the second mode even in the linear instability stage. Without the transient growth effect, the open symbols illustrate that the N-factor of the planar wave exceeds that of the oblique wave by approximately 0.8 at xe, leading to a five-times-larger ratio of Chu's energy. Furthermore, EChu of the oblique wave is 11 times that of the streak at xe, indicating that the streak has a negligible impact.
Figure 5. Streamwise development of the (a) N-factor and (b) growth rate of Chu's energy dNs/
${\rm d}\kern 0.06em x$. Filled symbols, LPSE results initialised by resolvent-analysis profiles; open symbols, LPSE results initialised by LST profiles. Vertical dashed lines in (b): locations downstream of which the difference in dNs/
${\rm d}\kern 0.06em x$ between filled and open symbols is less than 1%, which can be viewed as the starting locations of purely modal instabilities for the planar and oblique waves.
4.2. The DNS results: instability waves and breakdown
Having examined the linear dynamics of optimal disturbances, we concentrate on the behaviour of the oblique and planar waves in the nonlinear regime. Firstly, the baseline case A1 and the case A2 with background noise are mainly discussed. Figure 6 shows the instantaneous and disturbed flow fields on the x–y plane. In addition to the outer oblique wave and inner planar wave, new streamwise flow scales characterising nonlinear sum or difference interactions appear upstream of xe = 0.2 m. The newly produced streamwise scales and downstream breakdown seem to occur earlier in the outer boundary layer than in the inner layer, as illustrated in figure 6(a). This observation may indicate the importance of oblique waves in the breakdown. These multiple streamwise scales are also displayed in figure 6(b) and a Rayleigh-scatter-like visualisation (Zhu et al. Reference Zhu, Zhang, Chen, Yuan, Wu, Chen, Lee and Elhak2016) in figure 6(e). Furthermore, the ‘rope-like’ structure of the density fluctuation, commonly observed for the second mode in experiments (Laurence, Wagner & Hannemann Reference Laurence, Wagner and Hannemann2016; Kennedy et al. Reference Kennedy, Laurence, Smith and Marineau2018b) and DNS (Pruett & Chang Reference Pruett and Chang1995; Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2020), is captured near the generalised inflection point in figure 6(d) and quantitatively shown in figure 6(f). Here, the generalised inflection point is the location where (∂/∂y)(ρb∂ub/∂y) = 0. In figure 6(f), a wavelength of around 5 mm is identified as the second-mode wavelength, since it is approximately twice the local boundary layer thickness δ 0.99 = 2.25 mm. However, another evident long wavelength of 18 mm is observed, which should not belong to the second mode. This differs from the previous conclusion in the literature that the ‘rope-like’ structure is attributed to the second mode (Laurence et al. Reference Laurence, Wagner and Hannemann2016). To determine the wavelength more accurately, narrow bandpass filtering is applied to the time series of ρ′, where the examined centre frequencies are ωo = 3ω 0 and ωp = 10ω 0. The corresponding streamwise wavelengths of the filtered data at x = 0.1 m are around 18.4 and 5.8 mm for the two frequencies, respectively. In § 4.3, figure 14(b) further demonstrates that the evident density fluctuation ρ′ near the rope-like structure is jointly contributed to by the oblique and planar waves. The different contributors to the rope-like structure in this paper and in the literature may arise from the differences in the geometry and flow environment. The present results here serve to show a possibility of the significance of low-frequency waves.
Figure 6. Flow fields on the plane z = Lz/2 for case A1 about (a) instantaneous u, full view, (b) instantaneous u, (c) disturbance u′, (d) disturbance ρ′ and (e) Rayleigh-scatter-like visualisation with a cut-off temperature T = 1.2, and (f) the disturbance ρ′ along the generalised inflection point in (d), with all non-dimensionalised by each free-stream base-flow value. Arrow in (a): the location where optimal forcings are imposed. Solid line in (d): location of the generalised inflection point. Notations in (b): the streamwise wavelengths of the oblique wave and of those generated by sum or difference nonlinear interactions.
The boundary layer also presents very different breakdown scenarios at various wall-normal heights, as shown in figure 7. The near-wall structures appear to be positive velocity streaks in figure 7(a). The inner layer consists of elongated structures with large velocity fluctuations first as an aligned pattern upstream and then as a staggered pattern downstream in figure 7(b). The outer layer mainly contains the staggered 
$\varLambda$-vortices with small velocity fluctuations in figure 7(c). Breakdown is found to have occurred at x > 0.35 m for all three slices. These observations may indicate the existence of fundamental, subharmonic and oblique breakdown in different streamwise or wall-normal locations. As a visualisation of vortices, figure 8 shows the Q-criterion iso-surface. The early spanwise distortion of the vortices and the formation of 
$\varLambda$-vortices and hairpin vortices are displayed.
Figure 7. Instantaneous dimensionless velocity fluctuation field u′ (compared with the base flow) on the x–z planes for case A1 at (a) y = 0.1 mm ( y/δ 0.99 ≈ 0.016 at x = 300 mm), (b) y = 1.3 mm ( y/δ 0.99 ≈ 0.21 at x = 300 mm) and (c) y = 3 mm ( y/δ 0.99 ≈ 0.49 at x = 300 mm).
Figure 8. The Q-criterion iso-surface (Lref/u ∞)2Q = 107 in the range 0.04 m < x < 0.3 m for case A1, coloured by the dimensionless pressure 
$p/{\rho _\infty }u_\infty ^2$.
To obtain a closer observation of the near-wall dynamics, figure 9 shows the instantaneous skin friction coefficient Cf for both cases A1 and A2. The streamwise travelling wave, oblique wave and large-scale 
$\varLambda$-vortex are seen sequentially. According to figure 9(b), the spanwise scale is mainly dominated by half the wavelength of the optimal oblique wave, i.e. λz ,o/2, which corresponds to the wavenumber 2β 0. This is mainly contributed by the streak (0, 2) characterising the oblique breakdown, as shown by the Fourier analysis next. In the streamwise direction, the planar-wave wavelength λx ,p is first dominant and the oblique-wave wavelength λx ,o is then highly pronounced. This can be understood by a subsequent modal analysis in which the planar wave saturates earlier than the oblique wave. What appears to be new is that staggered 
$\varLambda$-vortices are formed upstream of the eventual breakdown with a streamwise wavenumber that is not an integer multiple either of the oblique one or of the planar one. Hence, these 
$\varLambda$-vortices are inferred to be products of wave–wave interactions. Furthermore, the wavelength of the 
$\varLambda$-vortex is approximately five times the planar wavelength, corresponding to an angular frequency of 2ω 0 provided that phase locking occurs. The physical mechanisms behind this are focused on in the following sections. Finally, by comparing figures 9(b) and 9(c), the white noise effect is found only to affect local meticulous structures (e.g., in circled regions of figure 9(c)) and the group velocity of disturbances at approximately x > 0.22 m. This finding suggests that the background noise does not have an evident impact on the interactions of the optimal disturbances and their higher harmonics.
Figure 9. Instantaneous skin friction coefficient Cf of (a) full view for case A1, (b) zoomed-in view for case A1 (without the background noise effect) and (c) zoomed-in view for case A2 (with the background noise effect). The symbols λz,o and λx,o denote the spanwise and streamwise wavelengths of the oblique wave (3, 1), respectively, and λx,p and λx,(2,2) represent the streamwise wavelengths of the planar wave (10, 0) and mode (2, 2), respectively.
Statistical results are shown to identify the pre-transitional, transitional and fully developed turbulent regions. The results are ensured to be independent from the selection of temporal observation windows. Figure 10(a) shows that transition begins at x ≈ 0.13 m with a minimum Cf and ends at x ≈ 0.37 m with a Cf value collapsing onto the turbulent correlation. Moreover, a small deviation between cases A1 and A2 in Cf or the boundary layer thickness exists at approximately x > 0.26 m. This finding shows that the background noise effect is not very significant to the mean flow evolution. Figure 10(b) shows the streamwise velocity spectrum at x = 480 mm. The agreement with the law of ω −5/3 in the inertial subrange and the law of ω −7 in the dissipation scales suggests the establishment of a fully developed turbulent state at this station.
Figure 10. Quantitative results of (a) mean skin friction coefficient and boundary layer thickness for cases A1 (thick line) and A2 (thin line) and the van Driest II formula for Cf (dashed line), and (b) streamwise velocity spectrum 
${E_{11}} = \langle u^{\prime}u^{\prime}\rangle /u_\infty ^2$ at x = 480 mm for case A1, where 
$\langle \cdot \rangle$ denotes the ensemble average.
The bi-Fourier transform into the t–z space gives information on a series of Fourier modes (m, n). Convergence of spectra in the pre-transitional and transitional regions has been achieved based on the following settings: Hann window function is used with 50% overlapping, and the lowest angular frequency is 0.05ωp in each window. Firstly, for cases A1 and A2, figures 11(a) and 11(b) show the evolution of the maximum mass flow rate fluctuation and Chu's energy, respectively. The rapidly growing modes at x < 0.1 m include the oblique wave (3, 1) and the resulting streak (0, 2), the planar wave (10, 0), MFD (0, 0) and mode (10, 2) associated with the fundamental resonance. For these significant modes, the background noise only has a limited effect in the late stage of nonlinear instability and in the turbulent region. It is concluded that the background noise can neither select new significant modes with different frequencies or wavenumbers nor manipulate the main characteristics of the existing modes. In addition, modes (2, 0) and (2, 2) corresponding to the large-scale 
$\varLambda$-vortices in figure 9 and mode (5, 1) associated with the subharmonic resonance with the optimal planar wave start to be amplified rapidly in the vicinity of x = 0.1 m. In terms of the 
$\varLambda$-vortices, some neighbouring modes such as (1, 0) and (1.5, 0), with approaching streamwise scales to modes (2, 0) and (2, 2), are shown not to be pronounced in figure 11. Therefore, there should be a frequency selection process, which is discussed next. For brevity, evolution of some other significant modes is shown in Appendix B. The main interaction occurring in the oblique breakdown related to mode (3, 1), the fundamental resonance related to mode (10, 0) and the subharmonic resonance related to mode (10, 0) can be expressed by
\begin{gather}(3,1) - (3, - 1) \to (0,2),\end{gather}
\begin{gather}(10,\textrm{0}) + (0,\textrm{2}) \to (10,2),\end{gather}
\begin{gather}(10,\textrm{0}) - (5,\textrm{1}) \to (5, - \textrm{1}),\end{gather}respectively. The corresponding reverse interactions among the triads are also included, such as (0, 2) + (3, −1) → (3, 1) for the oblique breakdown.
Figure 11. Streamwise development of disturbances on (a) maximum mass flow rate 
${(\rho u)^{\prime}_{max}}$ for cases A1 (solid line) and A2 (dashed line), (b) Chu's energy EChu for cases A1 (solid line) and A2 (dashed line), (c) Chu's energy EChu for case A4 (initial forcing ratio of planar to oblique waves equals 10:1), (d) Chu's energy EChu for case A5 (initial forcing ratio of planar to oblique waves equals 1:10) and (e) contour for logarithmic fast Fourier transform of wall pressure fluctuation for case A1, ln(|p′|).
The successive appearance of highly amplified modes (3, 1), (0, 2), (10, 0), (10, 2) and detuned (2, 2) or (2, 0) in figures 11(a) and 11(b) indicates the occurrence of oblique and fundamental breakdowns in conjunction with an unconfirmed combination resonance. With regard to the modal amplitude, the streak mode (0, 2) is dominant in the early region x < 0.2 m. Comparing the downstream first and Mack second modes, those of the oblique wave (3, 1) are stronger than and saturate later than those of the planar wave (10, 0). These results are consistent with the observations of flow structures in figures 6–9. However, the wall pressure fluctuation shows a considerable magnitude for the second-mode angular frequency 10ω 0 at x < 0.2 m, as shown in figure 11(e). It can thus be inferred that the wall pressure measurement in experiments may not completely represent the possible modal competition in the outer layer of the boundary layer away from the wall, particularly for non-acoustic disturbances. The integral EChu, including the region away from the wall and containing kinetic and thermodynamic energy simultaneously, may characterise the modal growth and structure instabilities in figures 6–9 more comprehensively. In addition to the sum angular frequency 13ω 0 and difference angular frequency 7ω 0 in figure 11(e), the interaction between (3, 1) and (10, 0) generates a series of angular frequencies that are only integer multiples of the fundamental one ω 0. These are caused by the fact that the two identified optimal frequencies 10ω 0 and 3ω 0 under this freestream condition are coprime. Physically, the spectral broadening during the transition could be enhanced by these discrete modes.
In terms of a comparison with other well-known breakdown scenarios, further information is shown by figure 11(c) for the planar-wave-dominant case A4 and by figure 11(d) for the oblique-wave-dominant case A5. Several interesting aspects are discussed as follows.
(i) For case A4 where the initial forcing ratio of planar to oblique waves equals 10:1, the planar wave (10, 0) saturates early in the vicinity of x = 0.2 m before the MFD is strong enough. Meanwhile, modes (10, n) are found not to be amplified enough, including those unshown with n ≠ 2. This indicates that fundamental resonance does not trigger the eventual breakdown due to a relatively weak second-mode wave in the considered flat-plate boundary layer, which may differ from those scenarios for a cone by Fasel's group. Particularly, the flared-cone flow adding the background noise by Hader & Fasel (Reference Hader and Fasel2018) supported a constantly growing second mode with a slowly varying boundary layer thickness, which resulted in the spectral peak of the second mode at 300 kHz and its harmonics. Under our conditions, the optimal planar-wave frequency f ≈ 125.8 kHz, which is nearly the frequency of peak pressure fluctuations in the experiment by Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013), does not possess a long unstable streamwise region. By contrast, the oblique wave of case A4 undergoes a much longer distance of the instability stage in this state. The amplified oblique wave exceeds the planar wave finally even with a lower initial amplitude. The downstream scenario in case A4 is more likely to be a joint type of oblique breakdown, combination resonance and second-mode subharmonic breakdown. The significance of low-frequency oblique waves in the late stage is somewhat consistent with the conclusions under the flow and geometric conditions at Peking University (Zhu et al. Reference Zhu, Zhang, Chen, Yuan, Wu, Chen, Lee and Elhak2016; Chen et al. Reference Chen, Zhu and Lee2017).
(ii) The similarities between cases A1 and A5 upstream of x = 0.2 m suggest that the oblique breakdown is quite important in the early stage. Further downstream, with the existence of mode (10, 0), the possible detuned mode (2, 2) can be promoted. Among the comparative cases, case A1 reports the most prominent amplification of mode (2, 2). With regard to the signature of oblique breakdown, the streak (0, 2) of case A1 becomes weaker than that of case A5 in the saturation stage downstream of x = 0.2 m. Consequently, oblique breakdown plays a crucial role mainly in the early stage.
(iii) By comparing case A1 with cases A4 and A5, case A1 shows a slightly lower Chu's energy for both modes (10, 0) and (3, 1) upstream of x = 0.2 m. With equal initial forcings of the primary oblique and planar waves in case A1, more energy is probably transferred to accelerate the growth of secondary modes, such as the detuned one (2, 2) and subharmonic one (5, 1).
The modal growth results in figure 11 may not show straightforwardly that mode (2, 2) or (2, 0) is important in the late nonlinear stage. To provide more direct evidence, the maximum absolute contribution from mode (m, n) to the instantaneous skin friction is defined by
\begin{equation}\Delta {C_{f,(m,n)}}(x) = \mathop {\max }\limits_{(z,t)} |{{C_{f\textrm{,(}m,n\textrm{)},disturbed}} - {C_{f,laminar}}} |,\end{equation}
where Cf ,(m,n),disturbed and Cf ,laminar are the instantaneous skin friction induced by the laminar flow + mode (m, ±n) alone and the laminar flow alone, respectively. This indicator is defined to highlight and visualise the Fourier modes which are the source contributors to Cf in figure 9. Figure 12(a) demonstrates that mode (2, 2) is dominant in the vicinity of x = 0.3 m for both cases A1 and A2, excluding the MFD and mode (0, 2) which do not directly affect the streamwise scales. Since the large-scale 
$\varLambda$-vortices are observed in the same region in figure 9, it is clarified that mode (2, 2) is responsible for the formation of near-wall 
$\varLambda$-like structures for case A1. In addition to the baseline forcing settings, the effects of the amplitude ratio and the absolute amplitude of the initial forcings are of further interest. Figure 12(b–d) shows the results when the initial forcing ratio of the planar wave to the oblique one is set to 2:1 (case A3), 10:1 (case A4) and 1:10 (case A5), respectively. Figure 12(e) displays the result of case A6 when the absolute initial forcing decreases to 0.2 times the baseline value for both planar and oblique waves. It is observed that the pronounced region of mode (2, 2) is postponed compared with the baseline case, and that the corresponding proportion of the contribution to the skin friction becomes smaller. However, the presence of this detuned mode is found to be independent of the effects of the amplitude ratio and absolute amplitude under the considered parametric set-up.
Figure 12. Maximum absolute modal contribution to the instantaneous skin friction coefficient ΔCf for (a) cases A1 (solid line) and A2 (dashed line), (b) case A3, (c) case A4, (d) case A5 and (e) case A6.
Figure 13 further shows a snapshot of the skin friction coefficient for cases A3–A6. Note that in cases A4 and A5, the initial planar and oblique waves are dominant, respectively. Consequently, early fundamental and oblique breakdown scenarios characterised by the aligned and staggered vortices are shown in figures 13(b) and 13(c), respectively. In terms of the combination resonance, the signature of the large-scale 
$\varLambda$-vortices becomes weaker if the amplitude ratio or the absolute amplitude changes. As shown in figure 13(a–d), the spanwise length scale of the 
$\varLambda$-vortices may be visibly affected because of other mixed-up modes, and accordingly the head angle of the vortex is altered. However, the signature of the large streamwise length scale corresponding to the detuned frequency 2ω 0 is still maintained in the late stage of the transition, although it is postponed depending on the active region of the detuned mode (2, 2) in figure 12. These findings possibly indicate the general existence of the combination resonance mechanism if multiple instability waves are seeded. The detailed strength of this resonant scenario depends on the flow conditions, the forcing components, amplitudes, etc. The features of strength dependence and existence independence of the combination resonance have also been reported by Fezer & Kloker (Reference Fezer and Kloker2000) in supersonic and hypersonic states. In their studies, the included inflow modes (1, ±4) and (1/2, ±3) with different (relatively) low amplitudes at Mach numbers 2 and 6.8 could generate a rapidly growing ‘combination mode’ (1/2, ±7). Differently, the combination resonance under the condition of Fezer & Kloker seemed not to be pronounced enough to exceed the subharmonic and oblique breakdown scenarios somewhere and present evident flow phenomena. In addition, this mode (1/2, ±7) did not introduce new streamwise scales compared with the inflow modes. The emergence of mode (1/2, ±7) is also less complicated than that of mode (2, 2) in this paper (see (4.6)).
Figure 13. Instantaneous skin friction coefficient Cf of (a) case A3, (b) case A4, (c) case A5 and (d) case A6. The symbols λz,o and λx,o denote the spanwise and streamwise wavelengths of the oblique wave (3, 1), respectively, and λx,p and λx,(2,2) represent the streamwise wavelengths of the planar wave (10, 0) and mode (2, 2), respectively. The contour levels are identical to those in figure 9.
4.3. Weakly nonlinear stability analysis
In this section, the NPSE serves to reveal part of the physical mechanisms to explain the flow phenomena in § 4.2. As a simplified physical model for the weakly nonlinear stability analysis, only sum and difference modes of the seeded inflow disturbances are additionally generated in the NPSE simulation. Figure 14 compares the evolution of Chu's energy and profiles of the density fluctuation between the DNS and NPSE before the NPSE fails to converge. Generally, the agreement is fine. Meanwhile, the approach of the density peak between mode (3, 1) and (10, 0) supports that the evident density fluctuation of the ‘rope-like’ structure in figure 6(f) is jointly contributed to by the oblique and planar waves.
Figure 14. Comparison between DNS (solid line) and NPSE (symbol) of (a) Chu's energy and (b) profiles of density fluctuation |ρ′| at x = 0.12 m.
In addition to the self-interactions, mutual sum or difference interactions between the planar and oblique waves enable the occurrence of numerous modes. Figure 15 shows the weakly nonlinear effect on either the planar or oblique wave, where the sum or difference interaction is artificially cut off in the NPSE code. In figure 15, case B0 (see table 2) is the baseline case which is compared with DNS data. Case B0 considers the mutual interactions and self-interactions between modes (10, 0) and (3, ±1) (as well as MFD, not mentioned below). Compared with the baseline B0, case B1 interrupts the difference interaction between planar and oblique waves, which generates (7, ±1), while case B2 removes their sum interaction, which yields (13, ±1). Obviously, by comparing cases B1 and B2 with the baseline case, the difference interaction between the oblique and planar waves destabilises both oblique and planar waves, while the sum interaction stabilises both waves. The slight shift in EChu indicates that the direct interaction between the oblique and planar waves is weak in the initial stage. However, the two generate a series of Fourier modes via multiple generations of nonlinear interactions, which are important in the eventual breakdown. Among them, the possible detuned mode (2, 2) is of particular interest.
Figure 15. Effects of the absence of sum or difference interactions for planar waves (thick lines) and oblique waves (thin lines) on EChu.
Table 2. Settings for NPSE cases on the effects of sum or difference interactions.
The selection mechanism of the frequency and spanwise wavenumber for mode (2, 2) can be straightforwardly understood from the procedure of NPSE simulations. Corresponding to the dominant wavenumber 2β 0, no other preferential frequency in the neighbourhood except for 2ω 0 is found to be largely amplified. The detailed generation process of mode (2, 2) in the NPSE is described by
\begin{gather}(3,1) + (3, - 1) \to (6,0),\end{gather}
\begin{gather}(10,0) - (3, - 1) \to (7,1),\end{gather}
\begin{gather}(7,1) - (6,0) \to (1,1),\end{gather}
\begin{gather}(1,1) + (1,1) \to (2,2).\end{gather}Aimed at revealing the physical mechanism of mode (2, 2), we design several NPSE cases with information given in table 3. The modes are seeded at x = 0.045 m in the NPSE with DNS-initialised profiles. Note that the primary waves have evidently larger amplitudes than the secondary instability modes in the early stage. With these primary instability waves added into the background flow, secondary instability is likely to occur. To simplify the physical problem, only primary waves (10, 0) and (3, ±1) and the concerned mode (2, 2) are included in the inflow profiles. The purpose is to identify the necessary primary instability waves which support the rapid growth of mode (2, 2).
Table 3. Settings for NPSE cases on the generation mechanism of mode (2, 2).
In table 3, cases C0 and C1 simulate the development of the small-amplitude mode (2, 2) under the large-amplitude planar and oblique waves and the base flow. Different from case C0, case C1 interrupts the sum and difference interactions between (10, 0) and (3, ±1) downstream of the inflow station x = 0.045 m. Accordingly, case C1 resembles a Floquet analysis on the secondary instability (Ng & Erlebacher Reference Ng and Erlebacher1992), where the original base flow, mode (10, 0) and mode (3, ±1) constitute a new periodic background flow. In the frame of reference moving with the group velocity of the primary waves, the combination of the primary mode (10, 0) and mode (3, ±1) varies slowly as a background flow, while small-amplitude higher-order harmonics can undergo secondary super-exponential growths. In comparison, case C2 or C3 only contains mode (10, 0) or (3, ±1) as the background primary wave, respectively. To approach the Floquet theory, a phase-locked assumption can be made. Here, we follow the usage of the terminology ‘phase-locked’ by Chang et al. (Reference Chang, Malik, Erlebacher and Hussaini1993) in NPSE studies. Actually, ‘phase-locked’ usually refers to a state whereby two types of waves have nearly the same phase velocity (Chen et al. Reference Chen, Zhu and Lee2017), i.e. a lock of the phase velocity. Once the phase-locked state occurs, the generated sum or difference mode can constitute a resonant triad with the original two waves, which can give rise to a super-exponential growth. Specifically, the iterative scheme of the streamwise wavenumber except for mode (10, 0) is changed from (2.17) to
\begin{gather}{\mathop{\rm Im}\nolimits} ({\alpha _{mn,new}}) \!=\! {\mathop{\rm Im}\nolimits} ({\alpha _{mn,old}}) - \textrm{iRe}\left[ {\frac{1}{{{E_{mn}}}}\int_0^\infty {\bar{\rho }\left( {\hat{u}_{mn}^ \ast \frac{{\partial {{\hat{u}}_{mn}}}}{{\partial x}} \!+\! \hat{v}_{mn}^ \ast \frac{{\partial {{\hat{v}}_{mn}}}}{{\partial x}} + \hat{w}_{mn}^ \ast \frac{{\partial {{\hat{w}}_{mn}}}}{{\partial x}}} \right)\textrm{d}y} } \right],\end{gather}
\begin{gather}\textrm{Re}({\alpha _{mn,new}}) = \textrm{Re}({\alpha _{(10,\textrm{0})}}) \times m/10,\end{gather}while mode (10, 0) maintains the calculation manner in (2.17). Provided that wave–wave resonance occurs, the synchronisation of the phase velocity allows the real streamwise wavenumber to be proportional to the angular frequency, as described by (4.7b).
In this complicated problem, multiple Fourier modes with approaching strength are able to affect the shape of mode (2, 2) at different wall-normal heights. Therefore, the growth of Chu's energy rather than the mode shape is chosen as the indicator to identify the dominant effect in the underlying resonant state. Figure 16 illustrates that both cases C0 and C1 can achieve an exponential growth for mode (2, 2) with a slope close to that of the DNS counterpart. The initial difference between the NPSE and DNS may come from the simplification manipulation of the NPSE in the governing equations and interaction models. Nevertheless, the agreement in the slope in the exponential amplification stage confirms that the secondary instability mechanism is mainly responsible for the growth of mode (2, 2). Furthermore, the difference between cases C0 and C1 suggests that the interaction between modes (10, 0) and (3, ±1) has a lesser effect on the amplification of mode (2, 2). The interaction between modes (10, 0) and (3, ±1) generates additional modes and constantly feeds energy to mode (2, 2) through the indirect route via these new modes in addition to a direct interaction. By comparing cases C0 or C1 with C2 and C3, it is clear that the planar or oblique primary wave alone cannot induce the secondary instability of mode (2, 2) effectively. Hence, the combination resonance of the detuned mode (2, 2) requires participation of both mode (10, 0) and mode (3, ±1). This resonant mechanism is different from conventional mechanisms, where only one primary wave is needed. We have also ruled out the rapid growth of mode (2, 2) caused only by interactions between mode (3, ±1) and other 3-D modes by additional unshown cases. In these cases, the seeded inflow modes include (3, ±1), (2, ±2) and one of the following modes to constitute an interaction triad: (1, ±1), (1, ±3), (5, ±1) or (5, ±3). The growth of mode (2, 2) in these unshown cases resembles that in cases C2 and C3. Thus, the inclusion of the primary mode (10, 0) is required to support the significant amplification of mode (2, 2). Figure 16 also shows that the phase-locked assumption essentially does not alter the results, which agrees with the conclusion by Chang et al. (Reference Chang, Malik, Erlebacher and Hussaini1993).
Figure 16. Comparison between the NPSE results of EChu for cases C0 to C3 and the DNS results. Dashed line: results without the phase-locked assumption.
4.4. Resonant state and energy transfer mechanisms
The resonant state can be examined by plots of wave vectors of the Fourier modes (Fezer & Kloker Reference Fezer and Kloker2000). If the three wave vectors (αr, β) of the considered triad constitute a closed form, the resonance condition will be satisfied (Craik Reference Craik1971):
\begin{equation}{\alpha _{r,1}} + {\alpha _{r,2}} = {\alpha _{r,3}},\quad {\beta _1} + {\beta _2} = {\beta _3},\quad {\omega _1} + {\omega _2} = {\omega _3},\end{equation}
where the angular frequency relation is automatically satisfied for the examined Fourier triad. The streamwise wavenumber is obtained via αr = d
$\varTheta$/
${\rm d}\kern 0.06em x$, where 
$\varTheta$ is the phase angle of Fourier transform of the wall pressure in DNS. Figures 17(a) and 17(b) depict the wave vectors of the four interactions in (4.6) at x = 0.2 m and x = 0.3 m, respectively. The resonance condition for the generation of mode (2, 2) is generally fulfilled at x = 0.2 m, although the self-interaction (1, 1) + (1, 1) → (2, 2) does not fully form a closed triad. Downstream at x = 0.3 m, better closed triads indicate an enhanced resonant state. The resonance of the planar wave, the oblique wave and the detuned mode is also shown by the well-closed vector chain (2, 2) + (3, 1) + (3, −1) + (2, −2) → (10, 0) in figure 17(c).
Figure 17. Wave vectors (αr, β) of the modes in (4.6a–d) at (a) x = 0.2 m and (b) x = 0.3 m, and (c) the vector chain (2, 2) + (3, 1) + (3, −1) + (2, −2) → (10, 0) at x = 0.3 m for case A2.
To understand the energy transfer in the nonlinear stage, a posterior energy budget analysis is presented first. With regard to mode (2, 2), the energy budget of the NPSE cases C0 to C3 in § 4.3 is considered here. Similar to the derivation by Chen et al. (Reference Chen, Zhu and Lee2017), the wall-normal integral of the fluctuation kinetic energy equation for mode (m, n) can be written as
\begin{equation}\int_0^\infty {\frac{1}{2}\bar{\rho }\frac{{D({{|{{u_{mn}}} |}^2} + {{|{{v_{mn}}} |}^2} + {{|{{w_{mn}}} |}^2})}}{{Dt}}\,\textrm{d}y} = \mathcal{P} + \mathcal{Q} + \mathcal{F} + \mathcal{D} + \mathcal{N},\end{equation}
where the right-hand side contains the terms in the form of a wall-normal integral, including the linear mean-shear production term 
$\mathcal{P}$, the non-parallel term 
$\mathcal{Q}$, the pressure diffusion and dilatation term 
$\mathcal{F}$, the viscous dissipation term 
$\mathcal{D}$ and the nonlinear term 
$\mathcal{N}$. Detailed expressions can be found in Chen et al. (Reference Chen, Zhu and Lee2017). In addition, 
$\mathcal{L}$ represents the sum of all linear terms on the right-hand side.
Figure 18 indicates that the rapid energy growth of mode (2, 2) in cases C0 and C1 is caused by the positive nonlinear term. In contrast, the sums of the linear terms as well as most linear terms alone are negative and stabilise mode (2, 2), although the mean shear production effect is slight destabilisation. In terms of cases C2 and C3, the term magnitude is very small and mode (2, 2) fails to become strong.
Figure 18. Energy budget of mode (2, 2) of (a) cases C0 (solid line) and C1 (dashed line) and (b) cases C2 (solid line) and C3 (dashed line). All the terms of the four cases are normalised by the maximum 
$\mathcal{N}$ of case C0.
In the present study, it is difficult for the energy budget analysis to provide a comprehensive interpretation of energy transfer between multiple discrete Fourier modes. We further apply the ‘amplitude correlation method’ first proposed by Crossley et al. (Reference Crossley, Uddholm, Duncan, Khalid and Rusbridge1992) and Duncan & Rusbridge (Reference Duncan and Rusbridge1993) in plasma dynamics. Xia & Shats (Reference Xia and Shats2003) used this technique to provide the first experimental evidence of an inverse energy cascade in turbulence with broadband spectra. Zhang & Shi (Reference Zhang and Shi2022) employed this method to analyse the modal interaction in the hypersonic boundary layer transition. The realisation is based on the cross-correlation function (CCF) in the form of the fourth-order covariance as
\begin{equation}K(\tau ) = \langle [x_I^2(t)][x_{II}^2(t + \tau )]/({\parallel} [x_I^2(t)]\!\parallel{\!\parallel} [x_{II}^2(t + \tau )]\!\parallel )\rangle .\end{equation}
Here, xI and xII are filtered signals of the same original sample, which are centralised on the frequencies 
$f_{1}$ and f2, respectively. The corresponding filters F1 and F2 are chosen as the Butterworth bandpass filter with the flattest frequency response and no phase shift. Moreover, [
${\cdot}$] denotes a low-frequency filtration by filter F3 first and an operation to subtract the expected value then to obtain the fluctuation. The symbol ||
${\cdot}$|| denotes the Euclidean norm for normalisation, and 
$\langle \cdot\rangle $ means the ensemble average. Finally, the CCF K(τ) is a function of the time lag τ. If the 
$\tau_{lag}$ that maximises the CCF is positive, the signal xII lags xI and energy flows from component f1 to f2 through the sideband nonlinear interaction. Otherwise, energy is considered to flow from f2 to f1 instead.
The physical explanation of this technique is briefly introduced here. Assume that the narrowband signals xI and xII contain major energy in the frequency ranges [f1 − Δf1, f1 + Δf1] and [f2 − Δf2, f2 + Δf2], respectively, similar to those in figure 11(e). Thus 
$x_I^2$ has two bands [2f1 − 2Δf1, 2f1 + 2Δf1] and [0, 2Δf1] (similarly for 
$x_{II}^2$). The low-frequency filter F3 removes the higher band and maintains the lower band. Thus the CCF characterises information transfer between frequencies f1 and f2 through the sidebands with lengths of 2Δf1 and 2Δf2, and thus contains information on the direction of energy transfer. In the present study, the width of filter F3 is varied from 0.1f0 to f0 to ensure the independence of Sgn(
$\tau_{lag}$), where Sgn(
${\cdot}$) is the sign function.
Figure 19 shows the signal samples of the wall pressure fluctuation 
$x_I^2$ (f1 = 2f0) and 
$x_{II}^2$ (f2 = 19f0) and the filtered ones by F3 as well as the CCF. The sampled station is at x = 300 mm, where the strong nonlinear interaction gives rise to stochastic time signals and the concerned 
$\varLambda$-vortices are evident. The pressure signal is chosen for the analysis to approach an experimental signal (Zhang & Shi Reference Zhang and Shi2022). The two example results of the CCFs in figure 19(c) indicate that the frequency f1 = 2f0 lags both frequencies f2 = 19f0 and f2 = 10f0, and thus energy flows from these two high frequencies to f1 = 2f0.
Figure 19. Signal sample of wall pressure fluctuations at x = 300 mm for case A2 of (a) 
$x_I^2$ for f = 2f0 and (b) 
$x_{II}^2$ for f = 19f0, and (c) the resulting CCF. Thick lines in (a,b): filtered signal by the low-frequency filter F3. Dashed-dotted line in (c): another CCF plot with 
$x_{II}^2$ belonging to f = 10f0.
By means of the approach above, how mode (2, 2) sustains its energy and supports the formation of 
$\varLambda$-vortices can be analysed. Let f1 = 2f0, and the characteristic time lag 
$\tau_{lag}$ is obtained via the calculation of CCF for each pair of (f1, f2). If Sgn(
$\tau_{lag}$) > 0, energy flows from f1 to f2; if Sgn(
$\tau_{lag}$) < 0, energy flows from f2 to f1; if Sgn(
$\tau_{lag}$) = 0, there is no evident peak of the CCF with a non-negligible cross-correlation, say, (CCF)max < 0.1. Figure 20 shows a plot of Sgn(
$\tau_{lag}$) and the resulting diagram of the energy flow direction. It is known from the energy budget analysis in figure 18 that mean shear slightly destabilises mode (2, 2) and other linear terms stabilise it. Therefore, mean shear feeds energy to mode (2, 2), and mode (2, 2) dissipates energy into other linear mechanisms. Meanwhile, the constantly growing MFD is the product of the sum interaction between each mode and its complex conjugate, i.e. (2, 2) + (−2, −2) → (0, 0). It is thus inferred that mode (2, 2) participates in energy transfer to the MFD. The remaining groups in figure 20(b) are categorised into nonlinear energy ‘contributors’ and ‘receivers’. In the studied case, the energy contributors for f1 = 2f0 consist of two neighbours f2 = f0 and f2 = 4f0, the optimal planar wave frequency f2 = 10f0 and additional higher modes from f2 = 12f0 to f2 = 20f0. Note that the low frequency f2 = f0 is also one of the source modes for mode (2, 2) via (1, 1) + (1, 1) → (2, 2) (see (4.6) for details). The energy receivers for mode (2, 2) are the optimal oblique wave frequency f2 = 3f0 and other components f2 = 8f0, f2 = 9f0 and f2 = 11f0 in the vicinity of the frequency of mode (10, 0). Corresponding to each energy flow direction from f2 to f1 or from f1 to f2, there are two types of potential interaction processes to allow nonlinear energy transfer: the sum interaction (f1 + f2 → f3 = f1 + f2) and difference interaction (f1 − f2 → f3 = |f1 − f2|) in conjunction with their reverse counterparts, such as (f1 + f2) − f1 → f2. In summary, figure 20(b) gives a physical image of energy transfer through linear mechanisms and through the narrow sidebands in nonlinear modal interactions.
Figure 20. The energy flow direction at x = 0.3 m for f1 = 2f0 of case A2: (a) sign function of the time lag of f2 compared to f1 and (b) diagram of the energy flow direction.
Note that the energy flow direction may depend on the signal type and wall-normal height for analysis. Amplitude correlation analysis is further performed based on the (ρu)′ signal at the height of its maximum for mode (3, 1) at x = 300 mm. What is different is that f1 = f0 and the high frequencies f1 = 13f0, 14f0 and 15f0 change from ‘contributors’ to ‘receivers’, while f1 = 3f0 changes from ‘receivers’ to ‘contributors’. However, the behaviour that multiple frequencies participate in nonlinear energy transfer as ‘contributors’ and ‘receivers’ is found to be consistent.
5. Conclusions
Concentrating on the linear and nonlinear dynamics of optimal disturbances, a Mach 6 flat-plate boundary layer is studied by resolvent analysis, PSE and DNS. Newly observed flow phenomena and physical interpretations of the modal interaction and breakdown scenario are presented.
The resolvent analysis demonstrates that the optimal planar wave mode (10, 0) and the optimal oblique wave mode (3, 1) are notably more significant than the optimal streak. Due to the strong transient growth, the integrated Chu's energy of the oblique wave is comparable with that of the planar wave, although the growth rate of the second mode exceeds that of the first mode. Introducing the two optimal disturbances jointly with the same magnitude of initial forcing, the DNS results show that the oblique wave is highly pronounced in the outer boundary layer. Quantitatively, the active performance of the non-acoustic oblique wave away from the wall may not be fully detected by conventional wall pressure measurements, such as PCB sensors.
Modal decomposition via a Fourier transform identifies the coexistence of the oblique breakdown, fundamental resonance triggered by mode (10, 0), combination resonance related to a detuned mode (2, 2) and subharmonic resonance triggered by mode (10, 0). Particularly, mode (2, 2) is responsible for the occurrence of large-scale staggered 
$\varLambda$-vortices near the wall before the eventual breakdown. More importantly, the presence of mode (2, 2) and the 
$\varLambda$-vortices are found not to depend on the (considered) amplitude ratio and absolute amplitude of the initial forcings corresponding to the oblique and planar waves. This suggests that the combination resonance can be a viable transition scenario with the presence of multiple instability waves. The NPSE analyses further confirm that the detuned mode (2, 2) is mainly a consequence of the secondary instability, where the background primary waves require participation of both modes (10, 0) and (3, 1). The nonlinear interaction between the two types of primary waves generates a series of discrete Fourier modes. These multiple components of spectra in the early nonlinear stage are different from the observation in a single common breakdown scenario. In addition, the resonant state associated with these multiple components is observed. Finally, the energy budget and amplitude correction techniques describe the energy flow direction. The linear production and dissipation effects as well as nonlinear interactions with ‘contributor’ modes and ‘receiver’ modes are involved in the physical image of energy transfer for mode (2, 2). The present study indicates that the boundary layer transition induced by multiple primary instability waves can contain rich flow physics and unexpected breakdown scenarios. Particularly, the significance of combination resonance and primary oblique wave in a hypersonic state is highlighted in this work. The multiple resonance mechanisms including the combination one, revealed by a specific initialisation of frequency and wavenumber pair in the present paper, may constitute a building block in a real natural transition with a wide and continuous frequency range.
Funding
This work is supported by the Hong Kong Research Grants Council (no. 15206519, no. 25203721 and no. 15216621) and the National Natural Science Foundation of China (no. 12102377).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of the NPSE solver
The NPSE part of the present code CHASES is validated through comparisons with an oblique breakdown NPSE case (Hein Reference Hein2005) and an NPSE case containing multiple breakdown scenarios including fundamental and subharmonic ones (Chang & Malik Reference Chang and Malik1994). The base flow of both cases is the flat-plate boundary layer at Mach 1.6 with an adiabatic wall. A pair of oblique waves (1, ±1) are added in the inflow of the case by Hein (Reference Hein2005). Both oblique waves (1, ±1) and (2, ±1) and the planar wave (2, 0) are included in the inflow condition by Chang & Malik (Reference Chang and Malik1994). Good agreement in the modal amplitude is shown in figure 21, which demonstrates the reliability of the present NPSE code.
Figure 21. Comparison of the NPSE results of the maximum streamwise velocity fluctuation between CHASES (solid line) and the literature (symbol) of (a) Hein (Reference Hein2005) and (b) Chang & Malik (Reference Chang and Malik1994). The Reynolds number Re is based on the length scale l = (
$\mu_{\infty}$x/ρ∞u∞)1/2.
Appendix B. Evolution of some other Fourier modes
There exist some other significant Fourier modes in scenarios such as the oblique breakdown, which are not shown in figure 11. For a comprehensive examination, evolutions of modes (1, 3), (3, 3), (6, 0) and (7, 1) of case A1 are depicted in figure 22, which are usually pronounced in the pure oblique breakdown. Clearly, these modes can be amplified rapidly, which indicates the significance of the conventional oblique breakdown and other potential mechanisms in the initial stage. However, in the vicinity of x = 0.3 m where the 
$\varLambda$-vortex appears, the detuned mode (2, 2) has an evident advantage in the energy magnitude over most of the remaining modes. Thus, the significant role of the combination resonance before the eventual breakdown could be further justified.
Figure 22. Streamwise development of Chu's energy of other significant Fourier modes for case A1.
Open access
