1. Introduction
Understanding and modelling high-speed turbulent flows are of fundamental importance in aerodynamic applications (Spina, Smits & Robinson Reference Spina, Smits and Robinson1994). The turbulent boundary layer influences surface drag and heat transfer dramatically. Thereby, developing accurate prediction models has long been the pursuit for reliable vehicle design. Compared with the incompressible case, our knowledge of compressible turbulent flows is rather limited, for two main reasons. First, the compressible flow suffers from coupling effects of various factors, which, in addition to the Reynolds number, include the Mach number, surface heat transfer, shock waves and high-enthalpy effects (Coleman, Kim & Moser Reference Coleman, Kim and Moser1995; Duan & Martín Reference Duan and Martín2011; Fu et al. Reference Fu, Karp, Bose, Moin and Urzay2021; Fu, Bose & Moin Reference Fu, Bose and Moin2022). Second, the available experimental and numerical databases are rather limited because of the restricted facilities and computational resources.
With the advancement of modern supercomputers, direct numerical simulations (DNS) (see e.g. Moin & Mahesh Reference Moin and Mahesh1998) are increasingly powerful for the research of compressible turbulent flows. For example, recent DNS data from Zhang, Duan & Choudhari (Reference Zhang, Duan and Choudhari2018) and Huang, Duan & Choudhari (Reference Huang, Duan and Choudhari2022) cover a freestream Mach number 
${{Ma}}_\infty$ range of up to 14 for flat-plate boundary layers with adiabatic and cold walls. In addition, the maximum friction Reynolds number 
${{Re}}_\tau$ reaches as high as 4000 for the boundary layer (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2013, with 
${{Ma}}_\infty =2$) and 2000 for the channel flow (Yao & Hussain Reference Yao and Hussain2020, with a bulk Mach number 
${{Ma}}_b=1.5$). Utilizing these DNS data, the effects of Mach number, wall cooling, pressure gradient, high-enthalpy, and so on, on turbulent scalings and structures, can be isolated and studied. For example, Duan, Beekman & Martín (Reference Duan, Beekman and Martín2011) found that many of the scaling relations, though derived based on the Morkovin hypothesis, remained valid for boundary-layer flows with 
${{Ma}}_\infty$ up to 12. In addition, recent work of Di Renzo, Fu & Urzay (Reference Di Renzo, Fu and Urzay2020), Di Renzo & Urzay (Reference Di Renzo and Urzay2021) and Passiatore et al. (Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022) demonstrated the applicability of these scaling relations on high-enthalpy boundary layers with thermochemical non-equilibrium effects.
Despite much progress, the parameter ranges of compressible turbulent DNS databases are still quite limited, especially for channel flows. The highest 
${{Ma}}_b$ reported for a channel flow is approximately 4 (Trettel & Larsson Reference Trettel and Larsson2016). The 
${{Re}}_\tau$ of that case is near 1000, but the centreline friction Reynolds number in semi-local units, 
${{Re}}_{\tau,c}^*$, is only approximately 200. As mentioned above, the highest 
${{Re}}_{\tau,c}^*$ of a channel flow case is reported by Yao & Hussain (Reference Yao and Hussain2020), while the 
${{Ma}}_b$ of that case is limited to 1.5. Any effort to extend current parameter ranges is extremely resource-demanding. For example, the grid number in the case of Yao & Hussain (Reference Yao and Hussain2020) is over 8 billion. If 
${{Ma}}_b$ is increased to 3 at the same 
${{Re}}_\tau$, then the required grid number is evaluated to be dauntingly over 20 billion. The substantial computational cost prohibits a systematic parameter study of various influencing factors.
Compared to instantaneous DNS data, the calculation of only turbulent mean flow can be substantially cheaper, such as through the Reynolds-averaged Navier–Stokes (RANS) or large-eddy simulations (LES) approaches. With only mean profiles in hand, one can still investigate some important turbulence characteristics using linear models, which will be introduced at length in § 1.2. Therefore, the combination of cheap reliable mean-flow calculation and the mean-flow-based analysis can be an effective way to explore the turbulence at high Mach and Reynolds numbers beyond the current capability of DNS. As a demonstration, Cossu, Pujals & Depardon (Reference Cossu, Pujals and Depardon2009) and Hwang & Cossu (Reference Hwang and Cossu2010b) used a curve-fitted incompressible mean flow 
${{Re}}_\tau >10^{4}$, and analysed the transient growth and the responses to forcing for boundary-layer and channel flows. Pickering et al. (Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) proposed a resolvent-based model with the mean flow from LES and captured successfully the characteristics of turbulent coherent structures in the shear flow. This type of method combination will be deployed in this work for canonical compressible channel flows. The focus is on the linear response of mean profiles and coherent structures, as an effort to extend the present knowledge frontier on the turbulence characteristics in high Mach and Reynolds number regimes beyond the scope of existing DNS databases. More specifically, the first aim is to utilize well-established universal relations to obtain turbulent mean profiles with a large parameter space. The second is to develop the linear response solvers to harmonic and stochastic forcing for compressible turbulent flows. The third is to investigate the response characteristics with varying Mach and Reynolds numbers, and provide supporting theoretical analysis.
1.1. Universal relations for compressible wall-bounded turbulent flow
It is well known that incompressible wall-bounded turbulent flows exhibit nearly universal mean streamwise velocity profiles versus the wall-normal coordinate when normalized with the wall viscous units. Therefore, a universal curve-fitted expression is available, and a prevailing one for the channel flow is from Cess (Reference Cess1958). This universal curve is also pursued in compressible turbulent flows for the mean velocity and temperature. For instance, the velocity transformation is designed to transform the compressible streamwise velocity profile to the incompressible counterpart. The pioneering work was accomplished by van Driest (Reference van Driest1951), and subsequent improved versions were proposed by Zhang et al. (Reference Zhang, Bi, Hussain, Li and She2012), Trettel & Larsson (Reference Trettel and Larsson2016), Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020), among others. Very recently, the total-stress-based transformation proposed by Griffin, Fu & Moin (Reference Griffin, Fu and Moin2021b) collapses successfully nearly all available DNS data, with and without heat transfer, for the boundary-layer, channel and pipe flows. The relative errors are typically less than 3 % when 
${{Re}}_\tau > 200$. Note that the performance of this transformation for non-canonical turbulent boundary-layer flows has been reported further in Bai, Griffin & Fu (Reference Bai, Griffin and Fu2022). In terms of the temperature profile, Crocco (Reference Crocco1932) and Busemann (Reference Busemann1931) derived a temperature–velocity relation showing that the mean temperature was nearly a quadratic function of the mean streamwise velocity. Later, a less-restricted relation was proposed by Walz (Reference Walz1969) through introducing the recovery temperature. Furthermore, Duan & Martín (Reference Duan and Martín2011) made an improvement by introducing the mean enthalpy and provided more universal fitting coefficients for the quadratic expression based on the data of wide parameter ranges. This relation is shown to be sufficiently accurate for both boundary-layer and channel flows, even with high-enthalpy effects (Passiatore et al. Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022). Also, Zhang et al. (Reference Zhang, Bi, Hussain and She2014) proposed a general Reynolds analogy for compressible wall-bounded turbulent flows by introducing a general recovery factor.
The success of these universal relations makes it possible to obtain the mean profiles by solving an inverse problem under proper flow conditions in the outer region. This idea was adopted recently by Griffin, Fu & Moin (Reference Griffin, Fu and Moin2022) to realize a highly accurate wall-modelled LES framework and to estimate the grid points required for DNS and LES (Griffin, Fu & Moin Reference Griffin, Fu and Moin2021a). For the channel flow with symmetric isothermal boundaries, Song et al. (Reference Song, Zhang, Liu and Xia2022) proposed an empirical relation between the centreline mean temperature and velocity. This can serve as an outer boundary condition to close the inverse problem, such that the mean-flow profiles can be obtained by solving a simple ordinary differential equation (ODE).
1.2. Linear analysis of turbulent wall-bounded flows
The linear models for analysing turbulent flows have received increasing interest, especially over the last decade (see McKeon Reference McKeon2017). This method of analysis shares some fundamental procedures with linear stability theory, and borrows some critical thoughts from control theory (Taira et al. Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020; Jovanović Reference Jovanović2021). The starting point is the linearized Navier–Stokes equation around a mean state for the Fourier-decomposed perturbation. The model equation can be written as
\begin{equation} \mathcal{L}_{mn}\,{\hat{\boldsymbol{q}}}_{mn} = {\hat{\boldsymbol{f}}}_{mn}, \end{equation}
where 
$\mathcal {L}$ is a linear operator related to the mean flow as well as the temporal and spatial wavenumbers, the subscripts 
$m$ and 
$n$ are the Fourier indexes, and 
${\hat {\boldsymbol {q}}}$ denotes the perturbation shape function. The nonlinear term 
${\hat {\boldsymbol {f}}}$ was ignored in early works, and all the eigenmodes of 
$\mathcal {L}$ were found to be stable for the incompressible turbulent profile (Malkus Reference Malkus1956). Instead of discarding it totally, Reynolds & Hussain (Reference Reynolds and Hussain1972) modelled the nonlinear term into the linear operator by using the Boussinesq approximation and the eddy viscosity 
$\mu _t$. These eddy-viscosity-enhanced linear models have been adopted widely ever since. Better behaviour of artificially induced waves was reproduced in their results by considering 
$\mu _t$, but the mean profile was still stable with 
${{Re}}_\tau \sim 1000$.
As no unstable modes were found, linear model research was quiet for a while until two sets of mathematical tools were applied successfully. The first is the non-modal instability theory with pioneering work from Trefethen et al. (Reference Trefethen, Trefethen, Reddy and Driscoll1993). Specifically, as the eigenvectors of 
$\mathcal {L}$ are not orthogonal to each other, 
${\hat {\boldsymbol {q}}}$ can experience strong transient growth even though each mode is asymptotically stable. This was identified by Butler & Farrell (Reference Butler and Farrell1993) for incompressible channel flows. Furthermore, del Álamo & Jiménez (Reference del Álamo and Jiménez2006) (also Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009) found that by using the eddy-viscosity-enhanced model, two peaks of the transient growth with different spanwise wavenumbers were well related to the near-wall motions and the outer-layer large-scale motions, respectively. The second methodology originates from control theory (Schmid & Henningson Reference Schmid and Henningson2001; Kim & Bewley Reference Kim and Bewley2006; Zare, Georgiou & Jovanović Reference Zare, Georgiou and Jovanović2020), and is called the input–output or resolvent analysis. By rewriting (1.1) as 
${\hat {\boldsymbol {q}}}_{mn} = \mathcal {L}_{mn}^{-1}\,{\hat {\boldsymbol {f}}}_{mn}$ (this is an oversimplified interpretation and the theory details are presented in § 3), 
${\hat {\boldsymbol {f}}}$ is regarded as the input or forcing, and 
${\hat {\boldsymbol {q}}}_{mn}$ is the output or response. Their relation is established through the transfer matrix 
$\mathcal {L}_{mn}^{-1}$, known as the resolvent. Based on artificially selected forcing, the response characteristics of turbulent mean flow can be obtained through singular-value and Karhunen–Loève (KL) decompositions. Hwang & Cossu (Reference Hwang and Cossu2010a,Reference Hwang and Cossub) analysed the linear responses to both harmonic and stochastic forcing for incompressible Couette and channel flows. Similar to the transient growth, the two most amplified structures are the near-wall and large-scale motions of infinite streamwise wavelength. As the most amplified modes can occupy a large portion of the response energy, they represent the characteristic flow structures, which can be deployed for reduced-order modelling. For example, through analysing the most amplified response, McKeon & Sharma (Reference McKeon and Sharma2010) emphasized the vital role of the critical layer for incompressible flows, where the phase velocity of the response equalled the mean velocity. The viscosity effects are significant in the critical layer, so their influence domain can extend well beyond the immediate vicinity of the wall. In the follow-up work, Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013) designed a low-rank model based on the principal response mode. Different classes of scalings for the streamwise turbulence intensity were found.
The results of the above non-modal growth and resolvent analyses are also closely related. The strong resemblance was demonstrated among the amplified structures of transient growth, optimal harmonic and stochastic responses (Hwang & Cossu Reference Hwang and Cossu2010b). Meanwhile, the linearly amplified modes all exhibit self-similarity in the mid-wavelength range related to the logarithmic region (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013; McKeon Reference McKeon2019; Vadarevu et al. Reference Vadarevu, Symon, Illingworth and Marusic2019), which provides mathematical evidence to the well-known attached-eddy model (AEM) initially proposed by Townsend (Reference Townsend1976) (also see the recent review of Marusic & Monty Reference Marusic and Monty2019).
In most of the research mentioned above, the nonlinear forcing term 
${\hat {\boldsymbol {f}}}$ is modelled ideally, instead of an extraction from real flow data. Specifically, 
${\hat {\boldsymbol {f}}}$ is either assumed to be a stochastic signal or taken in an optimal form allowing for the strongest perturbation growth. Consequently, the results using the artificially designed forcing are more qualitative than quantitative, compared with the analysis based on DNS data (Zare, Jovanović & Georgiou Reference Zare, Jovanović and Georgiou2017). The ‘colour’ of the forcing can lead to different weights of the amplified modes, potentially changing the relative importance of different modes. The weights of modes can be obtained by partly matching the DNS data. The resulting energy spectra and other statistics from the resolvent analysis are well consistent with those from DNS (Moarref et al. Reference Moarref, Jovanović, Tropp, Sharma and McKeon2014; McMullen, Rosenberg & McKeon Reference McMullen, Rosenberg and McKeon2020). As another way of improvement, different optimization processes were designed by Zare et al. (Reference Zare, Jovanović and Georgiou2017) and Hwang & Eckhardt (Reference Hwang and Eckhardt2020) to obtain the cospectrum of the forcing signal, which is closer to realistic flows. To gain a complete understanding of the forcing statistics for incompressible turbulence, considerable efforts have been made recently to analyse the forcing based on the DNS data (Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021; Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021; Nogueira et al. Reference Nogueira, Morra, Martini, Cavalieri and Henningson2021). It is revealed that the forcing exhibits a relatively high level of spatial–temporal coherence. The low-rank nature of the forcing will benefit the construction of more reliable reduced-order models of turbulent flow. In addition, Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021) clarified why the introduction of 
$\mu _t$ in the linear operator leads to an improved response prediction compared to that without 
$\mu _t$. That is because the weights or the combination coefficients of the principal resolvent modes are closer to those from DNS. In this sense, the usage of 
$\mu _t$ is also a partial modelling of the forcing colour, consistent with its low-rank nature.
The above results are mainly for incompressible turbulence. The extension and application of linear models to compressible turbulent flows are still limited. Alizard et al. (Reference Alizard, Pirozzoli, Bernardini and Grasso2015) first presented an analysis of linear transient growth on the turbulent boundary layers with 
${{Ma}}_\infty$ up to 4. The inner- and outer-layer modes, analogous to the incompressible ones, are identified, and the weak effects of compressibility are revealed on these characteristic structures. In addition to the similarities, an emerging problem is the dynamics of the suggested streaky temperature perturbation and its coupling with the velocity components. Recently, Bae, Dawson & McKeon (Reference Bae, Dawson and McKeon2020a,Reference Bae, Dawson and McKeonb) applied the resolvent formulation to supersonic turbulent boundary layers with both adiabatic and cold walls. Their results highlight distinct features below and above the relative sonic line of the response. Moreover, there are implications for the geometrically self-similar structure of response modes and the analogous rules for the incompressible counterpart. Dawson & McKeon (Reference Dawson and McKeon2020) investigated the same cases with emphasis on predicting the shapes of the resolvent mode at different Mach numbers. Up to now, to the authors’ knowledge, the analysis framework of stochastic forcing and response has not been applied to supersonic turbulent flows, although there were successful attempts on the stability and receptivity of compressible laminar flows (Alizard et al. Reference Alizard, Gibis, Selent, Rist and Wenzel2022; Madhusudanan & McKeon Reference Madhusudanan and McKeon2022), which demonstrates the feasibility and value of this linear framework in interpreting the physics of compressible flows. Therefore, a systematic study on the linear responses of supersonic turbulent flows to harmonic and stochastic forcing is required, which is the main focus of the present work.
1.3. Streaky motions and compressibility effects
Streamwise streaks are the most common coherent structures across different layers in wall-bounded turbulent flows. Their behaviours in incompressible flows have been studied in great detail, and are introduced briefly below (see e.g. Marusic et al. Reference Marusic, McKeon, Monkewitz, Nagib, Smits and Sreenivasan2010b; Jiménez Reference Jiménez2013, for more details). The streaks are narrow regions elongated in the streamwise direction, containing mainly the fluctuations of streamwise velocity. In the near-wall region, their mean streamwise length and spanwise spacing are approximately 1000 and 100 wall units in the buffer layer, first measured in experiments (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Smith & Metzler Reference Smith and Metzler1983). Furthermore, they are dominant in carrying turbulence energy. The near-wall streaks are also known to be able to self-sustain along with the streamwise vortices, where the linear mechanisms of transient growth and instability play an essential role (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995; Schoppa & Hussain Reference Schoppa and Hussain2002). The streaky motions also exist at larger scales in the outer region, and are more often termed the large-scale or very-large-scale motions (Kim & Adrian Reference Kim and Adrian1999; del Álamo & Jiménez Reference del Álamo and Jiménez2003). For the channel flow, the peak of the pre-multiplied spectra of streamwise velocity corresponds to characteristic streamwise and spanwise scales 
$5h\unicode{x2013}6h$ and 
$h\unicode{x2013}2h$, respectively (del Álamo & Jiménez Reference del Álamo and Jiménez2003; Lee & Moser Reference Lee and Moser2015), where 
$h$ is the channel half-height. These large-scale motions are extremely energetic. Meanwhile, they are least affected by the presence of the wall; on the other hand, they have a modulation effect on the near-wall turbulence intensity (Hutchins & Marusic Reference Hutchins and Marusic2007; Marusic, Mathis & Hutchins Reference Marusic, Mathis and Hutchins2010a).
In compressible flows, additional temperature and density streaks are present due to the wall-normal gradients of the mean flow (see e.g. Coleman et al. Reference Coleman, Kim and Moser1995). Superficial similarities of near-wall streaks were demonstrated in boundary-layer flows for the low, moderate and high Mach number (up to 12) cases (Duan et al. Reference Duan, Beekman and Martín2011). Meanwhile, two crucial features of compressibility effects on the near-wall motion are noted. The first is the variation of the spacing of streaks, on which different trends were reported. For example, wall cooling is shown to increase the streak spacing (Coleman et al. Reference Coleman, Kim and Moser1995; Duan, Beekman & Martin Reference Duan, Beekman and Martin2010), while the increase of 
${{Ma}}$ has a decreasing effect with an adiabatic wall condition (Duan et al. Reference Duan, Beekman and Martín2011). Further analysis shows that the spacing variation is consistent with the change of 
${{Re}}_\tau ^*$. Therefore, the spacing experiences much less variation in semi-local units (Morinishi, Tamano & Nakabayashi Reference Morinishi, Tamano and Nakabayashi2004; Patel, Boersma & Pecnik Reference Patel, Boersma and Pecnik2016; Yao & Hussain Reference Yao and Hussain2020). The second feature is that wall cooling enhances the coherence of the near-wall motions. Meanwhile, the vortical structures become less chaotic (Coleman et al. Reference Coleman, Kim and Moser1995; Morinishi et al. Reference Morinishi, Tamano and Nakabayashi2004; Duan et al. Reference Duan, Beekman and Martin2010; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022). In comparison, the large-scale motions in the outer layer are less influenced by the compressibility effects from both experimental and numerical observations (Spina et al. Reference Spina, Smits and Robinson1994; Pirozzoli, Bernardini & Grasso Reference Pirozzoli, Bernardini and Grasso2008; Williams et al. Reference Williams, Sahoo, Baumgartner and Smits2018; Cheng & Fu Reference Cheng and Fu2022a).
The linear models of turbulence have long been used to investigate streaky motions. They are especially suitable for studying the structure, coherence and instability of streaks (see e.g. Schoppa & Hussain Reference Schoppa and Hussain2002; del Álamo & Jiménez Reference del Álamo and Jiménez2006; Hwang & Cossu Reference Hwang and Cossu2010b; McKeon & Sharma Reference McKeon and Sharma2010). Nevertheless, as introduced in § 1.2, these linear models have rarely been deployed for compressible turbulent flows. The combination of the ODE mean-flow solver and the linear response analysis developed in this work enables a systematic study of the compressibility effects on the streaky motions in the channel flow. The parameter ranges considered are 
${{Ma}}_b$ from 0 to 5, and 
${{Re}}_\tau$ up to 
$10^4$, both of which are far beyond the current DNS database.
The remainder of the article is organized as follows. Section 2 describes the governing equations and the approaches to calculating the mean turbulent flow. Section 3 provides the formulation of the linear response analysis to both harmonic and stochastic forcing. The main results are discussed in § 4, and some theoretical proofs are provided in § 5. Finally, the work is summarized in § 6.
2. Governing equations and mean-flow calculations
2.1. Compressible Navier–Stokes equations
We consider a canonical compressible turbulent channel flow with symmetric isothermal boundaries. The flow illustration and coordinate set-up are shown in figure 1. With the assumption of a calorically perfect gas, the non-dimensional Navier–Stokes equations are written as
$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho {\boldsymbol{u}}\right) = 0, \end{gather}$$
$$\begin{gather}\rho \left(\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}\right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot}\left[\frac{\mu}{{{Re}}_b} \left( {\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} {{\boldsymbol{u}}^{\rm T}}} \right) \right] - \frac{2}{3}\,\boldsymbol{\nabla} \left( \frac{\mu}{{{Re}}_b}\,{\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}} \right), \end{gather}$$
\begin{gather} \rho \left(\frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} T \right) - {{Ec}} \left(\frac{\partial p}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} p \right) = \frac{{{Ec}}\,\mu}{{{Re}}_b} \left[ \boldsymbol{\nabla} \boldsymbol{u} \boldsymbol{:} \left( {\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} {{\boldsymbol{u}}^{\rm T}}} \right) - \frac{2}{3} \left( {\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}} \right)^2 \right]\nonumber\\ + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\kappa}{{{Re}}_b}\, \boldsymbol{\nabla} T\right), \end{gather}
\begin{gather} p = \frac{\rho T}{\gamma {{Ma}}_b^2}, \end{gather}
where 
$\rho$, 
$\boldsymbol {u}=[u,v,w]^{\rm T}$, 
$p$ and 
$T$ are the fluid density, velocity, pressure and temperature, respectively, and 
$\mu$ and 
$\kappa$ are the viscosity and conductivity. For deriving (2.1), the following non-dimensionalizations are employed:
\begin{equation} \left.\begin{gathered} {\boldsymbol{x}} = \dfrac{{\boldsymbol{x}}^d}{h^d}, \quad \boldsymbol{u} = \dfrac{\boldsymbol{u}^d}{U_b^d}, \quad t = \dfrac{t^d U_b^d}{h^d}, \quad \rho = \dfrac{\rho^d}{\rho_b^d}, \quad p = \dfrac{p^d}{\rho_b^d U_b^{d2}}, \quad T = \dfrac{T^d}{T_w^d}, \\ \mu = \dfrac{\mu^d}{\mu_w^d}, \quad \kappa = \dfrac{\kappa^d}{\mu_w^d c_p^d}. \end{gathered}\right\} \end{equation}The resulting Mach, Reynolds, Prandtl and Eckert numbers are defined as
\begin{align} {{Ma}}_b = \dfrac{U_b^d}{a_w^d}, \quad {{Re}}_b = \dfrac{\rho_b^d U_b^d h^d}{\mu_w^d}, \quad {{Pr}} = \dfrac{\mu^d c_p^d}{\kappa^d}=0.72, \quad {{Ec}} = \dfrac{U_b^{d2}}{c_p^d T_w^d} = (\gamma-1)\,{{Ma}}_b^2, \end{align}
where the superscript 
$d$ stands for dimensional quantities, and the subscript 
$w$ denotes quantities at the wall (
$y=\pm 1$). Also, 
$h$ is the channel half-height, and 
$c_p$ is the isobaric specific heat. The speed of sound is 
$a^d=(\gamma R^d T^d)^{1/2}$, with 
$R$ the gas constant, and the specific heat ratio 
$\gamma =1.4$. The bulk density 
$\rho _b^d$ and bulk velocity 
$U_b^d$ are selected to meet the requirement
\begin{equation} \int_{{-}1}^1 {\bar{\rho} \, {\mathrm{d}} y} = \int_{{-}1}^1 {\bar{\rho} \bar{u} \, {\mathrm{d}} y} = 2, \end{equation}where the overbar denotes a mean variable. The viscosity is obtained through Sutherland's law, where the fitting constant is 110.4 K and the reference temperature is 293.15 K.
Schematic of the coordinate set-up for the channel flow.
In wall-bounded turbulent flows, it is common to express the variables in wall viscous units with a superscript 
$+$ as
\begin{equation} {\boldsymbol{x}}^+ = \dfrac{{\boldsymbol{x}}^d}{\delta_\nu^d}, \quad \boldsymbol{u}^+ = \dfrac{\boldsymbol{u}^d}{u_\tau^d}, \quad \rho^+ = \dfrac{\rho^d}{\rho_w^d}, \quad T^+ = \dfrac{T^d}{T_w^d} = T, \quad \mu^+ = \dfrac{\mu^d}{\mu_w^d} = \mu, \end{equation}
where 
$u_\tau ^d = \sqrt {\tau _w^d/\rho _w^d}$ is the friction velocity, 
$\tau _w$ is the wall shear, and 
$\delta _\nu ^d=\mu _w^d/(\rho _w^d u_\tau ^d)$ is the length unit. The corresponding friction Reynolds number is 
${{Re}}_\tau =h^d/\delta _v^d$. Furthermore, Huang, Coleman & Bradshaw (Reference Huang, Coleman and Bradshaw1995) suggested a set of semi-local units for a better collapse of the turbulence statistics with the incompressible counterparts. The velocity and length units are based on the local mean density and viscosity as 
$u_\tau ^{*d}=\sqrt {\tau _w^d/\bar {\rho }^d}$ and 
$\delta _\nu ^{*d}=\bar {\mu }^d/(\bar {\rho }^d u_\tau ^{*d})$, hence 
$y^*=y_d/\delta _v^{*d}$ and 
${{Re}}_\tau ^*=h^d/\delta _v^{*d}$. Hereinafter, the quantities using semi-local units are expressed with a superscript *.
The basic variable of (2.1) is 
$\boldsymbol {q}=[\rho,u,v,w,T]^{\textrm {T}}$, which is decomposed into a mean part and a perturbed part as
\begin{equation} \boldsymbol{q} \left(x,y,z,t\right) = \bar{\boldsymbol{q}} \left(y\right) + \tilde{\boldsymbol{q}} \left(x,y,z,t\right), \end{equation}
where the mean component in the channel flow is taken as 
${\bar {\boldsymbol {q}}}=[\bar {\rho }(y), \bar {U}(y),0,0,\bar {T}(y)]^{\textrm {T}}$. The calculation of 
${\bar {\boldsymbol {q}}}$ is described in § 2.2, and the operations with 
$\tilde {\boldsymbol {q}}$ will be elucidated in § 3.
2.2. Calculation of turbulent mean flow
For the incompressible channel flow, a semi-empirical expression proposed by Cess (Reference Cess1958) is used widely to describe the mean velocity profile 
$\bar {U}_{inc}^+$ at different 
${{Re}}_\tau$ as
\begin{equation} \left.
\begin{array}{l@{}} \displaystyle \bar{U}_{inc}^+
(y^+;{{Re}}_\tau) = \int_{0}^{y^+} {\dfrac{1 -
y^+{/}{{Re}}_\tau}{1 +
\nu_t}\,{\mathrm{d}} y^+} ,\\
\displaystyle \nu_t = \dfrac{1}{2} \sqrt{1 + f_{\nu_t}^2} -
\dfrac{1}{2}, \quad f_{\nu_t} = \dfrac{\kappa_c\,{{Re}}_\tau}{3}\,(1-y^2)
(1+2y^2) \left[1 - \exp
\dfrac{{{Re}}_\tau\,(|y|-1)}{A_{\nu_t}}
\right] , \end{array} \right\}
\end{equation}
where 
$\nu _t=\mu _t/\bar {\rho }$ is the kinematic eddy viscosity, and 
$\kappa _c$ and 
$A_{\nu _t}$ are two fitting parameters. Based on the mean profile of 
${{Re}}_\tau =2003$, the values recommended by del Álamo & Jiménez (Reference del Álamo and Jiménez2006) are 
$\kappa _c=0.426$ and 
$A_{\nu _t}=25.4$, respectively. By including the DNS dataset of 
${{Re}}_\tau$ up to 5200, the values adopted in the present work are 
$\kappa _c=0.418$ and 
$A_{\nu _t}=25.1$ for all 
${{Re}}_\tau$ (see more details in Appendix A). The analytical equation (2.7) has been utilized widely in the linear and resolvent analyses of incompressible flows over wide parameter ranges (see e.g. Butler & Farrell Reference Butler and Farrell1993; del Álamo & Jiménez Reference del Álamo and Jiménez2006; Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013). However, no such expressions are currently available for the compressible channel flow.
As introduced in § 1.1, there are already adequate universal relations for solving the inverse problem. Consequently, the mean profile can be obtained quickly by solving an ODE. The algorithm involves four universal relations: (1) the incompressible velocity profile as shown in (2.7); (2) the compressible velocity transformation from Trettel & Larsson (Reference Trettel and Larsson2016) or Griffin et al. (Reference Griffin, Fu and Moin2021b); (3) the algebraic relation from Duan & Martín (Reference Duan and Martín2011) between the temperature and velocity distributions; (4) the outer boundary condition from Song et al. (Reference Song, Zhang, Liu and Xia2022) connecting the centreline temperature and velocity. Only two parameters are required as the inputs, i.e. 
${{Ma}}_b$ and 
${{Re}}_\tau$ (or equivalently, 
${{Re}}_b$, 
${{Re}}_\tau ^*$ or others), and the solution takes less than one second on a laptop. The algorithm and the solver verification are elaborated in Appendix A.
With this efficient solver in hand, the mean turbulent flow of wide ranges of 
${{Re}}_\tau$ and 
${{Ma}}_b$ can be obtained. As a demonstration, figure 2 depicts the contours of 
${{Re}}_{\tau,c}^*$ and 
$\bar {T}_c$ as functions of 
${{Re}}_\tau$ and 
${{Ma}}_b$, with subscript 
$c$ denoting the channel centre. The black contour with 
${{Re}}_{\tau,c}^*=200$ denotes the lower available boundary of the solver. As can be seen, at a fixed 
${{Re}}_\tau$, 
${{Re}}_{\tau,c}^*$ drops rapidly with the increase of 
${{Ma}}_b$, which heavily restricts the range of 
${{Re}}_{\tau,c}^*$ of the DNS database. The present solver enables us to explore the flow characteristics far beyond the parameter range of the DNS database. Although its accuracy cannot be directly proved rigorously, the following three facts further strengthen our confidence in high 
${{Ma}}_b$ and 
${{Re}}_\tau$ regimes.
(i) The velocity transformation (see (A2)) has been verified for flows with
${{Ma}}_b$ up to 14, and tends to be more accurate at higher ${{Re}}_{\tau,c}^*$ (Trettel & Larsson Reference Trettel and Larsson2016; Griffin et al. Reference Griffin, Fu and Moin2021b).(ii) The accuracy of the velocity–temperature relation (see (A3)) has been proved for flows with
${{Ma}}_b$ up to 14 (Duan & Martín Reference Duan and Martín2011).(iii) The
$\bar {T}_c$ shown in figure 2(b) is insensitive to ${{Re}}_\tau$ when ${{Re}}_{\tau,c}^*>200$. Meanwhile, $\bar {T}_c$ increases at a comparable speed with the adiabatic wall temperature (Song et al. Reference Song, Zhang, Liu and Xia2022).
Contours of (a) the friction Reynolds number in semi-local units 
${{Re}}_{\tau,c}^*$, and (b) the centreline temperature. The black line in (b) is the contour line of 
${{Re}}_{\tau,c}^*=200$.
It is worth mentioning that the ratio of 
$\bar {T}_c/T_w$ increases quickly with 
${{Ma}}_b$ rise. For example, 
$\bar {T}_c/T_w$ is over 5.2 at 
${{Ma}}_b=5$. This results in a relatively large temperature gradient, and hence strong heat transfer at the wall. Therefore, the supersonic channel flow case can be considered as a highly cooled wall case in analogy to the boundary-layer flow. In this sense, the increase of 
${{Ma}}_b$ in the channel flow leads to a wall-cooling effect simultaneously, which will be discussed further in § 4.3. It is also important to note that since 
$T_w<\bar {T}_c$, 
${{Re}}_{\tau }^*$ decreases with 
${{Ma}}_b$ rise, as shown in figure 2(a), contrary to the trend of an adiabatic or moderately cooled wall boundary-layer case with 
$T_w>T_\infty$.
3. Formulation of the linear response analysis
3.1. Linearized Navier–Stokes equation
To study the characteristics of 
$\tilde {\boldsymbol {q}}$, (2.1) is linearized around the mean state 
$\bar {\boldsymbol {q}}$. The resulting equation of 
$\tilde {\boldsymbol {q}}$ is written into a matrix form as
\begin{align} &\boldsymbol{\mathsf{F}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial t} + \boldsymbol{\mathsf{A}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial x} + \boldsymbol{\mathsf{B}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial y} + \boldsymbol{\mathsf{C}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial z} + \boldsymbol{\mathsf{D}} \tilde{\boldsymbol{q}} \nonumber\\ &\quad= \boldsymbol{\mathsf{H}}_{xx}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x^2} + \boldsymbol{\mathsf{H}}_{yy}\,\frac{{\partial ^2} \tilde{\boldsymbol{q}}}{\partial y^2} +\boldsymbol{\mathsf{H}}_{zz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial z^2} + \boldsymbol{\mathsf{H}}_{xy}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x\,\partial y} + \boldsymbol{\mathsf{H}}_{yz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial y\,\partial z} + \boldsymbol{\mathsf{H}}_{xz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x\,\partial z} + {\tilde{\boldsymbol{N}}}, \end{align}
where the first eleven terms are linear terms of 
$\tilde {\boldsymbol {q}}$, and 
${\tilde {\boldsymbol {N}}}$ collects all the nonlinear terms. The coefficients 
$\boldsymbol{\mathsf{F}}$, 
$\boldsymbol{\mathsf{A}}$, 
$\boldsymbol{\mathsf{B}}$, 
$\boldsymbol{\mathsf{C}}$, 
$\boldsymbol{\mathsf{D}}$ and 
$\boldsymbol{\mathsf{H}}$ are all 
$5\times 5$ matrices related only to 
${\bar {\boldsymbol {q}}}$. Their specific expressions are listed in Appendix C.
In a turbulent flow, 
$\tilde {\boldsymbol {q}}$ and 
${\tilde {\boldsymbol {N}}}$ can have large amplitudes. To facilitate the prediction capability of the linear model, turbulence models are used for the linearization of 
${\tilde {\boldsymbol {N}}}$, such as the widely used eddy-viscosity-enhanced linear model. The introduction of 
$\mu _t$ leads to better results in predicting the response growth, because it partially models the ‘colour’ of the nonlinear forcing term (Moarref & Jovanović Reference Moarref and Jovanović2012; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019, Reference Morra, Nogueira, Cavalieri and Henningson2021). In compressible flows, however, the expression of 
${\tilde {\boldsymbol {N}}}$ is far more complicated than its incompressible counterpart. Following Alizard et al. (Reference Alizard, Pirozzoli, Bernardini and Grasso2015) (also Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021), two substitutions in (3.1) are made for the turbulence modelling, i.e.
\begin{equation} \bar{\mu} \rightarrow \bar{\mu} + \mu_t, \quad \bar{\kappa} \rightarrow \bar{\kappa} + {\kappa}_t = \frac{\bar{\mu}}{{{Pr}}} + \frac{\mu_t}{{{Pr}}_t}. \end{equation}
Here, 
$\kappa _t$ and 
${{Pr}}_t$ are the turbulent heat conductivity and Prandtl number, and 
$\mu _t$ and 
$\kappa _t$ are non-dimensionalized in the same way as 
$\mu$ and 
$\kappa$ in (2.2). The assumptions behind (3.2a,b) are the Boussinesq assumption and the classic Reynolds analogy. The usage of the former is analogous to the incompressible counterpart. The deployment of the latter assumes implicitly that subject to the modelled part of the forcing, the diffusion of the internal energy of the response is proportional to the Reynolds stress part. We give more remarks on the usage of (3.2a,b). The two assumptions above may not be very accurate throughout the flow field. Nevertheless, the introduction of 
$\mu _t$ and 
$\kappa _t$ can improve the prediction performance of the linear model as long as the modelling points to the ‘right direction’. This is satisfied here for the canonical channel flow, which is away from the destruction of separation bubbles, shock/boundary-layer interactions, shedding vortex, etc. In this relatively ‘clean’ turbulence, the RANS models have proven their capability extensively (Wilcox Reference Wilcox2006). Furthermore, as noted by Cossu et al. (Reference Cossu, Pujals and Depardon2009) and will be shown in § 3.4, the results of the most amplified response are not very sensitive to the shapes of 
$\mu _t$ and 
${{Pr}}_t$.
3.2. Responses to harmonic and stochastic forcing
As the mean flow is homogeneous in the streamwise and spanwise directions, the following Fourier decomposition is applied on 
$\tilde {\boldsymbol {q}}$ as
\begin{equation} \tilde{\boldsymbol{q}} (x,y,z,t) = \iint_{-\infty}^{\infty} {{\hat{\boldsymbol{q}}}(y, t) \exp \left[{\mathrm{i}} (k_x x + k_z z)\right] {\mathrm{d}} k_x \,{\mathrm{d}} k_z}, \end{equation}
where 
$k_x$ and 
$k_z$ are the streamwise and spanwise wavenumbers, respectively, and 
${\hat {\boldsymbol {q}}}$ is the shape function. The same decomposition is performed on 
${\tilde {\boldsymbol {N}}}$. The resulting 
${\hat {\boldsymbol {N}}}$ involves the convolution of 
${\hat {\boldsymbol {q}}}$ components due to the nonlinearity.
After substituting (3.3) into (3.1), the following equation is arrived at for a single mode with 
$k_x\ne 0$ or 
$k_z\ne 0$:
\begin{equation} \boldsymbol{\mathsf{F}}\,\frac{\partial {\hat{\boldsymbol{q}}}}{\partial t} + \hat{\boldsymbol{D}} {\hat{\boldsymbol{q}}} + {\hat{\boldsymbol{B}}}\,\frac{\partial {\hat{\boldsymbol{q}}}}{\partial y} + {\hat{\boldsymbol{C}}}\,\frac{\partial^2 {\hat{\boldsymbol{q}}}}{\partial y^2} = {\hat{\boldsymbol{N}}}, \end{equation}where the rearranged matrices are
\begin{equation} \left.
\begin{array}{l@{}} {\hat{\boldsymbol{D}}} =
\boldsymbol{\mathsf{D}} + {\mathrm{i}} k_x
\boldsymbol{\mathsf{A}} + {\mathrm{i}} k_z
\boldsymbol{\mathsf{C}} + k_x^2
\boldsymbol{\mathsf{H}}_{xx} + k_x k_z
{\boldsymbol{\mathsf{H}}_{xz}} +
k_z^2\boldsymbol{\mathsf{H}}_{zz},\\
{\hat{\boldsymbol{B}}} = \boldsymbol{\mathsf{B}} -
{\mathrm{i}} k_x
\boldsymbol{\mathsf{H}}_{xy} - {\mathrm{i}}
k_z \boldsymbol{\mathsf{H}}_{yz},\\
{\hat{\boldsymbol{C}}}
={-}\boldsymbol{\mathsf{H}}_{yy}. \end{array} \right\}
\end{equation}
Equation (3.4) can be rewritten into a standard form as
\begin{equation} \frac{\partial {\hat{\boldsymbol{q}}}}{\partial t} = \boldsymbol{L} {\hat{\boldsymbol{q}}} + {\hat{\boldsymbol{f}}}, \quad \boldsymbol{L} ={-}\boldsymbol{\mathsf{F}}^{{-}1}\hat{\boldsymbol{D}} - \boldsymbol{\mathsf{F}}^{{-}1}\boldsymbol{\hat B}\, \frac{\partial}{\partial y} - \boldsymbol{\mathsf{F}}^{{-}1}\hat{\boldsymbol{C}}\,\frac{\partial^2}{\partial y^2}, \quad {\hat{\boldsymbol{f}}} = \boldsymbol{\mathsf{F}}^{{-}1} \hat{\boldsymbol{N}}, \end{equation}
where 
$\boldsymbol {L}$ is a linear operator, and 
${\hat {\boldsymbol {f}}}=[\hat {f}_\rho,\hat {f}_u,\hat {f}_v,\hat {f}_w,\hat {f}_T]^{\textrm {T}}$ is the forcing term.
Before discussing the energy amplification of responses, an energy norm of 
${\hat {\boldsymbol {q}}}$ (E) needs to be defined. For compressible flows, a widely used form was proposed by Chu (Reference Chu1965) as
\begin{align} \|{\hat{\boldsymbol{q}}}\|^2 = ( {\hat{\boldsymbol{q}}}, {\hat{\boldsymbol{q}}} )_E = \int_{{-}1}^{1} { \left( \bar{\rho} {\hat{\boldsymbol{u}}}^{H} {\hat{\boldsymbol{u}}} + \frac{R \bar{T}}{\bar{\rho}}\,\hat{\rho}^{\dagger} \hat{\rho} + \frac{\bar{\rho} c_v}{{{Ec}}\, \bar{T}}\,\hat{T}^{\dagger} \hat{T} \right) {\mathrm{d}} y} = \int_{{-}1}^{1} {{\hat{\boldsymbol{q}}}^{H} \boldsymbol{M} {\hat{\boldsymbol{q}}} \, {\mathrm{d}} y} \equiv 2 \hat{E}, \end{align}
under the current definition of non-dimensionalization. Here, the superscript 
${\dagger}$ stands for the complex conjugate, 
$c_v=c_v^d/c_p^d=1/\gamma$, 
$\boldsymbol {M}$ is a diagonal matrix, and 
$\hat {E}$ is the energy norm. Compared to the kinetic energy used in incompressible flows, (3.7) additionally considers the contribution from the acoustic and entropy components of 
${\hat {\boldsymbol {q}}}$ (see e.g. Chen, Wang & Fu Reference Chen, Wang and Fu2022). More importantly, from the numerical point of view, 
$\boldsymbol {M}$ is strictly positive definite, which guarantees the invertibility of the transfer matrix (see (3.9a,b)). For later use, we re-express (3.7) after the wall-normal discretization (see § 3.3 for details) in a vector form as
\begin{equation} ( {\hat{\boldsymbol{q}}}, {\hat{\boldsymbol{q}}} )_E = \hat{\boldsymbol{Q}}^{H} \boldsymbol{G}_{MI}\,\hat{\boldsymbol{Q}} \equiv \hat{\boldsymbol{Q}}_{MI}^{H}\,\hat{\boldsymbol{Q}}_{MI}, \quad \hat{\boldsymbol{Q}}_{MI} = \boldsymbol{G}_{MI}^{1/2}\,\hat{\boldsymbol{Q}}, \end{equation}
which is convenient for applying standard matrix functions. Here, the global vector 
$\hat {\boldsymbol {Q}}=[{\hat {\boldsymbol {q}}}_1,{\hat {\boldsymbol {q}}}_2,\ldots,{\hat {\boldsymbol {q}}}_{N_y}]^{\textrm {T}}$ contains all the wall-normal components, the global matrix 
$\boldsymbol {G}_{MI} = \textrm {blkdiag}\{ c_1 \boldsymbol {M}_1, \ldots, c_{N_y} \boldsymbol {M}_{N_y} \}$, with 
$c_j$ the coefficient for numerical integration, and 
$\boldsymbol {G}_{MI}^{1/2}$ is from the Cholesky decomposition of 
$\boldsymbol {G}_{MI}$.
In the following, the responses to harmonic and stochastic forcing are accounted for, in turn. A harmonic forcing takes the form 
${\hat {\boldsymbol {f}}}(y,t)={\check {\boldsymbol {f}}}(y) \exp (-{\mathrm {i}} \omega t)$, where 
$\omega$ is the circular frequency. After a sufficiently long time, the response is also harmonic as 
${\hat {\boldsymbol {q}}}(y,t)={\check {\boldsymbol {q}}}(y) \exp (-{\mathrm {i}} \omega t)$ if the system is linearly stable. From (3.6a–c), the transfer matrix 
$\boldsymbol {H}$ between 
${\check {\boldsymbol {f}}}$ and 
${\check {\boldsymbol {q}}}$ is
\begin{equation} {\check{\boldsymbol{q}}} = \boldsymbol{H} {\check{\boldsymbol{f}}}, \quad \boldsymbol{H} =(-{\mathrm{i}}\omega \boldsymbol{I} - \boldsymbol{\mathsf{L}})^{{-}1}, \end{equation}
where 
$\boldsymbol {H}$ is also known as the resolvent. The energy amplification factor between the forcing and the response is defined as 
$R(\omega,k_x,k_z)=\|\check {\boldsymbol {q}}\|^2 / \|\check {\boldsymbol {f}}\|^2$; then the optimal response is
\begin{equation} R_{max}(k_x,k_z) = R(\omega_{opt};k_x,k_z) = \max_{\omega} \frac{\|\check{\boldsymbol{q}}\|^2}{\|\check{\boldsymbol{f}}\|^2} = \|\boldsymbol{H}\|_\infty^2, \end{equation}
where 
$\omega _{opt}$ is the optimal frequency. Note that 
$R_{max}$ is also called the 
$H_\infty$ norm of the transfer matrix (Zhou, Doyle & Glover Reference Zhou, Doyle and Glover1996). Mathematically, 
$\|\boldsymbol {H}\|_\infty$ is obtained through a singular value decomposition of 
$\boldsymbol {H}$, thus 
$R_{max}=\sigma _1^2$, where 
$\sigma _1$ is the maximum singular value. Consequently, the relation between 
$\check {\boldsymbol {q}}$ and 
$\check {\boldsymbol {f}}$ is
\begin{equation} \check{\boldsymbol{q}} = \sum_j { \sigma_j \boldsymbol{\cdot} ({\check{\boldsymbol{\phi}}}_j, \check{\boldsymbol{f}} )_E \boldsymbol{\cdot} {\check{\boldsymbol{\psi}}}_j}, \end{equation}
where 
$\sigma _1\ge \sigma _j\ge 0$ are the singular values in descending order, and 
${\check {\boldsymbol {\psi }}}_j$ and 
${\check {\boldsymbol {\phi }}}_j$ are the corresponding response and forcing modes (McKeon & Sharma Reference McKeon and Sharma2010), respectively, which satisfy the orthogonal relations
\begin{equation} (\check{\boldsymbol{\psi}}_i,\check{\boldsymbol{\psi}}_j)_E = \delta_{ij}, \quad (\check{\boldsymbol{\phi}}_i,\check{\boldsymbol{\phi}}_j)_E = \delta_{ij}, \end{equation}
where 
$\delta _{ij}$ denotes the Kronecker delta. For the optimal response, the unit forcing is 
$\check {\boldsymbol {\phi }}_1$ and the response is 
$\check {\boldsymbol {\psi }}_1$. As the channel mean flow is symmetric around the centreline, 
$\check {\boldsymbol {\psi }}_j$ and 
$\check {\boldsymbol {\phi }}_j$ are either symmetric or antisymmetric. For the mode of small wavelength, the singular values even appear in pairs (i.e. 
$\sigma _{2j}=\sigma _{2j-1}$), so the shape function is not unique. For uniqueness, the shape function is normalized such that 
$\check {u}_1$ is real at its peak amplitude (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013).
When the system (3.6a–c) is driven by a stochastic forcing, the response is also stochastic. Note that the colour of the nonlinear forcing has been modelled partially by introducing 
$\mu _t$ and 
${{Pr}}_t$, so the remaining part of the forcing, i.e. 
${\hat {\boldsymbol {f}}}$ here, is assumed to be white-in-time. Following previous works (e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Hwang & Cossu Reference Hwang and Cossu2010b), 
${\hat {\boldsymbol {f}}}$ is a 
$\delta$-correlated Gaussian white noise process with zero mean:
\begin{equation} {\hat{\boldsymbol{f}}} = \boldsymbol{B} {\hat{\boldsymbol{f}}}_0, \quad \langle {\hat{\boldsymbol{f}}}_0 \rangle = 0, \quad \langle {\hat{\boldsymbol{f}}}_0(t_1)\, {\hat{\boldsymbol{f}}}_0^{H}(t_2) \rangle = \boldsymbol{I}\,\delta_d (t_1 - t_2), \end{equation}
where 
$\delta _d$ is the Dirac function, 
${\hat {\boldsymbol {f}}}_0$ is the unmasked forcing, and 
$\boldsymbol {B}$ is a mask matrix allowing for the calculation of separate responses to different forcing components (Jovanović & Bamieh Reference Jovanović and Bamieh2005; Madhusudanan & McKeon Reference Madhusudanan and McKeon2022). Specifically, 
$\boldsymbol {B}$ is an identity matrix in normal cases. If, for example, only 
$\hat {f}_T$ – i.e. the temperature component of the forcing – is considered, then 
$\boldsymbol {B}=\boldsymbol {B}_T=\textrm {diag}([0\ 0\ 0\ 0\ 1])$.
The variance of the response 
$V=\langle \|{\hat {\boldsymbol {q}}}\|^2\rangle$ is focused on as a measure of energy amplification, where 
$\langle {\cdot }\rangle$ denotes the ensemble average. Mathematically, 
$V$ is also referred to as the 
$H_2$ norm of the transfer matrix, representing a temporal integral of 
$\textrm {tr}(\boldsymbol {H}\boldsymbol {H}^{H})$ in terms of 
$\omega$ (Zhou et al. Reference Zhou, Doyle and Glover1996), where 
$\textrm {tr}({\cdot })$ is the matrix trace. Equivalently, we have
\begin{equation} V = \langle\hat{\boldsymbol{Q}}_{MI}^{H}\,\hat{\boldsymbol{Q}}_{MI}\rangle = \text{tr}(\boldsymbol{X}_{MI}), \quad \boldsymbol{X}_{MI} = \langle\hat{\boldsymbol{Q}}_{MI}\,\hat{\boldsymbol{Q}}_{MI}^{H}\rangle. \end{equation}
Here, the global covariance tensor 
$\boldsymbol {X}_{MI}$ is obtained by solving the algebraic Lyapunov equation
\begin{align} \boldsymbol{G}_{L,MI}\,\boldsymbol{X}_{MI} + \boldsymbol{X}_{MI}\, \boldsymbol{G}_{L,MI}^{H} + \boldsymbol{G}_{B,MI}\,\boldsymbol{G}_{B,MI}^{H} = {\boldsymbol{0}}, \quad \langle \hat{\boldsymbol{Q}}\,\hat{\boldsymbol{Q}}^{H} \rangle = \boldsymbol{G}_{MI}^{{-}1/2}\,\boldsymbol{X}_{MI}\,\boldsymbol{G}_{MI}^{{-}1/2{H}}, \end{align}
where 
$\boldsymbol {G}_{L,MI}=\boldsymbol {G}_{MI}^{1/2}\,\boldsymbol {G}_{L}\, \boldsymbol {G}_{MI}^{-1/2}$ and 
$\boldsymbol {G}_{L}$ is the global matrix of 
$\boldsymbol {L}$, and a similar definition for 
$\boldsymbol {G}_{B,MI}$. As a result, 
$V$ is obtained using (3.14a,b) without calculating the temporal integral. We can proceed to solve the eigenvalue problem of 
$\langle \hat {\boldsymbol {Q}} \hat {\boldsymbol {Q}}^{H} \rangle$ as
\begin{equation} \langle \hat{\boldsymbol{Q}} \hat{\boldsymbol{Q}}^{H} \rangle\,\boldsymbol{G}_{MI} {\hat{\boldsymbol{Y}}}_j = \theta_j {\hat{\boldsymbol{Y}}}_j, \end{equation}
where 
$\theta _j$ and 
${\hat {\boldsymbol {Y}}}_j$ are the eigenvalue and eigenfunction, and 
${\hat {\boldsymbol {Y}}}_j$ satisfies the same orthogonal relation as (3.12a,b). Note that 
${\hat {\boldsymbol {Y}}}_j$ is also called the KL decomposition or the ‘proper orthogonal decomposition’ mode, representing the leading energy-containing structure of the response. The forcing mode is solved using the back Lyapunov equation. More details can be found in Farrell & Ioannou (Reference Farrell and Ioannou1993), Zhou et al. (Reference Zhou, Doyle and Glover1996) and Bamieh & Dahleh (Reference Bamieh and Dahleh2001).
The shape function of the harmonic forcing will be denoted as 
${\hat {\boldsymbol {q}}}$ not 
$\check {\boldsymbol {q}}$ below, where there is no ambiguity. This is intended for consistency with the notation of the stochastic forcing and thus convenience for writing.
3.3. Numerical consideration
The procedures in § 3.2 require a discretization in the wall-normal direction. The differential matrices are from the Chebyshev collocation point method (Trefethen Reference Trefethen2001). By default, 
$N_y=501$ points are used, which is abundant to ensure grid independence (Bae et al. Reference Bae, Dawson and McKeon2020b). It is worth mentioning that the 
$c_j$ for numerical integration come from those of Clenshaw & Curtis (Reference Clenshaw and Curtis1960), rather than the direct inverse of the differential matrix, to ensure the positive definiteness of 
$\boldsymbol {G}_{MI}$.
The problems of singular value decomposition, Lyapunov equation and eigenvalue, as described in (3.11), (3.15a,b) and (3.16), respectively, are solved using the svd, eig and lyap functions of the Matlab software. In terms of the boundary condition, a no-slip and isothermal wall is assumed on both sides. Therefore, the perturbations at the wall satisfy
\begin{equation} \hat{u}_w = \hat{v}_w = \hat{w}_w = \hat{T}_w = 0 \quad \text{at}\ y={\pm} 1, \end{equation}
and 
$\hat {\rho }_w$ is solved through the perturbed continuity equation. For channel flows, 
$\hat {\rho }_w$ has the same finite amplitudes at both wall sides, and oscillating unphysical modes may appear. These unphysical modes are recognized and removed in the program as they cannot converge with increasing 
$N_y$ (Bewley & Liu Reference Bewley and Liu1998). Verification of the linear response solvers with previous results is provided in Appendix B.
3.4. Eddy viscosity and thermal conductivity
If the turbulent statistics are not available, then 
$\mu _t$ can be evaluated from the mean streamwise momentum (2.1b) as
\begin{equation} 0 ={-}\frac{{\mathrm{d}} \bar{P}}{\mathrm{d}\kern0.06em x} + \frac{{\mathrm{d}}}{{\mathrm{d}} y} \left[ (\bar{\mu} + \mu_t)\, \frac{{\mathrm{d}} \bar{U}}{{\mathrm{d}} y} \right] \quad \rightarrow \quad \mu_t = \frac{1 - y^+{/}{{Re}}_\tau}{{{\mathrm{d}} \bar{U}^+}/{{\mathrm{d}} y^+}} - \bar{\mu}, \end{equation}
the same as that in (2.7). For 
${{Pr}}_t$, its expression by definition is
\begin{equation} {{Pr}}_t = \frac{\overline{\tilde{u}\tilde{v}}\,\partial \bar{T} / \partial y}{\overline{\tilde{v}\tilde{T}} \, \partial \bar{U} / \partial y}, \end{equation}
which is shown to be approximately unity in most regions in the channel flow (Huang et al. Reference Huang, Coleman and Bradshaw1995; Modesti & Pirozzoli Reference Modesti and Pirozzoli2016). Therefore, 
${{Pr}}_t$ can be assumed to be simply a constant. The linear response results of different 
${{Pr}}_t$ are compared below.
To highlight the influence of 
${{Pr}}_t$, the case of a relatively high 
${{Ma}}_b$ of 4 is selected, and 
${{Re}}_\tau ={6000}$. Figure 3(a) gives the distribution of 
$\kappa _t$ with 
${{Pr}}_t=0.8$, 1.0 and 1.2, respectively. Despite the relatively large difference in 
${{Pr}}_t$, the amplification factor of the optimal harmonic response is little affected, as shown in figure 3(b). A visible difference exists only in the region 
$k_zh<1$. The same conclusion applies to the stochastic response (not shown here). For the convenience of calculation, 
${{Pr}}_t=1$ is selected for later use.
(a) Eddy thermal conductivity and (b) energy amplification factor of the optimal harmonic response (
$k_x=0$) based on the mean flow of different turbulent Prandtl number. The case parameters are 
${{Ma}}_b=4$ and 
${{Re}}_\tau ={6000}$.
3.5. Linear response calculation based on the ODE solver
Although the ODE solver can provide an accurate mean flow, the corresponding linear response results may have large deviations from those using the DNS mean flow since there are first- and second-order derivatives of the mean flow in (3.4). Here, the reliability of the ODE-based mean flow in predicting the linear response is examined further.
Two representative cases are selected. The first is the incompressible 
${{Re}}_\tau =5200$ case from Lee & Moser (Reference Lee and Moser2015), whose 
${{Re}}_\tau$ is very high for a channel flow. This helps to examine the present solver when extended to the high-Reynolds-number regime. Note that 
${{Ma}}_b$ is set to 0.01 in the present solver for this case, to avoid singularities. The second case is the 
${{Ma}}_b=1.5$, 
${{Re}}_\tau =1910$ case from Yao & Hussain (Reference Yao and Hussain2020) (no. (15) in table 2 in Appendix A), which has both a high 
${{Re}}_\tau$ and obvious compressibility effects. The energy amplification factors of both optimal harmonic and stochastic responses for the two cases are plotted in figure 4. As can be seen, the results based on the mean flow from the ODE solver match well with those using the DNS mean flow. A small visible difference exists near 
$k_zh\sim 10^{2}$ in figure 4(b). This is attributed to the difference of 
$\bar {U}$ at the junction of the buffer and logarithmic layers, which originates from the deviation of the incompressible curve fit in (2.7) (see Appendix A for more details). In short, figure 4 demonstrates the reliability of the present ODE solver in calculating the linear response of the turbulent mean flow.
Energy amplification factors of the (a,b) optimal harmonic and (c,d) stochastic forcing based on the mean flows from DNS and the ODE solver. Two cases are (a,c) the 
${{Re}}_\tau =5200$ case from Lee & Moser (Reference Lee and Moser2015), and (b,d) the 
${{Ma}}_b=1.5$, 
${{Re}}_\tau =1910$ case from Yao & Hussain (Reference Yao and Hussain2020).
4. Results and discussions
4.1. Basic characteristics of responses
The basic characteristics of the responses to the harmonic and stochastic forcing are analysed first. The case parameters are selected to be 
${{Ma}}_b=2.0$ and 
${{Re}}_\tau =6000$, which is a benchmark case in this work. A parameter study will be conducted in § 4.3.
First, the curves of 
$R_{max}$ and 
$V$ with different 
$k_x$ and 
$k_z$ are plotted in figures 5(a) and 5(b). An overall feature is that the 
$R_\textrm {max}$ and 
$V$ values in the large-scale region (in terms of the spanwise wavelength 
$\lambda _z=2{\rm \pi} /k_x$) are orders of magnitude higher than those in the small-scale region, reflecting the energetic nature of the former (Kim & Adrian Reference Kim and Adrian1999). Meanwhile, 
$R_{max}$ and 
$V$ of the large-scale motion are strongly dependent on 
$k_x$ when 
$\lambda _z>0.3h$. The modes with the largest 
$R_{max}$ and 
$V$ both have an infinite streamwise wavelength (i.e. 
$k_x=0$); for the harmonic forcing, the most amplified mode is also steady (
$\omega _{opt}=0$). In contrast, 
$R_{max}$ and 
$V$ in the small-scale region are very insensitive to 
$k_x$, where the curves of different 
$k_x$ nearly collapse with each other.
(a,b) Energy amplification factors, and (c,d) corresponding pre-multiplied factors of the (a,c) optimal harmonic and (b,d) stochastic forcing in the 
${{Ma}}_b=2$, 
${{Re}}_\tau =6000$ case.
Hwang & Cossu (Reference Hwang and Cossu2010b) found, in their incompressible flow case, that there exist a 
$k_z^{-2}$ scaling and a 
$k_z^{-1}$ scaling in the mid-
$k_z$ range for 
$R_{max}$ and 
$V$, respectively. We also apply these scalings for the compressible case here for examination. A theoretical analysis on the linear operator related to these scalings will be provided in § 5. The pre-multiplied energy amplification factors 
$k_z^2 R_{max}$ and 
$k_z V$ are plotted in figures 5(c) and 5(d). The variation of 
$k_z^2 R_{max}$ bears a strong resemblance to that in the incompressible case, where the classic bimodal structure is observed (del Álamo & Jiménez Reference del Álamo and Jiménez2006; Cossu et al. Reference Cossu, Pujals and Depardon2009). These two peaks correspond to the large-scale motions in the outer layer (termed ‘outer peak’) and the small-scale motions in the near-wall region (termed ‘inner peak’), respectively, as also noted by Alizard et al. (Reference Alizard, Pirozzoli, Bernardini and Grasso2015) using the transient growth analysis. The 
$\lambda _z^+$ of the inner peak, denoted 
$\lambda _{z,i}^+$, is approximately 100, while the 
$\lambda _z$ of the outer peak, denoted 
$\lambda _{z,o}$, decreases with 
$k_x$ rise and has maximum value approximately 
$3.7h$ at 
$k_x=0$. For 
$k_z V$, it also peaks near 
$\lambda _{z,o}=3.8h$ in the large-scale region, but exhibits distinct features from the incompressible case in the small-scale region (see figure 19(b) in Appendix B for reference). The slope of 
$V$ does not decrease but increases at 
$k_zh>400$, so the inner peak is missing in 
$k_z V$. For reference, the variations of 
$\theta _1$ and 
$k_z \theta _1$ are also plotted in figures 5(b) and 5(d). As stated in § 3.2, 
$\theta _1$ is the largest eigenvalue of the covariance matrix, and represents the mode that receives the largest energy gain. As can be seen, 
$k_z \theta _1$ has the same trend of variation as 
$k_z V$ in the large-scale region, but exhibits the inner peak in the small-scale region. In addition, the inner peak also locates at approximately 
$\lambda _z^+=100$, very close to that of 
$k_z^2 R_{max}$. The different behaviour of 
$k_z V$ here from the incompressible case indicates that the compressibility effects can be the cause. This is discussed further below.
The mask matrix 
$\boldsymbol {B}$ (see (3.13a–c)) is used to study the responses to different forcing components. Therefore, 
$\boldsymbol {B}$ is separated into two parts, i.e. the kinematic part and the thermodynamic part (Madhusudanan & McKeon Reference Madhusudanan and McKeon2022), as
\begin{equation} \boldsymbol{B} = \boldsymbol{B}_{\boldsymbol{u}} + \boldsymbol{B}_{\rho T}, \quad \boldsymbol{B}_{\boldsymbol{u}}=\text{diag}([0\ 1\ 1\ 1\ 0]), \quad \boldsymbol{B}_{\rho T}=\text{diag}([1\ 0\ 0\ 0\ 1]). \end{equation}
The variances of the responses to the two forcing parts are denoted 
$V_{\boldsymbol {u}}$ and 
$V_{\rho T}$, respectively. By definition, 
$V_{\boldsymbol {u}}$ results from the forcing applied only to three velocity equations, i.e. 
${\hat {f}}_\rho ={\hat {f}}_T=0$, the same as that in an incompressible flow. Likewise, 
$V_{\rho T}$ results from the forcing components present only in the density and temperature equations, i.e. 
${\hat {f}}_u={\hat {f}}_v={\hat {f}}_w=0$. It is worth mentioning that 
${\hat {f}}_T=0$ does not mean 
$\hat {T}=0$ because of the variable coupling. Figure 6(a) gives the distribution of 
$V$, 
$V_{\boldsymbol {u}}$ and 
$V_{\rho T}$. Approximate linear superposition 
$V \approx V_{\boldsymbol {u}} + V_{\rho T}$ is satisfied with relative error less than 2 %. The shape of 
$k_z V_{\boldsymbol {u}}$ is quite similar to that in the incompressible flow (see figure 19 for reference). It tends to decrease at 
$k_zh>400$, so the inner peak related to the near-wall motion can be recognized, also located at 
$\lambda _z^+\approx 100$. In comparison, 
$k_z V_{\rho T}$ increases monotonically with 
$k_z$. It is one order of magnitude smaller than 
$k_z V_{\boldsymbol {u}}$ around the outer peak, but rises to be the larger one with 
$k_z$ higher than the inner peak (
$\lambda _z^+<100$), which results simultaneously in the quick increase of 
$k_z V$. Therefore, it is the thermodynamic part of the forcing – i.e. the components in the equations of 
$\hat {\rho }$ and 
$\hat {T}$ – that leads to the different trend of 
$k_z V$ in the compressible flow case.
(a) Pre-multiplied energy growth factors, and (b) contribution of the two most energetic KL modes of the stochastic response to different forcing components (
$k_x=0$).
Furthermore, the KL modes of the three covariance tensors corresponding to 
$V$, 
$V_{\boldsymbol {u}}$ and 
$V_{\rho T}$ are computed to see the contribution from the most energetic modes. The contribution is measured by the ratio 
$r_V=(\theta _1+\theta _2)/V$. The second mode (
$\theta _2$) is also included because in the small-
$\lambda _z$ region, the 
$\theta _j$ appear in pairs, i.e. 
$\theta _1=\theta _2$. As shown in figure 6(b), the two maxima of 
$r_{V_{\boldsymbol {u}}}$ are 89 % and 42 %, near the outer and inner peaks, respectively, indicating that a large portion of 
$V_{\boldsymbol {u}}$ is contributed by one mode pair (or two modes). This energetic mode is representative in describing the local turbulent characteristics. In comparison, 
$r_{V_{\rho T}}$ decreases continuously with 
$k_z$, and is less than 3 % at 
$k_zh>380$ (
$\lambda _z^+<100$). This suggests that no single mode dominates in terms of the energy contribution, i.e. the response is of little coherence in the small-
$\lambda _z$ region, though 
$V_{\rho T}$ is large. As 
$r_{V_{\boldsymbol {u}}}$ is much larger than 
$r_{V_{\rho T}}$, 
$\theta _{1,\boldsymbol {u}}$ dominates 
$\theta _{1,\rho T}$ in most regions. Therefore, 
$\theta _1$ nearly equals 
$\theta _{1,\boldsymbol {u}}$, or in other words, the most energetic mode to the total variance is contributed mainly by 
$\theta _{1,\boldsymbol {u}}$. The 
$r_V$ at the outer and inner peaks are 81 % and 21 %, respectively, pulled down by 
$r_{V_{\rho T}}$. The above comparison demonstrates clearly the distinct feature between the responses to the kinematic and thermodynamic components of the forcing. The former is analogous to the incompressible counterpart, exhibiting characteristic large- and small-scale motions. The latter significantly affects the small-scale motion, and is of little coherence.
4.2. Response structure
The shape functions of the modes at the outer and inner peaks (
$k_x=0$) in figure 5 are plotted in figure 7, where both the response and forcing are provided. The wall-normal and spanwise velocities are the two largest components in the input signal, and so are the streamwise velocity and temperature components in the output signal. One significant feature is that the shape functions of the optimal harmonic and stochastic forcing/response are nearly identical, in both shape and amplitude, as also noted by Hwang & Cossu (Reference Hwang and Cossu2010b) in the incompressible case. This is because 
${\hat {\boldsymbol {q}}}$ of the harmonic response is very insensitive to 
$\omega$ at 
$k_x=0$, so 
${\hat {\boldsymbol {q}}}$ is nearly unchanged from a harmonic forcing to a broadband stochastic one. From a physical point of view, the convection term 
${\bar {\boldsymbol {U}}}\boldsymbol {\cdot }\boldsymbol {\nabla }{{\hat {\boldsymbol {q}}}}$ vanishes with 
$k_x=0$, hence the response shape is determined primarily by the local mean-flow gradient. In terms of different components, the peak locations of 
$|\hat {v}|$ and 
$|\hat {w}|$ are different, forming a streamwise vortex to be shown later. The maxima of 
$|\hat {u}|$ and 
$|\hat {T}|$ of the inner-peak mode are at 
$y^+=11$ and 
$y^+=10$, both of which lie within the buffer layer. For the outer-peak mode, both 
$\hat {u}$ and 
$\hat {T}$ penetrate deep into and below the logarithmic region, as labelled by the dashed line at 
$y^+=100$ in figure 5(d) (recall that 
${{Re}}_\tau =6000$). This is indicative of an amplitude modulation effect of the outer-layer motions on the near-wall velocity and temperature (Marusic et al. Reference Marusic, Mathis and Hutchins2010a; Cheng, Shyy & Fu Reference Cheng, Shyy and Fu2022; Yu & Xu Reference Yu and Xu2022). Specifically, the forcing imposed in the outer region (especially 
$\hat {v}$) induces large amplitudes of streamwise velocity and temperature near the wall as the response. Meanwhile, 
$\hat {T}$ penetrates deeper towards the wall than 
$\hat {u}$, owing to the relatively large temperature gradient near the wall.
Shape functions of different components and the energy norm of the (a,b) forcing and (c,d) response at 
$k_x=0$ for the mode at the (a,c) inner peak (
$\lambda _z^+=100$) and (b,d) outer peak (
$\lambda _z=3.7h$) in figure 5.
As discussed in figure 6, the most energetic modes at the outer and inner peaks result mainly from the kinematic parts of the forcing. In this sense, 
$\hat {T}$ behaves more like a passive scalar (noted by Coleman et al. Reference Coleman, Kim and Moser1995), resulting from the interaction between velocity components of the forcing and the mean temperature gradient. This is demonstrated in figure 8, where the contours of 
$\tilde {u}$ and 
$\tilde {T}$ in the 
$y$–
$z$ plane are depicted for optimal harmonic response. The shape of the stochastic response is nearly the same (see figure 7) and is thus not shown here. The velocity vector 
$[\tilde {w},\tilde {v}]$ of the forcing is also plotted to identify the in-plane motion. As can be seen, the classic combination of quasi-streamwise vortices and streamwise streaks is recognized for the modes of both inner and outer peaks. Moreover, the streamwise vortex as the input strongly promotes the streaks of streamwise velocity and temperature as the output through the lift-up effect (Ellingsen & Palm Reference Ellingsen and Palm1975; Schoppa & Hussain Reference Schoppa and Hussain2002). For both modes here, the temperature streaks look more stuck to the wall, as already shown in figure 7.
Contours of the (a,c) streamwise velocity and (b,d) temperature of the optimal harmonic response for the modes of the (a,b) inner peak (
$\lambda _z^+=100$) and (c,d) outer peak (
$\lambda _z=3.7h$) in the 
$y$–
$z$ plane. The velocity vector is based on the spanwise and wall-normal velocities of the forcing.
More discussion is provided on the structure of the outer-peak mode. As shown in figures 8(c) and 8(d), the contour of the response almost fills in the whole spatial domain, and extends to as long as 
$\lambda _{z,o}=3.7h$ in the spanwise direction. This value is very close to that in incompressible channel flows (e.g. Cossu et al. Reference Cossu, Pujals and Depardon2009; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009). Actually, as will be shown in § 4.3, 
$\lambda _{z,o}$ is nearly independent of 
${{Ma}}_b$ and 
${{Re}}_\tau$, standing as a pretty robust feature of the optimal response in the outer layer. However, as noted by previous researchers, such a large-scale structure has not been observed in the DNS data or the experiment. The outer peak of the pre-multiplied energy spectra from DNS indicates a spanwise scale 
$\lambda _z\approx h- 2h$ for both incompressible (Lee & Moser Reference Lee and Moser2015) and compressible (Yao & Hussain Reference Yao and Hussain2020) channel flows. Part of the reason is that the structures in figure 8 correspond to infinite 
$\lambda _x$ (
$k_x=0$), but in actual flows it is hard to extend to that long. According to figure 5(c), 
$\lambda _{z,o}$ decreases at a lower 
$\lambda _x$, to be closer to the DNS data. In addition, nonlinear effects of the response growth can also lead to the difference (Cossu et al. Reference Cossu, Pujals and Depardon2009).
In the incompressible flow case, the 
$\hat {u}$ response with different 
$k_z$ values between those of the inner and outer peaks exhibits a geometrically self-similar structure (Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010b; McKeon Reference McKeon2019). Specifically, the normalized 
$|\hat {u}|$ of different 
$k_z$ can be collapsed approximately using a wall-normal coordinate 
$y^+/\lambda _z^+$. The amplitude peaks lie within the logarithmic layer, so this self-similar structure stands as a supporting proof for the well-known AEM (Townsend Reference Townsend1976). In addition, the outer-layer large-scale motions and near-wall small-scale motions can be regarded as the largest and smallest attached eddies, respectively (Hwang Reference Hwang2015; Cheng et al. Reference Cheng, Li, Lozano-Durán and Liu2019, Reference Cheng, Shyy and Fu2022). Although abundant evidence has appeared to support the AEM in the incompressible turbulence, its behaviour and validity in compressible turbulence are far from clear. Here, we provide some understanding of the compressible AEM by examining the geometrical self-similarity of the response.
The streamwise velocity, temperature and energy norm of the optimal harmonic response are plotted in figure 9, with 
$\lambda _z^+$ ranging from 800 to 4000 (
$\lambda _z/h=0.13- 0.67$). They all exhibit a streamwise-elongated streaky structure, as observed in the high-speed experiment (Ganapathisubramani, Clemens & Dolling Reference Ganapathisubramani, Clemens and Dolling2006). The peaks of 
$|\hat {u}|$ of these modes lie within the logarithmic layer, but those of 
$|\hat {T}|$ are much closer to the wall due to large mean temperature gradient there. To collapse the highly scattering distribution of the shape function, the wall-normal coordinate 
$y^+$ is normalized by 
$\lambda _z^+$. As the amplitude of 
${\hat {\boldsymbol {q}}}$ is determined by (3.12a,b) based on the integral form in (3.7), the amplitude is normalized correspondingly using 
$\lambda _z^+$. As shown in figures 9(d)–9(f), the normalized curves of different 
$\lambda _z^+$ are much more concentrated, and are collapsed to a first approximation. Also, the maximum of the normalized 
$|\hat {T}|$ experiences larger variations with 
$\lambda _z^+$ than those of 
$|\hat {u}|$ and 
$\hat {E}$. For comparison, the shapes of the inner-peak mode (
$\lambda _z^+=100$), as shown in figure 7, are also plotted in black dashed lines. As can be seen, the curve of 
$\lambda _z^+=100$ apparently deviates from others of higher 
$\lambda _z^+$. This is reasonable because its amplitude peak is already within the buffer layer (see figure 7), out of the logarithmic region, which violates the assumptions of the AEM. Another wall-normal coordinate widely used in compressible turbulence is the semi-local one, 
$y^*$. But here, the ratio 
$y^*/\lambda _z^*$ equals 
$y^+/\lambda _z^+$, thus leading to no difference in collapsing curves. As the unit 
$\delta _\nu ^*$ and thus 
$\lambda _z^*$ are not invariant at different 
$y^+$, other attempts at normalization can be deployed to reflect the difference between 
$y^+$ and 
$y^*$. For example, 
$\lambda _z^*(y^+)$ can be replaced with 
$\lambda _z^*(y_m^+)$, where 
$y_m^+$ is a fixed point for each 
$\lambda _z^+$. As a choice of 
$y_m^+$, the peak location of 
$|\hat {u}|$ (or 
$|\hat {T}|$, 
$\hat {E}$), denoted 
$y_{U_m}^+$, can be regarded as a characteristic wall-normal scale of the response. Nevertheless, using 
$y^*/\lambda _z^*(y_{U_m}^+)$ as the wall-normal coordinate results in no visible improvement in collapsing curves (hence not shown here).
Wall-normal distribution of the (a,d) streamwise velocity, (b,e) temperature, and (c,f) energy norm for the optimal harmonic response (
$k_x=0$) with different spanwise wavelengths. The black dashed lines are the corresponding results of the 
$\lambda _z^+=100$ mode. The coordinates in (d–f) are normalized by the spanwise wavelength.
Finally, the response modes of 
$\lambda _z^+<100$ are focused on because they have different behaviours from the incompressible case (see figure 6). Figures 10(a) and 10(b) display two representative cospectra, 
$\langle \hat {u}\hat {T}^{\dagger} \rangle$ and 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle$, for the modes of 
$\lambda _z^+=100$ and 
$\lambda _z^+=20$, respectively. The cospectrum is obtained by extracting the sub-elements of 
$\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle$ in (3.15a,b). Moreover, to study the contributors to these cospectra, the responses to the kinematic and thermodynamic parts of the stochastic forcing are calculated separately, denoted as 
$\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\boldsymbol {u}}$ and 
$\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\rho T}$, respectively. Their relative contributions to the cospectrum are plotted in figures 10(c) and 10(d) in percentages. As can be seen, over 95 % of 
$\langle \hat {u}\hat {T}^{\dagger} \rangle$ is contributed by 
$\langle \hat {u}\hat {T}^{\dagger} \rangle _{\boldsymbol {u}}$ for both modes, indicating the strong coupling between the streaks of streamwise velocity and temperature induced by the streamwise vortices as shown in figure 8. On the other hand, the contribution from the thermodynamic part of the forcing is negligible. With 
$\lambda _z^+$ decreasing, the amplitude of 
$\langle \hat {u}\hat {T}^{\dagger} \rangle$ quickly diminishes (see figure 6), to less than half of 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ at 
$\lambda _z^+=20$. This is reasonable because as the streak motion is closer to the wall, the 
$\hat {v}$ component is restricted by the wall boundary, thus the streamwise vortex is suppressed. Compared with 
$\langle \hat {u}\hat {T}^{\dagger} \rangle$, 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ has two hot zones. One is around the buffer layer as well, resulting from the coupling of the temperature and density streaks. The other is attached to the wall within the viscous layer, contributed mainly by 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle _{\rho T}$. As shown in figure 6, the response of the second hot zone is of little coherence, and is connected to the large mean density gradient there through the perturbed continuity equation. With 
$\lambda _z^+$ decreasing from 100 to 20, the peak amplitude of 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ shifts to the wall, reflecting the strengthening of the relative importance of 
$\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\rho T}$.
(a,b) Component cospectrum in the 
$y^+$–
$y^+$ plane, and (c,d) the relative contribution from the response to different components of the forcing. The two 
$y^+$ in the horizontal and vertical coordinates in (a,b) are for the two variables in 
$\langle {\cdot }\rangle$, respectively. The response mode is 
$k_x=0$, with (a,c) 
$\lambda _z^+=100$ and (b,d) 
$\lambda _z^+=20$. Each pair of plots (contour and bar) in (a,c) and (b,d) represents the same variable as labelled.
4.3. Effects of Reynolds and Mach numbers
The cases of different 
${{Ma}}_b$ and 
${{Re}}_\tau$ are investigated further, with an emphasis on compressibility effects. Figure 11 depicts the energy amplification factors of the optimal harmonic responses for different cases. Both the outer and inner units are deployed to study the characteristics of the outer and inner peaks. For the outer-region motion, one significant feature from figures 11(a) and 11(c) is that the outer peaks all correspond to 
$\lambda _{z,o}\approx 3.6h$ for different 
${{Re}}_\tau$ and 
${{Ma}}_b$. Here, 
${{Re}}_\tau$ is from 
$10^{3}$ to 
$10^{4}$, and 
${{Ma}}_b$ is from 0.1 to 5, both covering a relatively large range. The weak dependence of 
$\lambda _{z,o}$ on 
${{Re}}_\tau$ has already been reported in previous research (see e.g. Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009), attributed to the fact that the large-scale motion in the outer layer cannot directly ‘feel’ the presence of the wall, and is hence least affected (Pirozzoli Reference Pirozzoli2014). Furthermore, the weak dependence of 
$\lambda _{z,o}$ on 
${{Ma}}_b$ is demonstrated in the present work. This trend is consistent with the observation in experiment (Spina et al. Reference Spina, Smits and Robinson1994; Williams et al. Reference Williams, Sahoo, Baumgartner and Smits2018) and from DNS of the energy spectrum at different 
${{Ma}}_b$ (Modesti & Pirozzoli Reference Modesti and Pirozzoli2016; Yao & Hussain Reference Yao and Hussain2020; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022; Cheng & Fu Reference Cheng and Fu2022b). As discussed in figure 6, the large-scale structure corresponding to the outer peak results mainly from the kinematic part of the forcing, not the thermodynamic part. Thereby, its motion nature can partly explain the weak dependence on 
${{Ma}}_b$.
Pre-multiplied energy amplification factors of the optimal harmonic responses (
$k_x=0$) for the cases with (a,b) different 
${{Re}}_\tau$ (
${{Ma}}_b=2$), and (c,d) different 
${{Ma}}_b$ (
${{Re}}_\tau =6000$). Plots (a,c) are in the outer units, and (b,d) are in the inner units. Note that the inset in (b) gives the same results as (b) except that 
${{Ma}}_b=4$.
For the mode of the inner peak, it is observed that the pre-multiplied energy growths of different 
${{Re}}_\tau$ are nearly collapsed using the inner units for the Mach 2 case (see figure 11b), where 
$\lambda _{z,i}^+\approx 100$. This collapse is also achieved when 
${{Ma}}_b$ is increased to 4, as displayed in the inset, except that 
$\lambda _{z,i}^+$ rises to 180. Actually, 
$\lambda _{z,i}^+$ indeed increases with 
${{Ma}}_b$, as shown in figure 11(d). Meanwhile, with the increase of 
${{Ma}}_b$, the 
$k_z^{-2}$ scaling of 
$R_{max}$ in the mid-
$k_z$ range has more significant deviations, and 
$k_z^2 R_{max}$ of the inner peak rises in both outer and inner units. This trend will be discussed further in § 5 by analysing the perturbation equations. Through a comparison with the 
$k_z^2 R_{max}$ of the outer peak, the relative importance of the near-wall motion is enhanced, in terms of amplifying the response energy.
The variation of 
$\lambda _{z,i}^+$ with 
${{Ma}}_b$ is further quantified, as depicted in figure 12(a). The 
$\lambda _{z,i}^+$ evaluated from 
$k_z^+ \theta _1^+$ of the stochastic response is also plotted. As discussed in figure 5, the variation of 
$k_z \theta _1$ bears a strong resemblance to 
$k_z^2 R_{max}$. It is observed that 
$\lambda _{z,i}^+$ increases by over three times with 
${{Ma}}_b$ lifted up from 0.1 to 5, whether evaluated from 
$k_z^+ \theta _1^+$ or 
$k_z^{2+} R_{max}^+$. As noted in § 2.2, the increase of 
${{Ma}}_b$ results simultaneously in a wall-cooling effect for the channel flow. Wall cooling is shown in the DNS results to rapidly increase the spanwise spacing of the near-wall streaks for both boundary-layer and channel flows (Morinishi et al. Reference Morinishi, Tamano and Nakabayashi2004; Duan et al. Reference Duan, Beekman and Martin2010; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022), consistent with the trend in figure 12(a). For more quantitative comparison, available DNS results of channel flows are added. The data are from Yao & Hussain (Reference Yao and Hussain2020) with a maximum 
${{Ma}}_b$ of 1.5. The 
$\lambda _{z,i}^+$ plotted are calculated from the local maximum of their pre-multiplied energy spectra of the streamwise velocity at a fixed 
$y^*=15$. As the spectra data are not publicly available, the points are evaluated by the present authors from their published figures, so error bars are added to display the potential error. There is a systematic difference between the 
$\lambda _{z,i}^+$ from the linear response analysis and the DNS results, which has also been noted and interpreted in previous research, and will be discussed at length later.
Variation of the spanwise wavelength of the inner-peak mode with 
${{Ma}}_b$, expressed using the (a) inner units and (b) semi-local units for both the optimal harmonic and stochastic responses. The results evaluated from the DNS energy spectra at 
$y^*=15$ are also plotted (Yao & Hussain Reference Yao and Hussain2020).
As demonstrated by Yao & Hussain (Reference Yao and Hussain2020), if expressed using semi-local units, 
$\lambda _{z,i}^*$ values are all approximately 120 for the three cases with varying 
${{Ma}}_b$. To apply this scaling, the wall-normal height to determine 
$\delta _\nu ^*$ is selected to be 
$y_{E_m}^+$, i.e. the peak location of 
$\hat {E}$, as introduced in § 4.2. The resulting variation of 
$\lambda _{z,i}^*$ with 
${{Ma}}_b$ is plotted in figure 12(b). It is observed that both 
$\lambda _{z,i}^*$ from the optimal and stochastic responses are nearly invariant with 
${{Ma}}_b$ up to 5, consistent with the trend in the DNS data in terms of the spanwise spacing of the near-wall streaks (Morinishi et al. Reference Morinishi, Tamano and Nakabayashi2004; Patel et al. Reference Patel, Boersma and Pecnik2016; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022). Nevertheless, the current DNS database of channel flows cannot reach a 
${{Ma}}_b$ as high as 5. In this sense, the result in figure 12(b) is predictive, waiting for more reference points at higher 
${{Ma}}_b$ for further examination. From the variation of 
$\delta _v^*$, it is conjectured that the increase of 
$\lambda _{z,i}^+$ with 
${{Ma}}_b$ rise is due mainly to the 
$\bar {T}/T_w$ increase, or equivalently the 
$\bar {\rho }/\rho _w$ decrease.
More discussion is provided on the systematic deviation of 
$\lambda _{z,i}^+$ and 
$\lambda _{z,i}^*$ in figure 12 from the DNS data. The first reason is that only the principal mode (pair) is considered here. Though it explains a large portion of the total energy amplification (see figure 6), secondary modes may also play crucial roles. An approach for improvement is to consider several leading modes, such as the work by McMullen et al. (Reference McMullen, Rosenberg and McKeon2020), which impressively reproduced the energy spectra in an incompressible case. However, the determination of weights for different modes involves the usage of instantaneous DNS data, which limits the prediction capability, especially for compressible flows. Second, the energy amplification factor in figure 11 at a given 
$(k_x, k_z, \omega )$ is a global measure of the response at all wall-normal heights, so it is not a direct analogy to the energy spectra in the DNS, which is fixed at a specific height. As noted by Cossu et al. (Reference Cossu, Pujals and Depardon2009), the 
$\lambda _{z,i}^+$ of the inner peak denotes the ‘most probable’ value of the spanwise spacing, while that from DNS is more like a ‘mean’ value.
The energy norm of the inner-peak mode is displayed in figure 13 for the cases of varying 
${{Re}}_\tau$ and 
${{Ma}}_b$. The usage of 
$\hat {E}$ is a universal measure of the response structure rather than 
$\hat {u}$ or 
$\hat {T}$ alone across different 
${{Ma}}_b$. As discussed in figure 12, 
$y^*$ is used for 
$\hat {E}$ of various 
${{Ma}}_b$, and hence also for the cases of varying 
${{Re}}_\tau$ for consistency. As can be seen, the 
$\hat {E}$ profiles of different 
${{Re}}_\tau$ are nearly collapsed in figures 13(b) and 13(c). In comparison, the 
$\hat {E}$ profiles of different 
${{Ma}}_b$ exhibit visible differences; 
$y_{E_m}^+$ decreases continuously with 
${{Ma}}_b$ rises, which is due to the increasingly large temperature gradient of the mean flow.
Wall-normal distributions of the normalized energy norm of the inner-peak mode for the cases (a) 
${{Re}}_\tau =6000$ and varying 
${{Ma}}_b$, (b) 
${{Ma}}_b=2$ and varying 
${{Re}}_\tau$, and (c) 
${{Ma}}_b=4$ and varying 
${{Re}}_\tau$.
The geometrical self-similarity of the response in the mid-
$k_z$ range has been discussed in § 4.2 for the case 
${{Ma}}_b=2$. Here, the self-similar characteristics are examined further for the cases of different 
${{Ma}}_b$. Figure 14 gives the normalized energy norms of different 
$\lambda _z^+$ with 
${{Ma}}_b=0.1$, 3 and 5, in the same style as in figure 9. The minimum 
$\lambda _z^+$ displayed increases from 500 to 800 with 
${{Ma}}_b$ rise because, as shown in figure 9, the mode with 
$\lambda _z^+$ close to 
$\lambda _{z,i}^+$ has obvious deviations, and 
$\lambda _{z,i}^+$ increases at a higher 
${{Ma}}_b$ (see figure 12a). As can be seen, the energy norms of different 
$\lambda _z$ in the mid-
$\lambda _z$ range all approximately collapsed with 
${{Ma}}_b$ ranging from 0.1 to 5. The above observation demonstrates that the geometrical self-similarity of the response is little influenced by the compressibility effect, which is helpful in understanding the AEM for supersonic flows.
Wall-normal distributions of the normalized energy norm for the cases with 
${{Re}}_\tau ={6000}$ and (a) 
${{Ma}}_b=0.1$, (b) 
${{Ma}}_b=3$, and (c) 
${{Ma}}_b=5$, for the optimal harmonic response (
$k_x=0$) with different spanwise wavelengths.
As shown in figure 6, 
$V_{\boldsymbol {u}}$ and 
$V_{\rho T}$ have totally different behaviours at different 
$k_z$. Their variations with 
${{Ma}}_b$ are investigated further, as plotted in figure 15. The variation of 
$k_z V_{\boldsymbol {u}}$ with 
${{Ma}}_b$ is strongly analogous to 
$k_z^2 R_{max}$ (see figure 11d). Both the inner and outer peaks are present in the curves of different 
${{Ma}}_b$. Furthermore, the 
$\lambda _{z,o}^+$ of 
$V_{\boldsymbol {u}}$ is nearly unchanged corresponding to 
$\lambda _{z,o}=3.6h$, while 
$\lambda _{z,i}^+$ increases continuously with 
${{Ma}}_b$, as already described in figure 12(b). The 
$k^{-1}$ scaling in the mid-
$k_z$ range gradually deviates at high 
${{Ma}}_b$, and the relative importance of the inner peak is enhanced in amplifying the response. In comparison, 
$V_{\rho T}$ experiences a much stronger variation with 
${{Ma}}_b$ than 
$V_{\boldsymbol {u}}$. As 
${{Ma}}_b$ increases, 
$k_z^+ V_{\rho T}^+$ in the large-scale region tends to diminish, while that in the small-scale region increases dramatically. Consequently, for the case 
${{Ma}}_b>0.5$, 
$k_z^+ V_{\rho T}^+$ starts to monotonically increase with the rise of 
$k_z^+$. It is also observed that both 
$V_{\boldsymbol {u}}$ and 
$V_{\rho T}$ are subject to the 
$k^{-1}$ scaling when 
${{Ma}}_b=0.1$ in the mid-
$k_z$ range (approximately 
$k_z^+$ from 0.002 to 0.1), as labelled by the horizontal dashed lines in figure 15. With the increase of 
${{Ma}}_b$, this scaling gradually deviates, especially for 
$V_{\rho T}$. In short, the motion due to the kinematic part of the forcing is much less sensitive to the compressibility effect than the thermodynamic part. With 
${{Ma}}_b$ increasing, the thermodynamic components begin to dominate the total variance of the stochastic response in the small-scale region (roughly 
$\lambda _z^+<60$).
Pre-multiplied variances of the stochastic responses to the forcing of only the (a) velocity components and (b) thermodynamic components in the cases of varying 
${{Ma}}_b$ (
${{Re}}_\tau =6000$).
As discussed in § 4.1, the response corresponding to 
$V_{\rho T}$ is much less coherent. Therefore, for the response to full forcing (
$V$), the increase of 
$V_{\rho T}$ with 
${{Ma}}_b$ can decrease the dominance of the principal mode pair and thus the coherence of the flow. This is examined by plotting the ratio 
$r_V=(\theta _1+\theta _2)/V$ for the cases of different 
${{Ma}}_b$, as displayed in figure 16. The ratio corresponding to the optimal harmonic response, 
$r_\sigma =(\sigma _1^2+\sigma _2^2)/\sum \sigma _j^2$, is also plotted. As can be seen, 
$r_V$ and 
$r_\sigma$ are not affected by 
${{Ma}}_b$ in the large-scale region near the outer peak. The motion there, with 
$r_V$ and 
$r_\sigma$ both larger than 80 %, is dominated by that induced by the kinematic part of the forcing, also insensitive to 
${{Ma}}_b$ and 
${{Re}}_\tau$ (see figure 11). In the small-scale region, 
$r_V$ experiences different trends of variation with 
${{Ma}}_b$ in the ranges larger and smaller than 
$\lambda _z^+=90$, respectively. It is the same for 
$r_\sigma$ except that the dividing line is at 
$\lambda _z^+=130$. Specifically, 
$r_V$ (also 
$r_\sigma$) increases with 
${{Ma}}_b$ for 
$\lambda _z^+>90$, i.e. the principal mode pair has an increasingly large contribution to the total energy growth, which is attributed to the increase of 
$V_{\boldsymbol {u}}^+$. The DNS results also indicate that wall cooling tends to increase the coherence of the near-wall motion in this wavelength range. The interpretation is that the anisotropic motion is suppressed by wall cooling due to the reduction in energy exchange (Coleman et al. Reference Coleman, Kim and Moser1995; Morinishi et al. Reference Morinishi, Tamano and Nakabayashi2004; Duan et al. Reference Duan, Beekman and Martin2010; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022). In contrast, it is observed here that 
$r_V$ (also 
$r_\sigma$) quickly decreases with 
${{Ma}}_b$ at 
$\lambda _z^+<90$, indicating a reduction in coherence and also the dominance of the principal mode pair for the total energy growth. This trend is obviously due to the enhancement of the response induced by the thermodynamic part of the forcing, as shown in figure 15(b).
Contribution from the principal mode pair to the total energy amplification (in %) for the (a) stochastic and (b) optimal harmonic responses in the cases of varying 
${{Ma}}_b$.
5. Theoretical analysis of the linear operator
It is observed in figures 11(c) and 15 that the pre-multiplied energy amplification of the optimal harmonic and stochastic response gradually deviates from the 
$k_z^{-2}$ and 
$k_z^{-1}$ scalings in the mid-
$k_z$ range at higher 
${{Ma}}_b$. We appeal to a theoretical analysis on (3.1) to support this observation.
To simplify the form of the equations, the power law 
$\mu =T^n$ is adopted temporarily for derivation use instead of the original Sutherland's law. This results in only a minor difference because 
$\mu$ is much smaller than 
$\mu _t$. The mode of 
$k_x=0$ is considered, and the terms related directly to 
$\bar {\rho }$ and 
$\bar {T}$ are combined to cluster the effect of compressibility in the equation. For example, the continuity and streamwise momentum equations are written as
\begin{equation} \left.
\begin{array}{l@{}} \displaystyle \dfrac{\partial}{\partial t}
\dfrac{\hat{\rho}}{\bar{\rho}} -
\dfrac{\bar{T}_y}{\bar{T}}\,\hat{v} + \left( \hat{v}_y +
{\mathrm{i}} k_z \hat{w} \right) =
\dfrac{\hat{f}_\rho}{\bar{\rho}} ,\\ \displaystyle
\dfrac{\partial \hat{u}}{\partial t} + \bar{U}_y \hat{v} -
\dfrac{1}{{{Re}}_b} \left[
\dfrac{\tilde{\mu}_{t,y}}{\bar{\rho}}\,\hat{u}_y +
\dfrac{\tilde{\mu}_t}{\bar{\rho}} \left(\hat{u}_{yy} -
k_z^2 \hat{u}\right) \right]\\ \displaystyle \quad{}-
\dfrac{n}{{{Re}}_b}\,\dfrac{\bar{\mu}}{\bar{\rho}}
\left[ \left(\bar{U}_{yy} + (n-1)
\bar{U}_y\,\dfrac{\bar{T}_y}{\bar{T}}\right)
\dfrac{\hat{T}}{\bar{T}} +
\bar{U}_y\,\dfrac{\hat{T}_y}{\bar{T}} \right] = \hat{f}_u ,
\end{array} \right\} \end{equation}
where the subscript 
$y$ denotes the wall-normal derivative, and 
$\tilde {\mu }_t=\mu _t+\bar {\mu }$ is the total viscosity. As can be seen, the mean-flow coefficients are 
$\bar {U}$, 
$\bar {T}$, 
$\tilde {\mu }_t$, 
$\bar {\mu }$ and their wall-normal derivatives. The other three equations can be rearranged in the same way; their expressions are specified in Appendix C. Therefore, if the basic variable 
$\hat{\boldsymbol{q}}$ is replaced by 
$\hat{\boldsymbol{q}}_0$, then the operator form of the equation is
\begin{equation} \frac{\partial \hat{\boldsymbol{q}}_0}{\partial t} = \mathcal{L} \hat{\boldsymbol{q}}_0 + \hat{\boldsymbol{f}}_0, \quad \hat{\boldsymbol{q}}_0 = \left[\frac{\hat{\rho}}{\bar{\rho}}, \hat{u},\hat{v},\hat{w},\frac{\hat{T}}{\bar{T}}\right], \end{equation}and the linear operator is
\begin{equation} \mathcal{L} = \mathcal{L} \left( \bar{U}_y,\bar{U}_{yy}, \frac{\bar{T}_y}{\bar{T}},\frac{\bar{T}_{yy}}{\bar{T}}, \frac{\bar{T}}{{{Ma}}_b^2}, \frac{\tilde{\mu}_t}{\bar{\rho}}, \frac{\tilde{\mu}_{t,y}}{\bar{\rho}}, \frac{\bar{\mu}}{\bar{\rho}} \right). \end{equation}
The variation of these mean-flow coefficients reflects the influencing path of compressibility on the response behaviours. Note that 
$\kappa _t$ and 
$\bar {\kappa }$ are related to viscosity through the constants 
${{Pr}}$ and 
${{Pr}}_t$, so they are not included explicitly in the bracket. In the following, the behaviours of these mean-flow coefficients are analysed for the cases of different 
$k_z$.
The self-similarity of the shape function is in terms of the normalized coordinate 
$y^+/\lambda _z^+$, so a rescaled coordinate is introduced as 
$\check {y}=\bar {k}_z (y+1)$, where 
$\bar {k}_z=k_z/\check {k}_z$, and 
$\check {k}_z$ is an arbitrarily selected reference value. Note that the overhead check in this section represents the rescaled variable, different from that defined in § 3.2 for a shape function The rescaled form of (5.2a,b) is written as
\begin{equation} \frac{\partial \hat{\boldsymbol{q}}_0}{\partial \check{t}} = \tilde{\mathcal{L}} \hat{\boldsymbol{q}}_0 + \frac{1}{\bar{k}_z}\, \hat{\boldsymbol{f}}_0, \end{equation}
where 
$\check {t}=\bar {k}_z t$. If the rescaled operator 
$\tilde {\mathcal {L}}$ is independent of 
$\bar {k}_z$, i.e. 
$\tilde {\mathcal {L}}=\check {\mathcal {L}}$, then Hwang & Cossu (Reference Hwang and Cossu2010b) proved that there are a 
$k_z^{-2}$ scaling for 
$R_{max}$ and a 
$k_z^{-1}$ scaling for 
$V$ in the mid-
$k_z$ range (see figure 5). Conversely, the variable that cannot be made independent of 
$\bar {k}_z$ is responsible for the failure of the scaling.
As 
$\hat {\boldsymbol {q}}_0$ is a solution defined on the 
$\check {y}$-coordinate, its wall-normal derivative follows
\begin{equation} \frac{\partial}{\partial y} = \bar{k}_z\,\frac{\partial}{\partial \check{y}}. \end{equation}
However, this chain rule is not applicable for the derivatives of mean-flow quantities because they are solutions on the 
$y$-coordinate. Taking the streamwise momentum equation as an example, the rescaled form is written as
\begin{equation} \left.
\begin{array}{l@{}} \displaystyle \dfrac{\partial
\hat{u}}{\partial \check{t}} +
\dfrac{\bar{U}_y}{\bar{k}_z}\,\hat{v} -
\dfrac{1}{{{Re}}_b} \left[
\dfrac{\tilde{\mu}_{t,y}}{\bar{\rho}}\,\hat{u}_y +
\dfrac{\bar{k}_z \tilde{\mu}_t}{\bar{\rho}}
\left(\hat{u}_{\check{y}\check{y}} - \check{k}_z^2
\hat{u}\right) \right] \\ \displaystyle \quad{}-
\dfrac{n}{{{Re}}_b}\,\dfrac{\bar{k}_z
\bar{\mu}}{\bar{\rho}} \left[
\left(\dfrac{\bar{U}_{yy}}{\bar{k}_z^2} +
(n-1)\,\dfrac{\bar{U}_y}{\bar{k}_z}\,\dfrac{\bar{T}_y}{\bar{k}_z
\bar{T}}\right) \dfrac{\hat{T}}{\bar{T}} +
\dfrac{\bar{U}_y}{\bar{k}_z}\, \dfrac{\hat{T}_{\check
y}}{\bar{T}} \right] = \dfrac{\hat{f}_u}{\bar{k}_z}.
\end{array} \right\} \end{equation}
The mean-flow coefficients are analysed as follows. The peaks of the perturbation in figure 9 are within the logarithmic region, where the streamwise velocity from the current velocity transformation (see (A2b)) satisfies
\begin{equation} \bar{\mu}\,\frac{\partial \bar{U}^+}{\partial y^+} = \frac{{\mathrm{d}} \bar{U}_{TL}^+}{{\mathrm{d}} y^*} = \frac{1}{\kappa_c y^*} \quad \rightarrow \quad \bar{U}_{y} = \frac{u_\tau}{\kappa_c}\,\frac{\sqrt{\bar{T}}}{y+1}, \end{equation}
where 
$\kappa _c$ is the Kármán constant. Therefore, the rescaled velocity gradient is
\begin{equation} \frac{\bar{U}_y(y)}{\bar{k}_z} = \frac{u_\tau}{\kappa_c}\, \frac{\sqrt{\bar{T}(y)}}{\check{y}} = \check{\bar U}_{\check{y}}\, \frac{\sqrt{\bar{T}(y)}}{\sqrt{\bar{T}(\check{y})}} \equiv K_1(y,\check{y})\,\check{\bar U}_{\check{y}}. \end{equation}
As can be seen, only if 
$K_1=1$, i.e. 
$\bar {T}(y)$ is constant, is 
${\bar {U}_y}/{\bar {k}_z}=\check {\bar U}_{\check {y}}$ independent of 
$\bar {k}_z$ in the logarithmic region. Similarly, the temperature gradient from (A3) is
\begin{equation} \frac{\bar{T}_y}{\bar{T}} = \frac{ C_1 \left[2 (C_2 - C_3) \bar{U} + C_3 \bar{U}_c \right]}{1 + C_1 \bar{U} \left[(C_2 - C_3) \bar{U} + C_3 \bar{U}_c \right]}\,\bar{U}_y \equiv C_T(y)\,\bar{U}_y, \end{equation}so the rescaled temperature gradient is
\begin{equation} \frac{\bar{T}_y(y)}{\bar{k}_z\,\bar{T}(y)} = C_T(y)\, \frac{\bar{U}_y(y)}{\bar{k}_z} = K_2\,\frac{\check{\bar T}_{\check y}(\check{y})}{\check{\bar T}(\check{y})}, \quad \text{where } K_2(y,\check{y}) = \frac{\sqrt{\bar{T}(y)}\, C_T(y)}{\sqrt{\bar{T}(\check{y})}\,C_T(\check{y})}. \end{equation}For the eddy viscosity, (3.18) in the logarithmic region gives
\begin{equation} \frac{\tilde{\mu}_t}{\bar{\rho}} = \frac{\kappa_c\,{{Re}}_\tau}{\rho_w}\, \sqrt{\bar{T}}\,(y+1) \left[ 1 - (y+1) \right], \end{equation}so the rescaled version is
\begin{equation} \frac{\bar{k}_z\,\tilde{\mu}_t(y)}{\bar{\rho}(y)} \approx \frac{\kappa_c\, {{Re}}_\tau}{\rho_w}\,\sqrt{\bar{T}}\,\check{y} \left( 1 - \check{y} \right) = K_1\,\frac{\check{\tilde{\mu}}_t(\check{y})}{\check{\bar{\rho}}(\check{y})}, \end{equation}
where the approximation is due to the assumption that 
$(y+1)\ll 1$ (recall that the wall is at 
$y=-1$). The other two mean-flow coefficients in (5.7) are rescaled as
\begin{align} \left. \begin{array}{c@{}}
\displaystyle \dfrac{\bar{U}_{yy}(y)}{\bar{k}_z^2} =
\Bigg(\dfrac{\bar{T}_y}{2\bar{k}_z \bar{T}} -
\dfrac{1}{\bar{k}_z (y+1)}\Bigg)
\dfrac{\bar{U}_y(y)}{\bar{k}_z} = K_1 K_2 \Bigg[
\check{\bar U}_{\check{y}\check{y}}(\check{y}) +
\Bigg(1-\dfrac{1}{K_2}\Bigg)
\dfrac{\check{\bar{U}}_{\check{y}}(\check{y})}{\check{y}}
\Bigg], \\ \displaystyle
\dfrac{\tilde{\mu}_t(y)}{\bar{\rho}(y)} \approx K_1 K_2
\Bigg[
\dfrac{\check{\tilde{\mu}}_{t,\check{y}}(\check{y})}{\check{\bar{\rho}}}
+ \Bigg(1-\dfrac{1}{K_2}\Bigg)
\dfrac{\check{\tilde{\mu}}_t(\check{y})}{\check{y}\,
\check{\bar{\rho}}(\check{y})} \Bigg]. \end{array}
\right\}
\end{align}
Therefore, if 
$K_1=K_2=1$, then 
${\bar {U}_{yy}}/{\bar {k}_z^2}$ and 
${\tilde {\mu }_t}/{\bar {\rho }}$ are independent of 
$\bar {k}_z$. The rescaled expressions in (5.8), (5.10), (5.12) and (5.13) are all explicitly irrelevant of 
$\bar {k}_z$, and dependent on 
$K_1$ and 
$K_2$ instead. The remaining mean-flow coefficient in (5.7), 
$\bar {k}_z \bar {\mu }/\bar {\rho }$, is related only to 
$\bar {T}$ and cannot be effectively made independent of 
$\bar {k}_z$. The other four equations can be analysed similarly. As a result, all the mean-flow coefficients can be categorized into two groups: (i) the ones that can be rescaled using 
$K_1$ and 
$K_2$ in the logarithmic region, including 
$\bar {U}_y$, 
$\bar {U}_{yy}$, 
$\bar {T}_y/\bar {T}$, 
$\tilde {\mu }_t/\bar {\rho }$ and 
$\tilde {\mu }_{t,y}/\bar {\rho }$ (see (5.3)); (ii) the ones related only to 
$\bar {T}$ that cannot be made explicitly irrelevant to 
$\bar {k}_z$, including 
$\bar {\mu }/\bar {\rho }$ and 
$\bar {T}/{{Ma}}_b^2$.
In summary, three factors are primarily responsible for the potential failure of the 
$k_z^{-2}$ and 
$k_z^{-1}$ scalings. The first is related to the logarithmic region. If a large portion of the shape function is out of the logarithmic region, then (5.7) and equations following all fail. This is why the scaling applies only in the mid-
$\lambda _z$ range and is consistent with the assumption of the AEM. The second factor is the substantial variations of 
$K_1$ and 
$K_2$ at different wall-normal heights. The third factor is the enormous contribution from the terms related to 
$\bar {\mu }/\bar {\rho }$ and 
$\bar {T}/{{Ma}}_b^2$. The latter two are absent in the incompressible case, both related directly to the distribution of mean temperature.
To see the variations of 
$K_1$ and 
$K_2$, the distributions of the functions 
$\sqrt {\bar {T}}$ and 
$\sqrt {\bar {T}}\,C_T$ (see (5.8) and (5.10)) are plotted in figure 17 for the cases of different 
${{Ma}}_b$ and 
${{Re}}_\tau =6000$. As can be seen, with the rise of 
${{Ma}}_b$, the two functions, and thus 
$K_1$ and 
$K_2$ at fixed 
$y$ and 
$\check {y}$, experience much larger variations. Meanwhile, the other coefficients related to 
$\bar {T}$ also change rapidly with 
$y^+$. This explains the trend observed in figures 11(c) and 15 that the 
$k^{-2}$ and 
$k^{-1}$ scaling for 
$R_{max}$ and 
$V$ have larger deviations at higher 
${{Ma}}_b$. In short, for the high-
${{Ma}}_b$ case of the channel flow, the sizeable mean temperature gradient near the wall destroys the 
$k_z$ power scaling of the response in the mid-
$\lambda _z$ range.
Wall-normal distributions of the two mean-flow-related coefficients corresponding to (a) 
$K_1$ and (b) 
$K_2$, for the cases of different 
${{Ma}}_b$ (
${{Re}}_\tau =6000$).
The above analysis should also have implications for the geometrical self-similarity of the response as discussed before. For example, 
$\sqrt {\bar {T}}$ and 
$\sqrt {\bar {T}}\,C_T$ in figure 17 change more slowly in the logarithmic region, especially for the former, which benefits the self-similarity of shape functions. However, the operation of the matrix inversion in (3.9a,b) obscures further analytical derivation on the perturbation equations and the transfer matrix. In fact, the results in figure 14 indicate the insensitivity of the self-similarity of the energy norm to 
${{Ma}}_b$.
6. Summary
In this work, the linear responses of turbulent mean flow to both harmonic and stochastic forcing are investigated for supersonic channel flow. Well-established universal relations are utilized to obtain efficiently the mean profiles with a large parameter space, with 
${{Ma}}_b$ up to 5 and 
${{Re}}_\tau$ up to 
$10^{4}$. As a result, a systematic parameter study is feasible. Comparisons with the DNS data demonstrate the reliability of the present solver in both calculating the mean flow and serving as the base flow for the linear response analysis.
Both the optimal harmonic and stochastic forcing are considered. The response of 
$k_x=0$ has the largest 
$R_{max}$ and 
$V$. The pre-multiplied energy growth 
$k_z^2 R_{max}$ exhibits the classic bi-modal structure in terms of 
$k_z$. The two peaks correspond to the large-scale motion in the outer layer, and the near-wall small-scale motion in the buffer layer, respectively. For 
$k_z V$, however, the inner peak disappears when 
${{Ma}}_b>0.5$; 
$k_z V$ increases monotonically with 
$k_z$. A further decomposition of the forcing shows that 
$k_z V_{\boldsymbol {u}}$, resulting from the kinematic part of the forcing, is quite similar to the incompressible counterpart, exhibiting the inner peak. However, the thermodynamic part, 
$k_z V_{\rho T}$, increases quickly at a larger 
$k_z$, and the resulting motions are much less coherent than the former. Specifically, the contribution of the most energetic mode of 
$V_{\rho T}$ to the energy growth is less than 3 % at 
$\lambda _z^+<100$ (
${{Ma}}_b=2$).
The spanwise wavenumber of the outer peak is 
$\lambda _{z,o}\approx 3.6h$, nearly unchanged at different 
${{Ma}}_b$ and 
${{Re}}_\tau$. The motion takes the form that the input is the streamwise vortex, and the output is the streamwise streak of velocity and temperature. The temperature streak penetrates deeper towards the wall than the streamwise velocity. The spanwise wavenumber of the inner peak 
$\lambda _{z,i}^+$ is nearly independent of 
${{Re}}_\tau$, but increases quickly with 
${{Ma}}_b$, due mainly to the wall-cooling effects in terms of the smaller ratio 
$T_w/\bar {T}_c$. From 
${{Ma}}_b=0.1$ to 5, 
$\lambda _{z,i}^+$ is lifted by over three times. If using semi-local units, 
$\lambda _{z,i}^*$ is approximately unchanged with 
${{Ma}}_b$. Wall-cooling effects enhance the coherence of the flow for the stochastic response of 
$\lambda _{z,m}^+>90$ and the optimal harmonic response of 
$\lambda _{z,m}^+>130$, consistent with the DNS data. On the contrary, it is found that the coherence of the response with 
$\lambda _z^+<\lambda _{z,m}^+$ decreases with 
${{Ma}}_b$. This is due to the strengthening of the less-coherent motion related to 
$V_{\rho T}$. As a result, the cospectrum of 
$\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ has a hot zone attached to the wall within the viscous layer.
The response modes with 
$\lambda _z$ between the inner and outer peaks exhibit a geometrical self-similarity to a first approximation, with 
${{Ma}}_b$ from 0.1 to 5. The normalized energy norm is approximately collapsed using the wall-normal coordinate 
$y^+/\lambda ^+_z$. The self-similarity is related to the logarithmic region, relatively insensitive to 
${{Ma}}_b$ and thus the compressibility effects. Finally, a theoretical analysis is conducted on the perturbation equations. The behaviour of the mean-flow coefficients is analysed within the logarithmic region. More importantly, two crucial parameters are summarized related to the dependence of the linear operator on 
$k_z$. The two parameters are related directly to the mean temperature, experiencing stronger variations in the wall-normal direction with 
${{Ma}}_b$ increase, which explains the trend that the 
$k_z^{-2}$ and 
$k_z^{-1}$ scalings of 
$R_{max}$ and 
$V$ in the mid-
$k_z$ range have increasingly large deviations with 
${{Ma}}_b$ rise.
Compared with the DNS data, the present results are more qualitative than quantitative. Some interpretations of the differences have been provided earlier in the text. Future attention will be paid to the improvement of the modelling of nonlinear forcing term, and the development of a prediction model on the energy spectra and space–time properties of the fluctuations.
Acknowledgements
The authors are grateful to K.P. Griffin for fruitful discussions on compressible velocity transformations.
Funding
This research was supported by the Center for Ocean Research in Hong Kong and Macau, a joint research centre between QNLM and HKUST. L.F. acknowledges the fund from the Research Grants Council (RGC) of the Government of Hong Kong Special Administrative Region (HKSAR) with RGC/ECS Project (no. 26200222) and the fund from the Project of Hetao Shenzhen-Hong Kong Science and Technology Innovation Cooperation Zone (no. HZQB-KCZYB-2020083).
Declaration of interests
The authors report no conflict of interests.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A. Fitting parameters and the ODE mean-flow solver
This appendix provides the details of fitting (2.7) and calculating the turbulent mean flow. For the fitting of the incompressible velocity profile, the DNS database is from Lee & Moser (Reference Lee and Moser2015) (downloaded at https://turbulence.oden.utexas.edu) with 
${{Re}}_\tau$ from 180 to 5200. The 
$\kappa _c$ and 
$A_{\nu _t}$ values are looped to minimize the global relative error 
$\epsilon (\bar {U})$ defined as
\begin{equation} \epsilon(\bar{U}) = \frac{\displaystyle\int\nolimits_{{-}1}^{0} {| \bar{U}_{fit} - \bar{U}_{DNS} | \, {\mathrm{d}} y}}{\displaystyle\int\nolimits_{{-}1}^{0} {\bar{U}_{DNS} \, {\mathrm{d}} y}}. \end{equation}
The results, with three significant digits, are given in table 1. Here, 
$\epsilon _{sgl}$ is the locally optimal error considering only a single profile, and 
$\epsilon _{glb}$ is the error using 
$\kappa _c={0.418}$ and 
$A_{\nu _t}=25.1$, which are optimized considering all five profiles with equal weights. Furthermore, a trend is observed that (2.7) predicts better the mean profile at a higher 
${{Re}}_\tau$.
Optimal fitting parameters for (2.7) and the relative error of the turbulent mean streamwise velocity for the incompressible channel flow.
The compressible turbulent mean flow can be obtained through an iterative method (Griffin et al. Reference Bai, Griffin and Fu2022; Song, Zhang & Xia Reference Song, Zhang and Xia2023). Here, we use four universal relations to derive an ODE for the mean flow, as summarized in § 2.2. In addition to (2.7), the other three relations are specified below.
(i) The total-stress-based transformation proposed by Griffin et al. (Reference Griffin, Fu and Moin2021b) well collapses nearly all the available profiles in the ranges 
${{Re}}_\tau ^*\in [200,4000]$ and 
${{Ma}}\in [0,15]$. The expression is written as
\begin{equation} \frac{\partial \bar{U}_{GFM}^+}{\partial y^*} = \left.{\dfrac{1}{\bar{\mu}^+}\,\dfrac{\partial \bar{U}^+}{\partial y^*}}\right/\left({1 + \dfrac{1}{\bar{\mu}^+}\,\dfrac{\partial \bar{U}^+}{\partial y^*} - \bar{\mu}^+\,\dfrac{\partial \bar{U}^+}{\partial y^*}}\right). \end{equation}Also, the transformation proposed by Trettel & Larsson (Reference Trettel and Larsson2016) is proven to have good performance in channel flows, written as
\begin{equation} \frac{\partial \bar{U}_{TL}^+}{\partial y^*} = \bar{\mu}^+\,\dfrac{\partial \bar{U}^+}{\partial y^+}. \end{equation}
The transformed velocities 
$\bar {U}_{GFM}^+(y^*)$ and 
$\bar {U}_{TL}^+(y^*)$ are expected to match 
$\bar {U}_{inc}^+$ in (2.7). Note that the transformation is designed for the flow within and below the logarithmic layer, and hence may deviate in the outer region. This will be checked later through the comparison with the DNS data.
(ii) For the temperature profile, Duan & Martín (Reference Duan and Martín2011) summarized the DNS data from low- to high-enthalpy flows, and provided the following improved relation for a calorically perfect gas:
\begin{equation} \left. \begin{array}{l@{}} \displaystyle \dfrac{\bar{T}}{T_\delta} = \dfrac{T_w}{T_\delta} + \dfrac{T_{aw} - T_w}{T_\delta}\,f_u \left(\dfrac{\bar{U}}{U_\delta}\right) + \dfrac{T_\delta - T_{aw}}{T_\delta} \left(\dfrac{\bar{U}}{U_\delta}\right)^2,\\ f_u (x) = (1-C_T) x^2 + C_T x, \quad C_T = 0.8259, \end{array} \right\} \end{equation}
where 
$T_\delta$ and 
$U_\delta$ are the boundary-layer edge values, and 
$T_{aw}=T_\delta [1+r_c (\gamma -1)/2\, {{Ma}}_\delta ^2]$ is the adiabatic wall temperature with 
$r_c=0.9$ as the recovery factor. The original data are based on the Favre average, while the Reynolds-averaged quantity is used here. The difference between the two averages is generally small at moderate 
${{Ma}}$ (Huang et al. Reference Huang, Coleman and Bradshaw1995), and the accuracy will be examined later. Equation (A3) also works for channel flows, where 
$T_\delta$, 
$U_\delta$ and 
${{Ma}}_\delta$ are replaced by 
$\bar {T}_c$, 
$\bar {U}_c$ and 
${{Ma}}_c$.
(iii) In terms of the outer boundary condition, Song et al. (Reference Song, Zhang, Liu and Xia2022) proposed an empirical scaling to connect 
$\bar {T}_c$ and 
$\bar {U}_c$:
\begin{equation} \frac{\bar{T}_c}{T_w} =
1 + C_{{Pr}}\,\frac{\bar{U}_c}{U_b}\,
\frac{\gamma-1}{2}\,{{Ma}}_b^2, \quad
C_{{Pr}} = 1.034 {{Pr}}.
\end{equation}
Most of the relative errors of computed 
$\bar {T}_c$ were shown to be below 1.5 % for the cases with 
${{Re}}_b$ and 
${{Ma}}_b$ ranging from 3000 to 34 000 and 0.5 to 4.0, respectively.
Based on (A2b), the equation to be solved is
\begin{equation} \frac{\bar{\mu}}{u_\tau}\,\dfrac{\partial \bar{U}}{\partial y^+} = \frac{\partial \bar{U}_{TL}^+}{\partial y^*} = \frac{{\mathrm{d}} \bar{U}_{inc}^+(y^*;{{Re}}_\tau^*)}{{\mathrm{d}} y^*}. \end{equation}
Here, 
$u_\tau$ is also an unknown. The version using (A2a) can be derived similarly. The combination of (A3) and (A4a,b) results in the following quadratic relation between 
$\bar {T}$ and 
$\bar {U}$:
\begin{equation} \frac{\bar{T}}{T_w} = 1 + C_1\,\frac{\bar{U}}{U_b} \left[(C_2 - C_3)\, \frac{\bar{U}}{U_b} + C_3\,\dfrac{\bar{U}_c}{U_b}\right], \quad \left\{ \begin{array}{l@{}} C_1 = \dfrac{\gamma-1}{2}\,{{Ma}}_b^2,\\ C_2 = C_{{Pr}}\,\dfrac{U_b}{\bar{U}_c}, \\ C_3 = C_T \left( C_2 + r_c \right). \end{array} \right. \end{equation}
Note that 
$U_b=T_w=1$ under the current definition of non-dimensional quantities (see (2.2)). The pressure is assumed to be a constant along the wall-normal direction, so the mean density is expressed as 
$\bar {\rho } = {\rho _w T_w}/{\bar {T}}$. The two unknowns 
$u_\tau$ and 
$\rho _w$ are solved using the two constraints in (2.4). Therefore, the ODE set is
\begin{equation} \left. \begin{array}{l@{}} \displaystyle F_1 (\bar{U},u_\tau) = \bar{\mu}\, \dfrac{\partial \bar{U}}{\partial y^+} - u_\tau\,\dfrac{{\mathrm{d}} \bar{U}_{inc}^+(y^*;{{Re}}_\tau^*)}{{\mathrm{d}} y^*} = 0, \\ \displaystyle F_2 (\bar{U}) = \int_{{-}1}^{0} {\dfrac{\bar{U} - 1}{\bar{T}} {\mathrm{d}} y} = 0. \end{array} \right\} \end{equation}
The only boundary condition required is 
$\bar {U}|_w=0$. For numerical solutions, the Chebyshev collocation method is adopted to cope with the derivatives and integrals. The computational domain is 
$y \in [-1,0]$, and the grid number is 
$N_y=301$. The initial 
$\bar {U}$ for iteration uses the incompressible profile. Newton's iteration is adopted for quick convergence.
For verification, the results from the present ODE solver are compared with the published DNS data. The cases with 
${{Re}}_\tau ^*>200$ are selected and the results are summarized in table 2. The relative errors of 
$\bar {T}$ and 
$\bar {\rho }$ are defined similarly to (A1). As can be seen, for the fourteen cases with 
${{Re}}_\tau ^*>300$, 
$\epsilon (\bar {U})$, 
$\epsilon (\bar {T})$ and 
$\epsilon (\bar {\rho })$ are all smaller than 1.5 %, and the average error is only 0.60 %, which demonstrates the reliability of the present ODE solver. If (A2b) is replaced by (A2a), then the average error is generally 1.5 % larger, and is located mainly in the outer region. It is noted that both (A2a) and (A2b) are not specially designed for the flow in the outer region, but (A2b) has unexpected good performance there and is thus used for later calculations.
Case parameters and the relative errors between the mean profiles from the ODE solver and the published DNS data. The abbreviations for the data source are TL from Trettel & Larsson (Reference Trettel and Larsson2016), MP from Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016), and YH from Yao & Hussain (Reference Yao and Hussain2020).
Figure 18 gives the comparison of the mean profiles for two representative cases with moderate errors, i.e. nos. (11) and (16). The excellent agreement between those from the ODE solver and the DNS data has already been shown in table 2. A visible difference in the velocity profile is near the junction of the buffer and logarithm layer, which originates from the incompressible curve fit of (2.7) and diminishes with the rise of 
${{Re}}_\tau$. The differences in the density and temperature profiles are mainly because of the slight difference in 
$\bar {T}_c$ and a slightly varying pressure along the wall-normal direction.
Mean profiles of the (a,c) streamwise velocity and (b,d) temperature for cases (a,b) no. (11) from Trettel & Larsson (Reference Trettel and Larsson2016), and (c,d) no. (16) from Yao & Hussain (Reference Yao and Hussain2020), as in table 2.
Appendix B. Verification of the response solver
The present linear response solver is verified through comparison with published results. The first case is from Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013) for an incompressible turbulent channel flow with 
${{Re}}_\tau =2003$. Figure 19 provides the comparison of the first twenty singular values of the transfer matrix 
$\boldsymbol {H}$ (see (3.9a,b)) for the harmonic response of 
$\lambda _x^+=700$, 
$\lambda _z^+=100$ and 
$\omega /k_x=10$. As can be seen, excellent agreement is obtained. The second case is from Hwang & Cossu (Reference Hwang and Cossu2010b) for an incompressible turbulent channel flow with 
${{Re}}_\tau =10^{4}$. The pre-multiplied energy amplification of the stochastic response is calculated for the mode of various 
$k_x$ and 
$k_z$. The agreement with the published results is good as well. Note that figure 19(b) is the incompressible counterpart to the results in figures 5 and 15. The third case is from Dawson & McKeon (Reference Dawson and McKeon2020) for a compressible flat-plate turbulent boundary layer with 
${{Ma}}_\infty =2$ and 
${{Re}}_\tau =900$, with the mean flow from Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011). The component shape functions of the optimal response mode (
$k_x\delta _{99}={\rm \pi} /2$, 
$k_z\delta _{99}=2{\rm \pi}$ and 
$\omega /k_x=0.8U_\infty$, with 
$\delta _{99}$ the nominal thickness) are demonstrated in figure 19(c), which is again in close agreement with the reference data. The above comparisons demonstrate the reliability of the present solver in calculating the linear responses to both harmonic and stochastic forcing.
(a) The first 20 singular values of the transfer matrix for the harmonic response and (b) the pre-multiplied energy amplification of the stochastic response of various 
$k_x$ and 
$k_z$ for the incompressible turbulent channel flow cases. (c) The shape functions of the optimal response mode for the compressible turbulent boundary layer case. The reference data are from Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013), Hwang & Cossu (Reference Hwang and Cossu2010b) and Dawson & McKeon (Reference Dawson and McKeon2020), respectively.
Appendix C. Matrix coefficients and the perturbation equation
Non-zero elements of the matrix coefficients in (3.1) are listed below. Note that 
$\tilde {\mu }_t=\mu _t+\bar {\mu }$ and 
$\tilde {\kappa }_t=\kappa _t+\bar {\kappa }$, and the subscripts 
$y$ and 
$T$ denote the corresponding partial derivatives.
Matrix 
$F$:
\begin{equation} \boldsymbol{\mathsf{F}}_{1,1} = 1, \quad \boldsymbol{\mathsf{F}}_{2,2} = \bar{\rho}, \quad \boldsymbol{\mathsf{F}}_{3,3} = \bar{\rho}, \quad \boldsymbol{\mathsf{F}}_{4,4} = \bar{\rho}, \quad \boldsymbol{\mathsf{F}}_{5,5} = \bar{\rho} c_v, \quad \boldsymbol{\mathsf{F}}_{5,1} ={-}{{Ec}}\,R \bar{T}. \end{equation} Matrix 
$A$:
\begin{equation} \left.
\begin{array}{l@{}} \displaystyle
\boldsymbol{\mathsf{A}}_{1,1} = \bar{U}, \quad
\boldsymbol{\mathsf{A}}_{1,2} = \bar{\rho}, \quad
\boldsymbol{\mathsf{A}}_{2,1} = R \bar{T}, \quad
\boldsymbol{\mathsf{A}}_{2,2} = \bar{\rho} \bar{U},
\quad \boldsymbol{\mathsf{A}}_{2,3}
={-}\dfrac{\tilde{\mu}_{t,y}}{{{Re}}_b},
\\ \displaystyle
\boldsymbol{\mathsf{A}}_{2,5} = \bar{\rho} R, \quad
\boldsymbol{\mathsf{A}}_{3,2} =
\dfrac{2}{3}\,\dfrac{\tilde{\mu}_{t,y}}{{{Re}}_b},
\quad \boldsymbol{\mathsf{A}}_{3,3} = \bar{\rho}
\bar{U}, \quad \boldsymbol{\mathsf{A}}_{3,5}
={-}\dfrac{\bar{\mu}_T
\bar{U}_y}{{{Re}}_b},\\
\displaystyle \boldsymbol{\mathsf{A}}_{4,4} =
\bar{\rho} \bar{U}, \quad
\boldsymbol{\mathsf{A}}_{5,1}
={-}{{Ec}}\,R \bar{T} \bar{U}, \quad
\boldsymbol{\mathsf{A}}_{5,3} ={-}\dfrac{2\,{{Ec}}\,\tilde{\mu}_t
\bar{U}_y}{{{Re}}_b}, \quad
\boldsymbol{\mathsf{A}}_{5,5} = \bar{\rho} c_v
\bar{U}. \end{array} \right\}
\end{equation}
Matrix 
$B$:
\begin{equation} \left. \begin{array}{l@{}} \displaystyle \boldsymbol{\mathsf{B}}_{1,3} = \bar{\rho}, \quad \boldsymbol{\mathsf{B}}_{2,2} ={-}\dfrac{\tilde{\mu}_{t,y}}{{{Re}}_b}, \quad \boldsymbol{\mathsf{B}}_{2,5} ={-} \dfrac{\bar{\mu}_T \bar{U}_y}{{{Re}}_b}, \quad \boldsymbol{\mathsf{B}}_{3,1} = R \bar{T}, \quad \boldsymbol{\mathsf{B}}_{3,3} ={-}\dfrac{4}{3}\,\dfrac{\tilde{\mu}_{t,y}}{{{Re}}_b},\\ \displaystyle \boldsymbol{\mathsf{B}}_{3,5} = \bar{\rho} R, \quad \boldsymbol{\mathsf{B}}_{4,4} ={-}\dfrac{\tilde{\mu}_{t,y}}{{{Re}}_b}, \quad \boldsymbol{\mathsf{B}}_{5,2} ={-}\dfrac{2\,{{Ec}}\,\tilde{\mu}_t \bar{U}_y}{{{Re}}_b}, \quad \boldsymbol{\mathsf{B}}_{5,5} ={-}\dfrac{\tilde{\kappa}_{t,y} + \bar{\kappa}_T \bar{T}_y}{{{Re}}_b}. \end{array} \right\} \end{equation} Matrix 
$C$:
\begin{equation} \boldsymbol{\mathsf{C}}_{1,4} = \bar{\rho}, \quad \boldsymbol{\mathsf{C}}_{3,4} = \frac{2}{3}\,\frac{\tilde{\mu}_{t,y}}{{{Re}}_b}, \quad \boldsymbol{\mathsf{C}}_{4,1} = R \bar{T}, \quad \boldsymbol{\mathsf{C}}_{4,3} ={-}\frac{\tilde{\mu}_{t,y}}{{{Re}}_b}, \quad \boldsymbol{\mathsf{C}}_{4,5} = \bar{\rho} R. \end{equation} Matrix 
$D$:
\begin{equation} \left. \begin{array}{l@{}} \displaystyle \boldsymbol{\mathsf{D}}_{1,3} \!=\! \bar{\rho}_y, \quad \boldsymbol{\mathsf{D}}_{2,3} \!=\! \bar{\rho} \bar{U}_y, \quad \boldsymbol{\mathsf{D}}_{2,5} ={-}\dfrac{\bar{\mu}_T \bar{U}_{yy} \!+\! \bar{\mu}_{yT} \bar{U}_y}{{{Re}}_b}, \quad \boldsymbol{\mathsf{D}}_{3,1} = R \bar{T}_y, \quad \boldsymbol{\mathsf{D}}_{3,5} = \bar{\rho}_y R,\\ \displaystyle \boldsymbol{\mathsf{D}}_{5,3} = \bar{\rho} c_v \bar{T}_y - {{Ec}}\,\bar{T} R \bar{\rho}_y, \quad \boldsymbol{\mathsf{D}}_{5,5} ={-}\dfrac{{{Ec}}\,\bar{\mu}_T \bar{U}_y^2+ \bar{\kappa}_{yT} \bar{T}_y + \bar{\kappa}_T \bar{T}_{yy}}{{{Re}}_b}. \end{array} \right\} \end{equation} Matrix 
$H$:
\begin{equation} \left. \begin{array}{l@{}} \displaystyle \boldsymbol{\mathsf{H}}_{xx,2,2} = \dfrac{4}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{yy,2,2} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zz,2,2} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \\ \displaystyle\boldsymbol{\mathsf{H}}_{xy,2,3} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zx,2,4} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b},\\ \displaystyle \boldsymbol{\mathsf{H}}_{xx,3,3} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{yy,3,3} = \dfrac{4}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zz,3,3} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \\ \displaystyle\boldsymbol{\mathsf{H}}_{xy,3,2} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{yz,3,4} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b},\\ \displaystyle \boldsymbol{\mathsf{H}}_{xx,4,4} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{yy,4,4} = \dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zz,4,4} = \dfrac{4}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \\ \displaystyle\boldsymbol{\mathsf{H}}_{yz,4,3} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zx,4,2} = \dfrac{1}{3}\,\dfrac{\tilde{\mu}_t}{{{Re}}_b},\\ \displaystyle \boldsymbol{\mathsf{H}}_{xx,5,5} = \dfrac{\tilde{\kappa}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{yy,5,5} = \dfrac{\tilde{\kappa}_t}{{{Re}}_b}, \quad \boldsymbol{\mathsf{H}}_{zz,5,5} = \dfrac{\tilde{\kappa}_t}{{{Re}}_b}. \end{array} \right\} \end{equation}In addition to the two equations in (5.1), the other three perturbation equations are written as
\begin{align} &\frac{\partial
\hat{v}}{\partial t} + \frac{\bar{T}}{\gamma\,{{Ma}}_b^2}
\left(
\frac{\bar{T}_y}{\bar{T}}\,\frac{\hat{\rho}}{\bar{\rho}} +
\frac{\hat{\rho}_y}{\bar{\rho}} -
\frac{\bar{T}_y}{\bar{T}}\, \frac{\hat{T}}{\bar{T}} +
\frac{\hat{T}_y}{\bar{T}} \right)\nonumber\\ &\quad -
\frac{1}{{{Re}}_b} \left[
\frac{\tilde{\mu}_{t,y}}{\bar{\rho}}
\left(\frac{4}{3}\,\hat{v}_y - \frac{2
{\mathrm{i}} k_z}{3}\,\hat{w} \right) + \frac{\tilde{\mu}_t}{\bar{\rho}} \left(\frac{4}{3}\,
\hat{v}_{yy} - k_z^2 \hat{v} +
\frac{{\mathrm{i}} k_z}{3}\,\hat{w}_y\right)
\right] = \hat{f}_v,
\end{align}
\begin{align} &\frac{\partial
\hat{w}}{\partial t} + {\mathrm{i}}
k_z\,\frac{\bar{T}}{\gamma\, {{Ma}}_b^2}
\left( \frac{\hat{\rho}}{\bar{\rho}} +
\frac{\hat{T}}{\bar{T}} \right)\nonumber\\ &\quad -
\frac{1}{{{Re}}_b} \left[
\frac{\tilde{\mu}_{t,y}}{\bar{\rho}} \left(\hat{w}_y +
{\mathrm{i}} k_z \hat{v} \right) +
\frac{\tilde{\mu}_t}{\bar{\rho}} \left(\hat{w}_{yy} -
\frac{4}{3}\,k_z^2 \hat{w} + \frac{{\mathrm{i}}
k_z}{3}\,\hat{v}_y\right) \right] = \hat{f}_w,
\end{align}
\begin{align} &\frac{\partial}{\partial t} \frac{\hat{T}}{\bar{T}} + \frac{\bar{T}_y}{\bar{T}}\,\hat{v} + \left(\gamma-1\right)\left(\hat{v}_y + {\mathrm{i}} k_z \hat{w}\right) - \frac{\gamma}{{{Re}}_b} \left[ \frac{\kappa_{t,y}}{\bar{\rho}}\,\frac{\hat{T}_y}{\bar{T}} + \frac{\kappa_t}{\bar{\rho}} \left(\frac{\hat{T}_{yy}}{\bar{T}} - k_z^2\, \frac{\hat{T}}{\bar{T}}\right) \right]\nonumber\\ &\quad - \frac{\gamma - 1}{{{Re}}_b}\,\frac{\gamma\,{{Ma}}_b^2}{\bar{T}}\, \bar{U}_y \left(n \bar{U}_y\,\frac{\bar{\mu}}{\bar{\rho}}\, \frac{\hat{T}}{\bar{T}} + 2\,\frac{\tilde{\mu}_t}{\bar{\rho}}\,\hat{u}_y \right)\nonumber\\ &\quad - \frac{n\gamma}{{{Re}}_b}\,\frac{\bar{\kappa}}{\bar{\rho}} \left[ \left(\frac{\bar{T}_{yy}}{\bar{T}} + (n-1) \left(\frac{\bar{T}_y}{\bar{T}}\right)^2 \right) \frac{\hat{T}}{\bar{T}} + \frac{\bar{T}_y}{\bar{T}}\,\frac{\hat{T}_y}{\bar{T}} \right] = \frac{\hat{f}_T}{\bar{T}}. \end{align}
Open access
