# Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers

Riwan Hammachi, Estelle Piot, Hugues Deniau, Jean-Philippe Brazier, Guillaume Daviller, José I. Cardesa

► To cite this version:

Riwan Hammachi, Estelle Piot, Hugues Deniau, Jean-Philippe Brazier, Guillaume Daviller, et al. Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers. Journal of Physics: Conference Series, 2024, 2753, pp.012012. 10.1088/1742-6596/2753/1/012012. [hal-04596746]

HAL Id: hal-04596746

https://hal.science/hal-04596746

Submitted on 31 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers

To cite this article: Riwan Hammachi et al 2024 J. Phys.: Conf. Ser. 2753 012012

View the article online for updates and enhancements.

You may also like

- The Intrinsic Properties of Multiwavelength Energy Spectra for Fermi Teraelectronvolt Blazars

R. X. Zhou, Y. G. Zheng, K. R. Zhu et al.

- FIRST-EPOCH VLBA IMAGING OF 20 NEW TeV BLAZARS

B. Glenn Piner and Philip G. Edwards

- Multi-epoch VLBA Imaging of 20 New TeV Blazars: Apparent Jet Speeds

B. Glenn Piner and Philip G. Edwards

<figure><img src="image_1.png" /></figure>

The
Electrochemical
Society

Advancing solid state &
electrochemical science & technology

<figure><img src="image_2.png" /></figure>

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_3.png" /></figure>
PAPER • OPEN ACCESS

# Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers

To cite this article: Riwan Hammachi et al 2024 J. Phys.: Conf. Ser. 2753 012012

View the article online for updates and enhancements.

You may also like

- The Intrinsic Properties of Multiwavelength Energy Spectra for Fermi Teraelectronvolt Blazars
R. X. Zhou, Y. G. Zheng, K. R. Zhu et al.

- FIRST-EPOCH VLBA IMAGING OF 20 NEW TeV BLAZARS
B. Glenn Piner and Philip G. Edwards

- Multi-epoch VLBA Imaging of 20 New TeV Blazars: Apparent Jet Speeds
B. Glenn Piner and Philip G. Edwards

<figure><img src="image_2.png" /></figure>

The
Electrochemical
Society

Advancing solid state &
electrochemical science & technology

<figure><img src="image_3.png" /></figure>

DISCOVER
how sustainability
intersects with
electrochemistry & solid
state science research

<figure><img src="image_4.png" /></figure>
# Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers

Riwan Hammachi<fnref n="1" />,<fnref n="2" />, Estelle Piot<fnref n="1" />, Hugues Deniau<fnref n="1" />, Jean-Philippe Brazier<fnref n="1" />, Guillaume Daviller<fnref n="2" />, José I. Cardesa<fnref n="1" />

E-mail: riwan.hammachi@onera.fr

**Abstract.** The natural transition of hypersonic boundary layers (HBLs) is often expressed in terms of discrete modes and their linear stability. A frequent interpretation revolves around fast and slow acoustic modes interacting in the vicinity of the vortical/entropic branches of the continuous spectrum found from stability analyses. Yet several transition scenarios are contingent upon factors such as the spectral content of the free-stream disturbances, or the interactions between the discrete modes within the boundary layer and the free-stream disturbances near the leading edge which can be decomposed into vortical, acoustic and entropic nature based on the fluid-thermodynamic (FT) components. Yet the interpretations of linear stability applied to discrete modes can lead to semantic conflicts with the terminology of FT components. To clarify the current description of the processes involved, this study chooses an approach aimed at characterizing the dynamics of the second Mack mode in transitional HBLs through coherent structure tracking. The method involves decomposing the flow perturbations into acoustic, vortical and entropic content, and following their associated coherent structures over time. For this purpose, direct numerical simulations are carried out to investigate the dynamics of the second Mack mode instability in two-dimensional HBLs, considering a flow at Mach 6 over a cooled and an insulated wall. It is found that vortical structures coexist at different heights along the wall surface, forming alternating sign doublets around the critical layer and above the relative sonic line. These structures are found to merge in the region of maximum second Mack mode instability.

## 1. Introduction

The pursuit of hypersonic-capable vehicles represents a complex technological challenge and a profound scientific endeavor. Among the myriad hurdles, the comprehension and mastery of hypersonic aerodynamics stands as paramount. At these high velocities, the transition to turbulence in the near-wall flow introduces a substantial increase in surface skin friction and significant wall heating. Consequently, unraveling the transition mechanisms within the hypersonic boundary layer (HBL) takes on pivotal importance in shaping the design of high-performance hypersonic vehicles. According to our current understanding of the physics governing this phenomenon, transition is a complex process that varies depending on numerous parameters such as environmental conditions, mean flow properties, and vehicle geometry. In the vicinity of hypersonic vehicles, the boundary layer is susceptible to various disturbances that can impact its dynamics, ultimately leading to the onset of turbulence [1]. In light of the consensus that atmospheric conditions during hypersonic flight correspond to a low-disturbance

<figure><img src="image_4.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ONERA/DMPE, Université de Toulouse, Toulouse, F-31055, France</footnote>

<footnote n="2">²CERFACS, 42 avenue Gaspard Coriolis, Toulouse, 31057, France</footnote>
# Coherent structure tracking of the second Mack mode in transitional hypersonic boundary layers

Riwan Hammachi¹,², Estelle Piot¹, Hugues Deniau¹, Jean-Philippe Brazier¹, Guillaume Daviller², José I. Cardesa¹

¹ONERA/DMPE, Université de Toulouse, Toulouse, F-31055, France

²CERFACS, 42 avenue Gaspard Coriolis, Toulouse, 31057, France

E-mail: riwan.hammachi@onera.fr

**Abstract.** The natural transition of hypersonic boundary layers (HBLs) is often expressed in terms of discrete modes and their linear stability. A frequent interpretation revolves around fast and slow acoustic modes interacting in the vicinity of the vortical/entropic branches of the continuous spectrum found from stability analyses. Yet several transition scenarios are contingent upon factors such as the spectral content of the free-stream disturbances, or the interactions between the discrete modes within the boundary layer and the free-stream disturbances near the leading edge which can be decomposed into vortical, acoustic and entropic nature based on the fluid-thermodynamic (FT) components. Yet the interpretations of linear stability applied to discrete modes can lead to semantic conflicts with the terminology of FT components. To clarify the current description of the processes involved, this study chooses an approach aimed at characterizing the dynamics of the second Mack mode in transitional HBLs through coherent structure tracking. The method involves decomposing the flow perturbations into acoustic, vortical and entropic content, and following their associated coherent structures over time. For this purpose, direct numerical simulations are carried out to investigate the dynamics of the second Mack mode instability in two-dimensional HBLs, considering a flow at Mach 6 over a cooled and an insulated wall. It is found that vortical structures coexist at different heights along the wall surface, forming alternating sign doublets around the critical layer and above the relative sonic line. These structures are found to merge in the region of maximum second Mack mode instability.

## 1. Introduction

The pursuit of hypersonic-capable vehicles represents a complex technological challenge and a profound scientific endeavor. Among the myriad hurdles, the comprehension and mastery of hypersonic aerodynamics stands as paramount. At these high velocities, the transition to turbulence in the near-wall flow introduces a substantial increase in surface skin friction and significant wall heating. Consequently, unraveling the transition mechanisms within the hypersonic boundary layer (HBL) takes on pivotal importance in shaping the design of high-performance hypersonic vehicles. According to our current understanding of the physics governing this phenomenon, transition is a complex process that varies depending on numerous parameters such as environmental conditions, mean flow properties, and vehicle geometry. In the vicinity of hypersonic vehicles, the boundary layer is susceptible to various disturbances that can impact its dynamics, ultimately leading to the onset of turbulence [1]. In light of the consensus that atmospheric conditions during hypersonic flight correspond to a low-disturbance

<figure><img src="image_5.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
environment, the pathway toward turbulence follows the “natural transition” scenario which predominantly hinges on the modal growth of the relevant instabilities [2]. In the initial phase, the freetream disturbances enter the boundary layer and excite its wave modes through the receptivity process. Subsequently, unstable eigenmodes (when present) undergo exponential growth up to certain amplitudes past which secondary instabilities, nonlinear interactions and three-dimensional effects enter the picture, leading ultimately to breakdown and turbulence. It is important to note that in the hypersonic regime, multiple instability mechanisms coexist within the boundary layer. This is in contrast with two-dimensional (2D) low-speed boundary layers, which typically exhibit only the viscous Tollmien–Schlichting (T–S) wave instability [3]. A specific characteristic of a 2D hypersonic boundary layer lies in the existence of a family of inviscid unstable modes: the Mack modes. In the terminology of Mack [4], they are denoted as first, second and higher modes. In his study, Mack revealed the influence of wall temperature on their stability. Wall cooling destabilizes the second and higher mode instabilities while stabilizing the first one [5]. Moreover, an increase in the Mach number at the boundary layer edge ($M_e$) produces similar effects. In the case of insulated 2D boundary layers at hypersonic speeds, starting from about $M_e > 4$, the transition is primarily governed by second-mode waves [6]. The Mack second, third, and higher modes are convectively-amplified waveguide normal modes arising from an inviscid instability which is characterized by trapped acoustic waves [7, 8, 9]. The existence of these waves is sustained by the existence of a relative sonic line at a wall height $y_s$, where $\bar{M}(y_s) = \frac{\bar{u}(y_s)-c}{\bar{a}(y_s)} = -1$ [10, 5]. Here, $\bar{M}$, $\bar{a}$, and $c$ are the relative Mach number, the speed of sound and the disturbance propagation speed, respectively. The presence of such a line implies that relative to the flow speed, disturbances travel subsonically above it and supersonically below it. Hence the lower part of the boundary layer behaves as an acoustic waveguide, where acoustic rays are reflected off the wall and turn around near the relative sonic line towards the wall. Above this turning point at $y = y_s$, the disturbance generates a distinctive “rope-like” wave pattern, as observed by numerous researchers both experimentally and numerically [11, 12]. These structures are centered around $\bar{M}(y_c) = 0$, where $y_c$ denotes the critical layer height. A schematic view in Figure 1, adapted from [13], concisely summarizes previous authors explanations on the behavior of the second Mack mode.

<figure><img src="image_5.png" /><figcaption>Figure 1: Schematic view of the Mack second mode in a high-speed boundary layer, adapted from Knisely et al. [13]</figcaption></figure>

The term “acoustic” has been resorted to in order to classify the physical nature of the second and higher Mack modes. It is common practice in this field to sort the various boundary layer modes into categories that link them in one way or another to a type of physical process or quantity. Examples include entropy modes, vorticity modes, fast or slow acoustic modes. This habit feeds from results provided by linear stability analysis, a tool that has significantly shaped our understanding of hypersonic boundary layer transition mechanisms. The spectrum obtained from such analysis and the structure of the normal modes have led to a refinement in the terminology used to describe the unstable modes in the hypersonic boundary layer. Yet
environment, the pathway toward turbulence follows the “natural transition” scenario which predominantly hinges on the modal growth of the relevant instabilities [2]. In the initial phase, the freetream disturbances enter the boundary layer and excite its wave modes through the receptivity process. Subsequently, unstable eigenmodes (when present) undergo exponential growth up to certain amplitudes past which secondary instabilities, nonlinear interactions and three-dimensional effects enter the picture, leading ultimately to breakdown and turbulence. It is important to note that in the hypersonic regime, multiple instability mechanisms coexist within the boundary layer. This is in contrast with two-dimensional (2D) low-speed boundary layers, which typically exhibit only the viscous Tollmien–Schlichting (T–S) wave instability [3]. A specific characteristic of a 2D hypersonic boundary layer lies in the existence of a family of inviscid unstable modes: the Mack modes. In the terminology of Mack [4], they are denoted as first, second and higher modes. In his study, Mack revealed the influence of wall temperature on their stability. Wall cooling destabilizes the second and higher mode instabilities while stabilizing the first one [5]. Moreover, an increase in the Mach number at the boundary layer edge ($M_e$) produces similar effects. In the case of insulated 2D boundary layers at hypersonic speeds, starting from about $M_e > 4$, the transition is primarily governed by second-mode waves [6]. The Mack second, third, and higher modes are convectively-amplified waveguide normal modes arising from an inviscid instability which is characterized by trapped acoustic waves [7, 8, 9]. The existence of these waves is sustained by the existence of a relative sonic line at a wall height $y_s$, where $\bar{M}(y_s) = \frac{\bar{u}(y_s)-c}{\bar{a}(y_s)} = -1$ [10, 5]. Here, $\bar{M}$, $\bar{a}$, and $c$ are the relative Mach number, the speed of sound and the disturbance propagation speed, respectively. The presence of such a line implies that relative to the flow speed, disturbances travel subsonically above it and supersonically below it. Hence the lower part of the boundary layer behaves as an acoustic waveguide, where acoustic rays are reflected off the wall and turn around near the relative sonic line towards the wall. Above this turning point at $y = y_s$, the disturbance generates a distinctive “rope-like” wave pattern, as observed by numerous researchers both experimentally and numerically [11, 12]. These structures are centered around $\bar{M}(y_c) = 0$, where $y_c$ denotes the critical layer height. A schematic view in Figure 1, adapted from [13], concisely summarizes previous authors explanations on the behavior of the second Mack mode.

<figure><img src="image_6.png" /><figcaption>Figure 1: Schematic view of the Mack second mode in a high-speed boundary layer, adapted from Knisely et al. [13]</figcaption></figure>

The term “acoustic” has been resorted to in order to classify the physical nature of the second and higher Mack modes. It is common practice in this field to sort the various boundary layer modes into categories that link them in one way or another to a type of physical process or quantity. Examples include entropy modes, vorticity modes, fast or slow acoustic modes. This habit feeds from results provided by linear stability analysis, a tool that has significantly shaped our understanding of hypersonic boundary layer transition mechanisms. The spectrum obtained from such analysis and the structure of the normal modes have led to a refinement in the terminology used to describe the unstable modes in the hypersonic boundary layer. Yet
the nomenclature of Mack's first and second modes is not particularly descriptive and can be limiting, both in qualitatively describing the transition mechanisms associated with these instabilities and in interpreting the physical aspects of these phenomena. In fact, it turns out that the first and second modes are not truly distinct modes. Notably, the works of Fedorov and Khokhlov [14] have revealed the continuity of the spectrum between Mack's so-called first and second modes. Fedorov and Tumin [15] subsequently introduced a new terminology based on the behavior of truly distinct modes in the linear stability spectrum. Their new terminology aimed to provide an alternative interpretation of the physical and mathematical description of hypersonic boundary layer instabilities, introducing modes S and F1 alongside continuous modes already identified previously: entropy, vorticity, fast and slow acoustic. For sake of simplicity, we follow the notation introduced by Fedorov [9] and label mode F1 as F.

Near the leading edge, modes S and F initiate in the slow (S) and the fast (F) acoustic part of the spectrum, respectively. Their initial phase velocity close to that of acoustic waves promotes the receptivity of these two modes to acoustic disturbances in the flow. Moving downstream, the phase velocity of mode F decreases to a point where it may synchronize with the vortical and entropic modes of the continuous spectrum. At this point, mode F becomes susceptible to free-stream temperature spottiness and turbulence [16]. Further downstream, if the tendency of mode F to become slower and S to become faster is unaltered, a second synchronization point occurs. Both synchronization points impact the stability of modes S and F. Many studies employing direct numerical simulations have verified these observations [17, 18]. Figure 2 illustrates the different scenarios leading to the growth of Mack modes in relation to the receptivity mechanisms of the flow disturbances. It is important to highlight that the evolution of the S and F modes, their possible interaction and their stability are complex functions of the flow conditions. Depending on the Mach number and the wall temperature, the stability of the S and F branches can switch, or the synchronization point between S and F modes can vanish [15].

<figure><img src="image_6.png" /><figcaption>Figure 2: Scenarios of Mack second mode evolution for a high-speed boundary layer on a flat plate, adapted from Fedorov *et al.* [9]</figcaption></figure>

While linear stability analysis have shed light on the role of the boundary layer modes in the transition process, the method is far from complete when it comes to physical interpretation. Its shortcomings motivated extensions or alternative approaches. Kuehl [19] proposed a thermoacoustic analysis on the second mode using an inviscid, Lagrangian, cycle-averaged disturbance acoustic energy formulation, leading to the hypothesis that the physical origins of Mack's second mode are in the thermoacoustic Reynolds stress. Tian *et al.* [20] performed a stability analyses based on the rates of change of perturbations associated to the Mack second mode. They found that the variation in fluctuating internal energy is maintained by the advection of perturbed thermal energy near the critical layer (where the flow propagates at the velocity of the disturbance) and the fluctuation in dilatation near the wall. Additionally, streamwise velocity perturbations were attributed to the simultaneous influence of momentum transfer from the wall-normal velocity fluctuation and the streamwise gradient of the pressure perturbation near the wall. Chen *et al.* [21] conducted an analysis of the phase of each term
the nomenclature of Mack's first and second modes is not particularly descriptive and can be limiting, both in qualitatively describing the transition mechanisms associated with these instabilities and in interpreting the physical aspects of these phenomena. In fact, it turns out that the first and second modes are not truly distinct modes. Notably, the works of Fedorov and Khokhlov [14] have revealed the continuity of the spectrum between Mack's so-called first and second modes. Fedorov and Tumin [15] subsequently introduced a new terminology based on the behavior of truly distinct modes in the linear stability spectrum. Their new terminology aimed to provide an alternative interpretation of the physical and mathematical description of hypersonic boundary layer instabilities, introducing modes S and F1 alongside continuous modes already identified previously: entropy, vorticity, fast and slow acoustic. For sake of simplicity, we follow the notation introduced by Fedorov [9] and label mode F1 as F.

Near the leading edge, modes S and F initiate in the slow (S) and the fast (F) acoustic part of the spectrum, respectively. Their initial phase velocity close to that of acoustic waves promotes the receptivity of these two modes to acoustic disturbances in the flow. Moving downstream, the phase velocity of mode F decreases to a point where it may synchronize with the vortical and entropic modes of the continuous spectrum. At this point, mode F becomes susceptible to free-stream temperature spottiness and turbulence [16]. Further downstream, if the tendency of mode F to become slower and S to become faster is unaltered, a second synchronization point occurs. Both synchronization points impact the stability of modes S and F. Many studies employing direct numerical simulations have verified these observations [17, 18]. Figure 2 illustrates the different scenarios leading to the growth of Mack modes in relation to the receptivity mechanisms of the flow disturbances. It is important to highlight that the evolution of the S and F modes, their possible interaction and their stability are complex functions of the flow conditions. Depending on the Mach number and the wall temperature, the stability of the S and F branches can switch, or the synchronization point between S and F modes can vanish [15].

<figure><img src="image_7.png" /><figcaption>Figure 2: Scenarios of Mack second mode evolution for a high-speed boundary layer on a flat plate, adapted from Fedorov *et al.* [9]</figcaption></figure>

While linear stability analysis have shed light on the role of the boundary layer modes in the transition process, the method is far from complete when it comes to physical interpretation. Its shortcomings motivated extensions or alternative approaches. Kuehl [19] proposed a thermoacoustic analysis on the second mode using an inviscid, Lagrangian, cycle-averaged disturbance acoustic energy formulation, leading to the hypothesis that the physical origins of Mack's second mode are in the thermoacoustic Reynolds stress. Tian *et al.* [20] performed a stability analyses based on the rates of change of perturbations associated to the Mack second mode. They found that the variation in fluctuating internal energy is maintained by the advection of perturbed thermal energy near the critical layer (where the flow propagates at the velocity of the disturbance) and the fluctuation in dilatation near the wall. Additionally, streamwise velocity perturbations were attributed to the simultaneous influence of momentum transfer from the wall-normal velocity fluctuation and the streamwise gradient of the pressure perturbation near the wall. Chen *et al.* [21] conducted an analysis of the phase of each term
in the linear stability equation to identify the significant physical sources of boundary layer instabilities. Noteworthy among their observations is the role of wall cooling, which destabilizes the second mode by establishing a pronounced region of wall-normal internal energy transport beneath the second generalized inflection point (GIP).

Unnikrishnan et al. [22] caution against semantic misinterpretations made in the literature associated with interpreting the nature of stability modes in terms of acoustic, vortical, or entropic content through their trajectory in the eigenspectra, which prompts the questions: can these three contributions be isolated in the flow? Can the dynamics of these contributions be studied in order to elucidate their true role in the transition process? A first step towards this involves being able to carry out an energy decomposition where each component can be assigned a type of contribution: acoustic, vortical, or entropic. Such a classification for small energy disturbances in a uniform flow can be inferred from the work of Chu [23], yet its limitations have given rise to several studies on hypersonic flows which rely on the momentum potential theory (MPT) decomposition by Doak [24] in order to individually analyse the different contributions [25, 22, 26, 27, 28]. Based on Doak’s MPT, Unnikrishnan et al. [22] were able to study the evolution of fluid-thermodynamic (FT) components associated with the Mack second mode instability, revealing many interesting results. In particular, they showed that vortical components are the most prominent in the Mack second mode instability, followed by thermal and acoustic components. They also showed that the thermal component is the most significant source term for the growth of acoustic perturbations, an effect which becomes pronounced under wall cooling conditions. Their results highlight that the vortical structures manifest as a series of rope-shaped cells across the GIP, and that both the acoustic and thermal components undergo a change in structural shape downstream of the mode synchronization exhibiting a “trapped” structure. The acoustic component displays alternating monopole arrays between the wall and the critical layer, while the thermal component is confined to the region between the wall and a line straddling the GIP. The MPT-based FT decomposition was also fruitfully used by Long et al. [26] to study the “sound radiation” caused by the supersonic mode [29], and the influence of a porous coating on this instability. All these studies seem to converge on the significant role of the GIP and the critical layer in energy amplification mechanisms, as well as on the utility of the Doak decomposition in studying the unstable modes in hypersonic boundary layers.

In the present study, we report on our progress towards exploring the dynamics of the second-mode instability through FT components. The study begins with a linear stability analysis to identify the distinct linear modes within hypersonic 2D flat plate boundary layers under different wall temperature conditions. Following this, forced direct numerical simulations (DNS) are conducted, introducing linear harmonic disturbances in the boundary layer that excite the normal modes predicted by the stability analyses. We then employ the tracking method pioneered by Lozano-Durán and Jiménez [30], and apply it to a temporal sequence of the fluctuating vorticity. This preliminary step towards our final goal allows us to track coherent structures over time, providing both qualitative and quantitative insights into their temporal behaviour, their lifetimes, their geometrical characterization and their advection velocities in a hypersonic boundary layer. Finally, the Doak decomposition is applied to the DNS results in order to compute, for individual flow fields, the FT components of the perturbations. Comparisons are drawn between the components observed over insulated and cooled wall conditions.

## 2. Simulation conditions

We consider two cases of two-dimensional HBLs over a flat plate. One where the wall temperature $T_w$ is equal to the recovery temperature $T_r$ (adiabatic case), another where $T_w/T_r = 0.145$ with $T_w = 300$ K. The free-stream Mach and unit Reynolds numbers are given
in the linear stability equation to identify the significant physical sources of boundary layer instabilities. Noteworthy among their observations is the role of wall cooling, which destabilizes the second mode by establishing a pronounced region of wall-normal internal energy transport beneath the second generalized inflection point (GIP).

Unnikrishnan et al. [22] caution against semantic misinterpretations made in the literature associated with interpreting the nature of stability modes in terms of acoustic, vortical, or entropic content through their trajectory in the eigenspectra, which prompts the questions: can these three contributions be isolated in the flow? Can the dynamics of these contributions be studied in order to elucidate their true role in the transition process? A first step towards this involves being able to carry out an energy decomposition where each component can be assigned a type of contribution: acoustic, vortical, or entropic. Such a classification for small energy disturbances in a uniform flow can be inferred from the work of Chu [23], yet its limitations have given rise to several studies on hypersonic flows which rely on the momentum potential theory (MPT) decomposition by Doak [24] in order to individually analyse the different contributions [25, 22, 26, 27, 28]. Based on Doak’s MPT, Unnikrishnan et al. [22] were able to study the evolution of fluid-thermodynamic (FT) components associated with the Mack second mode instability, revealing many interesting results. In particular, they showed that vortical components are the most prominent in the Mack second mode instability, followed by thermal and acoustic components. They also showed that the thermal component is the most significant source term for the growth of acoustic perturbations, an effect which becomes pronounced under wall cooling conditions. Their results highlight that the vortical structures manifest as a series of rope-shaped cells across the GIP, and that both the acoustic and thermal components undergo a change in structural shape downstream of the mode synchronization exhibiting a “trapped” structure. The acoustic component displays alternating monopole arrays between the wall and the critical layer, while the thermal component is confined to the region between the wall and a line straddling the GIP. The MPT-based FT decomposition was also fruitfully used by Long et al. [26] to study the “sound radiation” caused by the supersonic mode [29], and the influence of a porous coating on this instability. All these studies seem to converge on the significant role of the GIP and the critical layer in energy amplification mechanisms, as well as on the utility of the Doak decomposition in studying the unstable modes in hypersonic boundary layers.

In the present study, we report on our progress towards exploring the dynamics of the second-mode instability through FT components. The study begins with a linear stability analysis to identify the distinct linear modes within hypersonic 2D flat plate boundary layers under different wall temperature conditions. Following this, forced direct numerical simulations (DNS) are conducted, introducing linear harmonic disturbances in the boundary layer that excite the normal modes predicted by the stability analyses. We then employ the tracking method pioneered by Lozano-Durán and Jiménez [30], and apply it to a temporal sequence of the fluctuating vorticity. This preliminary step towards our final goal allows us to track coherent structures over time, providing both qualitative and quantitative insights into their temporal behaviour, their lifetimes, their geometrical characterization and their advection velocities in a hypersonic boundary layer. Finally, the Doak decomposition is applied to the DNS results in order to compute, for individual flow fields, the FT components of the perturbations. Comparisons are drawn between the components observed over insulated and cooled wall conditions.

## 2. Simulation conditions

We consider two cases of two-dimensional HBLs over a flat plate. One where the wall temperature $T_w$ is equal to the recovery temperature $T_r$ (adiabatic case), another where $T_w/T_r = 0.145$ with $T_w = 300$ K. The free-stream Mach and unit Reynolds numbers are given
by:

$$M_{\infty} = \frac{U_{\infty}}{a_{\infty}} = 6, \qquad (1)$$

$$Re_{\infty} = \frac{\rho_{\infty}U_{\infty}}{\mu_{\infty}} = 2 \cdot 10^6 \text{ m}^{-1}, \qquad (2)$$

where $ρ_∞$, $U_∞$, and $μ_∞$ are the density, velocity, and dynamic viscosity of the free-stream, respectively. $a_∞ = \sqrt{γRT_∞}$ is the speed of sound based on the free-stream temperature $T_∞$. The flow is considered to be a calorically perfect gas, with specific heat ratio $γ = 1.4$ and gas constant $R = 287.05 \text{ J/(kg} \cdot \text{K)}$. The dynamic viscosity $μ$ is chosen to follow Sutherland's law, and the heat conductivity $κ$ is related to $μ$ through a constant Prandtl number $Pr = 0.72$.

<figure><img src="image_7.png" /><figcaption>Figure 3: Sketch of the physical model and computational domains.</figcaption></figure>

To initialize the forced simulations, a precursor simulation is run with the finite-volume flow solver CEDRE [31]. This first simulation computes the steady base flow over the flat plate with a second-order finite volume spatial discretization, along with a first-order implicit temporal scheme. The computational domain of the steady base flow is referred to as the “low-order” region in Figure 3. Subsequently, the steady flow field is applied as both the initial condition and as the inlet and outlet boundary conditions within a subdomain where we perform the high-fidelity simulations, denoted as the “high-order computational domain” in Figure 3. The details of the domain geometry for this study are summarized in Table 1. We used bilinear interpolation to project the low-order solution onto the high-order domain. The high-fidelity simulations are run with JAGUAR [32, 33, 34, 35], a high-order spectral difference flow solver collaboratively developed by CERFACS and ONERA. At first, JAGUAR is run without any forcing applied to it for long enough to evacuate initial transients. The latter are due to the fact that the two solvers have different truncation errors. A steady-state is then achieved with JAGUAR, which is in excellent agreement with the low-order base flow solution given by CEDRE. It should be noted that in the case of the cold flat plate, the wall temperature is set to $T_w = 300 \text{ K}$ for both CEDRE and JAGUAR. However, for the adiabatic wall, an isothermal wall boundary condition is used in JAGUAR to enforce the wall temperature to $T_w = T_r(x)$, where $T_r(x)$ is the adiabatic wall temperature obtained from CEDRE when enforcing an adiabatic boundary condition.

Once a newly converged steady solution is achieved with Jaguar, it serves as the initial steady state for perturbed unsteady simulations with the same solver. To initiate the Mach second mode instability in the boundary layer, a wall blowing-suction actuator is employed, located within the region between $x_{f1}$ and $x_{f2}$, as depicted in Figure 3. The blowing/suction velocity follows a sinusoidal space-time dependency given by:

$$v'(t) = A_f \sin \left( 2\pi \frac{x - x_{f2}}{x_{f1} - x_{f2}} \right) \sin(\omega t), \quad x_{f1} \le x \le x_{f2} \qquad (3)$$
by:

$$M_{\infty} = \frac{U_{\infty}}{a_{\infty}} = 6, \qquad (1)$$

$$Re_{\infty} = \frac{\rho_{\infty}U_{\infty}}{\mu_{\infty}} = 2 \cdot 10^{6} \text{ m}^{-1}, \qquad (2)$$

where $ρ_∞$, $U_∞$, and $μ_∞$ are the density, velocity, and dynamic viscosity of the free-stream, respectively. $a_∞ = \sqrt{γRT_∞}$ is the speed of sound based on the free-stream temperature $T_∞$. The flow is considered to be a calorically perfect gas, with specific heat ratio $γ = 1.4$ and gas constant $R = 287.05 \text{ J/(kg} \cdot \text{K)}$. The dynamic viscosity $μ$ is chosen to follow Sutherland's law, and the heat conductivity $κ$ is related to $μ$ through a constant Prandtl number $Pr = 0.72$.

<figure><img src="image_8.png" /><figcaption>Figure 3: Sketch of the physical model and computational domains.</figcaption></figure>

To initialize the forced simulations, a precursor simulation is run with the finite-volume flow solver CEDRE [31]. This first simulation computes the steady base flow over the flat plate with a second-order finite volume spatial discretization, along with a first-order implicit temporal scheme. The computational domain of the steady base flow is referred to as the “low-order” region in Figure 3. Subsequently, the steady flow field is applied as both the initial condition and as the inlet and outlet boundary conditions within a subdomain where we perform the high-fidelity simulations, denoted as the “high-order computational domain” in Figure 3. The details of the domain geometry for this study are summarized in Table 1. We used bilinear interpolation to project the low-order solution onto the high-order domain. The high-fidelity simulations are run with JAGUAR [32, 33, 34, 35], a high-order spectral difference flow solver collaboratively developed by CERFACS and ONERA. At first, JAGUAR is run without any forcing applied to it for long enough to evacuate initial transients. The latter are due to the fact that the two solvers have different truncation errors. A steady-state is then achieved with JAGUAR, which is in excellent agreement with the low-order base flow solution given by CEDRE. It should be noted that in the case of the cold flat plate, the wall temperature is set to $T_w = 300 \text{ K}$ for both CEDRE and JAGUAR. However, for the adiabatic wall, an isothermal wall boundary condition is used in JAGUAR to enforce the wall temperature to $T_w = T_r(x)$, where $T_r(x)$ is the adiabatic wall temperature obtained from CEDRE when enforcing an adiabatic boundary condition.

Once a newly converged steady solution is achieved with Jaguar, it serves as the initial steady state for perturbed unsteady simulations with the same solver. To initiate the Mach second mode instability in the boundary layer, a wall blowing-suction actuator is employed, located within the region between $x_{f1}$ and $x_{f2}$, as depicted in Figure 3. The blowing/suction velocity follows a sinusoidal space-time dependency given by:

$$v'(t) = A_f \sin \left( 2\pi \frac{x - x_{f2}}{x_{f1} - x_{f2}} \right) \sin(\omega t), \quad x_{f1} \le x \le x_{f2} \qquad (3)$$
The value of the frequency *f* is chosen to be 98 kHz while the forcing amplitude *A*<sub>*f*</sub> = 0.1 m/s is set to trigger a linear response of the Mack second mode instability towards the outlet of the high-order computational domain. The forcing frequency comes from the outcome of a linear stability analysis of our base flow with the MAMOUT code [36], which indicates this frequency to be linearly unstable and generating most of its amplification close to *x*<sub>*m*</sub>. The forcing amplitude results from preliminary tests that confirmed the linear response of the Navier-Stokes solution for this amplitude value. The Jaguar simulations use a constant number of cells in the wall-normal direction set to 180 and 1505 cells in the streamwise direction. Given that the solution vector within each cell is expressed through a *p*-order polynomial, which yields *p* + 1 internal degrees of freedom within each cell along each spatial direction, the total number of degrees of freedom along each direction can be obtained by multiplying these values by *p* + 1. The computational grid used for solving the flow is refined near the wall, with a minimum cell size of approximately *Δy* ≈ 25 µm in the wall-normal direction and a constant streamwise spacing of *Δx* ≈ 0.85 mm. With this resolution, there are approximately 80 degrees of freedom within a wavelength of Mack's second mode in the streamwise direction. This oversampling is chosen to ensure that coherent structures are defined by a sufficient number of computational cells in the subsequent coherent structure tracking step.

<table><caption>Table 1: Details on the sizing of the domain geometry in Figure 3.</caption><thead><tr><th>Variable</th><th>Value (m)</th><th>Variable</th><th>Value (m)</th></tr></thead><tbody><tr><td>x<sub>c</sub></td><td>-0.3</td><td>x<sub>m</sub></td><td>1.5</td></tr><tr><td>x<sub>0</sub></td><td>0.2</td><td>y<sub>s</sub></td><td>0.03</td></tr><tr><td>x<sub>f1</sub></td><td>0.2225</td><td>y<sub>r</sub></td><td>0.235</td></tr><tr><td>x<sub>f2</sub></td><td>0.2375</td><td>y<sub>m</sub></td><td>0.3425</td></tr></tbody></table>

### 3. Methodology and Mathematical Description

#### 3.1. Direct Numerical Simulations

The physical problem under investigation is governed by the non-reacting Navier-Stokes equations in 2D generalized coordinates (ξ, η), which consists of the mass conservation equation, the two-component momentum conservation equations, and the energy equation:

$$ \frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{E}}{\partial \xi} + \frac{\partial \mathbf{F}}{\partial \eta} = 0 \qquad (4) $$

Here, the vector **Q** represents the array of conservative variables and **E** and **F** are flux vectors in the ξ and η directions, respectively. These vectors are defined in relation to the corresponding vectors **Q**<sub>c</sub>, **E**<sub>c</sub>, and **F**<sub>c</sub> in a cartesian coordinate system.

$$ \mathbf{Q} = J \mathbf{Q}_c \qquad (5) $$

$$ \mathbf{E} = J \left( \mathbf{E}_c \frac{\partial \xi}{\partial x} + \mathbf{F}_c \frac{\partial \xi}{\partial y} \right) \qquad (6) $$

$$ \mathbf{F} = J \left( \mathbf{E}_c \frac{\partial \eta}{\partial x} + \mathbf{F}_c \frac{\partial \eta}{\partial y} \right) \qquad (7) $$
The value of the frequency *f* is chosen to be 98 kHz while the forcing amplitude *A*<sub>*f*</sub> = 0.1 m/s is set to trigger a linear response of the Mack second mode instability towards the outlet of the high-order computational domain. The forcing frequency comes from the outcome of a linear stability analysis of our base flow with the MAMOUT code [36], which indicates this frequency to be linearly unstable and generating most of its amplification close to *x*<sub>*m*</sub>. The forcing amplitude results from preliminary tests that confirmed the linear response of the Navier-Stokes solution for this amplitude value. The Jaguar simulations use a constant number of cells in the wall-normal direction set to 180 and 1505 cells in the streamwise direction. Given that the solution vector within each cell is expressed through a *p*-order polynomial, which yields *p* + 1 internal degrees of freedom within each cell along each spatial direction, the total number of degrees of freedom along each direction can be obtained by multiplying these values by *p* + 1. The computational grid used for solving the flow is refined near the wall, with a minimum cell size of approximately *Δy* ≈ 25 µm in the wall-normal direction and a constant streamwise spacing of *Δx* ≈ 0.85 mm. With this resolution, there are approximately 80 degrees of freedom within a wavelength of Mack's second mode in the streamwise direction. This oversampling is chosen to ensure that coherent structures are defined by a sufficient number of computational cells in the subsequent coherent structure tracking step.

<table><caption>Table 1: Details on the sizing of the domain geometry in Figure 3.</caption><thead><tr><th>Variable</th><th>Value (m)</th><th>Variable</th><th>Value (m)</th></tr></thead><tbody><tr><td>x<sub>c</sub></td><td>-0.3</td><td>x<sub>m</sub></td><td>1.5</td></tr><tr><td>x<sub>0</sub></td><td>0.2</td><td>y<sub>s</sub></td><td>0.03</td></tr><tr><td>x<sub>f1</sub></td><td>0.2225</td><td>y<sub>r</sub></td><td>0.235</td></tr><tr><td>x<sub>f2</sub></td><td>0.2375</td><td>y<sub>m</sub></td><td>0.3425</td></tr></tbody></table>

### 3. Methodology and Mathematical Description

#### 3.1. Direct Numerical Simulations

The physical problem under investigation is governed by the non-reacting Navier-Stokes equations in 2D generalized coordinates (ξ, η), which consists of the mass conservation equation, the two-component momentum conservation equations, and the energy equation:

$$ \frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{E}}{\partial \xi} + \frac{\partial \mathbf{F}}{\partial \eta} = 0 \qquad (4) $$

Here, the vector **Q** represents the array of conservative variables and **E** and **F** are flux vectors in the ξ and η directions, respectively. These vectors are defined in relation to the corresponding vectors **Q**<sub>c</sub>, **E**<sub>c</sub>, and **F**<sub>c</sub> in a cartesian coordinate system.

$$ \mathbf{Q} = J \mathbf{Q}_c \qquad (5) $$

$$ \mathbf{E} = J \left( \mathbf{E}_c \frac{\partial \xi}{\partial x} + \mathbf{F}_c \frac{\partial \xi}{\partial y} \right) \qquad (6) $$

$$ \mathbf{F} = J \left( \mathbf{E}_c \frac{\partial \eta}{\partial x} + \mathbf{F}_c \frac{\partial \eta}{\partial y} \right) \qquad (7) $$
The symbol $J$ denotes the determinant of the transformation jacobian $J = |\partial(x, y)/\partial(\xi, \eta)|$. The vector of conservative variables and the flux vectors can be expressed in the physical coordinate system as:

$$ \mathbf{Q}_c = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho e \end{pmatrix} \qquad (8) $$

$$ \mathbf{E}_c = \begin{pmatrix} \rho u \\ \rho u^2 + p - \tau_{xx} \\ \rho uv - \tau_{xy} \\ u(\rho e + p) - (u\tau_{xx} + v\tau_{yx} + \kappa \partial_x T) \end{pmatrix} \qquad (9) $$

$$ \mathbf{F}_c = \begin{pmatrix} \rho v \\ \rho uv - \tau_{xy} \\ \rho v^2 + p - \tau_{yy} \\ v(\rho e + p) - (u\tau_{xy} + v\tau_{yy} + \kappa \partial_y T) \end{pmatrix}, \qquad (10) $$

where $\rho$ is the fluid density, $u$ and $v$ are the velocity components along physical coordinates $x$ and $y$, respectively, $p$ is the static pressure, $T$ the temperature, $e = p/\rho(\gamma - 1) + (u^2 + v^2)/2$ is the specific total energy and $\kappa$ the heat conductivity. The components of the viscous stress tensor are labeled $\tau_{ij}$.

The governing equations are solved with Jaguar employing a spectral difference scheme for the spatial discretization. This method relies on a polynomial representation of the solution vector and the flux components inside each cell. The scheme's order $p+1$ is adjustable, with its accuracy level determined by the degree $p$ of the interpolating polynomial. For this study, a fourth-order spectral difference scheme is used while the temporal integration is achieved using a fourth-order five-stage Runge-Kutta time-stepping scheme.

Due to complex nature of the hypersonic boundary layer flow, which encompasses both supersonic ($M > 1$) and subsonic ($M < 1$) regions (see Figure 3), distinct boundary treatments are applied at the inflow and outflow boundaries. Non-reflecting boundary conditions based on Navier-Stokes characteristic boundary conditions [37, 33] are set in the subsonic inflow and outflow regions. For the super/hypersonic inflow ones, Dirichlet conditions are imposed for $u$, $v$, $p$ and $T$, while for the outflow the extrapolated solution from the inner domain is used to define the boundary condition.

## 3.2. Momentum potential decomposition

Doak's momentum potential theory (MPT) enables the decomposition of momentum density $\rho u_i$ into a linear superposition of rotational, acoustic, and thermal components [24]. Following the Helmholtz decomposition, one can first split $\rho u_i$ into a rotational (solenoidal) component $\mathbf{B}$ and an irrotational component expressed as $\nabla \psi$:

$$ \rho \mathbf{u} = \mathbf{B} - \nabla \psi, \quad \nabla \cdot \mathbf{B} = 0. \qquad (11) $$

Assuming the flow to be statistically time-stationary, each instantaneous flow quantity consists of a mean ($\bar{\cdot}$) part and a fluctuation part ($\cdot'$). In this context, the mean scalar potential, denoted as $\bar{\psi}$, is both irrotational and solenoidal and can thus be assumed to be zero, leading to:

$$ \overline{\rho \mathbf{u}} = \overline{\mathbf{B}}, \quad (\rho \mathbf{u})' = \mathbf{B}' - \nabla \psi'. \qquad (12) $$
The symbol $J$ denotes the determinant of the transformation jacobian $J = |\partial(x, y)/\partial(\xi, \eta)|$. The vector of conservative variables and the flux vectors can be expressed in the physical coordinate system as:

$$ \mathbf{Q}_c = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho e \end{pmatrix} \qquad (8) $$

$$ \mathbf{E}_c = \begin{pmatrix} \rho u \\ \rho u^2 + p - \tau_{xx} \\ \rho uv - \tau_{xy} \\ u(\rho e + p) - (u\tau_{xx} + v\tau_{yx} + \kappa \partial_x T) \end{pmatrix} \qquad (9) $$

$$ \mathbf{F}_c = \begin{pmatrix} \rho v \\ \rho uv - \tau_{xy} \\ \rho v^2 + p - \tau_{yy} \\ v(\rho e + p) - (u\tau_{xy} + v\tau_{yy} + \kappa \partial_y T) \end{pmatrix}, \qquad (10) $$

where $\rho$ is the fluid density, $u$ and $v$ are the velocity components along physical coordinates $x$ and $y$, respectively, $p$ is the static pressure, $T$ the temperature, $e = p/\rho(\gamma - 1) + (u^2 + v^2)/2$ is the specific total energy and $\kappa$ the heat conductivity. The components of the viscous stress tensor are labeled $\tau_{ij}$.

The governing equations are solved with Jaguar employing a spectral difference scheme for the spatial discretization. This method relies on a polynomial representation of the solution vector and the flux components inside each cell. The scheme's order $p+1$ is adjustable, with its accuracy level determined by the degree $p$ of the interpolating polynomial. For this study, a fourth-order spectral difference scheme is used while the temporal integration is achieved using a fourth-order five-stage Runge-Kutta time-stepping scheme.

Due to complex nature of the hypersonic boundary layer flow, which encompasses both supersonic ($M > 1$) and subsonic ($M < 1$) regions (see Figure 3), distinct boundary treatments are applied at the inflow and outflow boundaries. Non-reflecting boundary conditions based on Navier-Stokes characteristic boundary conditions [37, 33] are set in the subsonic inflow and outflow regions. For the super/hypersonic inflow ones, Dirichlet conditions are imposed for $u$, $v$, $p$ and $T$, while for the outflow the extrapolated solution from the inner domain is used to define the boundary condition.

### 3.2. Momentum potential decomposition

Doak's momentum potential theory (MPT) enables the decomposition of momentum density $\rho u_i$ into a linear superposition of rotational, acoustic, and thermal components [24]. Following the Helmholtz decomposition, one can first split $\rho u_i$ into a rotational (solenoidal) component $\mathbf{B}$ and an irrotational component expressed as $\nabla \psi$:

$$ \rho \mathbf{u} = \mathbf{B} - \nabla \psi, \quad \nabla \cdot \mathbf{B} = 0. \qquad (11) $$

Assuming the flow to be statistically time-stationary, each instantaneous flow quantity consists of a mean ($\bar{\cdot}$) part and a fluctuation part ($\cdot'$). In this context, the mean scalar potential, denoted as $\bar{\psi}$, is both irrotational and solenoidal and can thus be assumed to be zero, leading to:

$$ \overline{\rho \mathbf{u}} = \overline{\mathbf{B}}, \quad (\rho \mathbf{u})' = \mathbf{B}' - \nabla \psi'. \qquad (12) $$
The fluctuating solenoidal $\mathbf{B}'$ field represents the hydrodynamic component associated with the momentum density field. The continuity equation can then be re-expressed as follows:

$$ \frac{\partial \rho'}{\partial t} + \nabla \cdot (\rho \mathbf{u})' = 0, \quad \nabla \cdot (\overline{\rho \mathbf{u}}) = 0. \qquad (13) $$

By inserting equation (12) in (13), a Poisson equation is derived for the fluctuating scalar potential field:

$$ \nabla^2 \psi' = \frac{\partial \rho'}{\partial t}. \qquad (14) $$

The irrotational field $\psi'$ can now be expressed as the linear superposition of an isentropic (acoustic) $\psi'_A$ and isobaric (thermal, or entropic) $\psi'_T$ contribution as follows:

$$ \psi' = \psi'_A + \psi'_T, \quad \nabla^2 \psi'_A = \left(\frac{1}{a^2}\right) \frac{\partial p'}{\partial t}, \quad \nabla^2 \psi'_T = \left(\frac{\partial \rho}{\partial S}\right) \frac{\partial S'}{\partial t}. \qquad (15) $$

The fluctuating momentum density is therefore depicted as a combination of a vortical component $\mathbf{B}'$, an acoustic component $-\nabla \psi'_A$, and a thermal component $-\nabla \psi'_T$. Here, $a = \sqrt{(\frac{\partial p}{\partial \rho})_S}$ represents the local speed of sound, and $S$ denotes the entropy of the fluid.

In this study, the Poisson equations of the total irrotational component fluctuation $\psi'$ and the acoustic component $\psi'_A$ are solved. Following this, $\psi'_T$ is determined by substitution. On the outer boundaries, thermal fluctuations are considered to be negligible, while strictly vortical disturbances are assumed (and introduced) within the forcing strip. As a result, identical boundary conditions are enforced for both Poisson problems. The Dirichlet boundary conditions for solving the Poisson problems are formulated by integrating along the boundaries, following the procedure detailed in [38, 39, 22]. The Poisson equations are solved using the open-source library Bcube written in the Julia language, employing finite elements (FE) with third-order Lagrange polynomials. The consistency in order across the different interpolations (SD-based Jaguar and FE-based Bcube) was favored with respect to other lower order, finite-volume based Poisson solvers. The impact of the different interpolating schemes on the resulting $\psi'$ fields has not been studied further. The time derivatives appearing on the right-hand sides of the Poisson equations are computed in Jaguar through an upstream fourth-order finite difference scheme.

## 4. Identification of coherent structures

In the present work, the tracking method considered is based on the work by Lozano-Durán et al. [30]. This section presents a definition of coherent structures and outlines the chosen procedure for their identification and extraction.

It is important to recall that the disturbances of interest to us are of the harmonic type. Hence, their characteristic size can be expressed in terms of their wave length $L_\lambda$. For the time tracking of such structures of size roughly equal to $L_\lambda/2$ in length, we estimate that within one forcing period $T_f$ a sampling frequency of $T_f/n$ is sufficient for time-tracking when $n > 2$. We set $n = 18$ as a conservative choice, since the structures may be smaller depending on the threshold. To be more precise, we conduct our analysis choosing a total of $N_s = 365$ flow field snapshots that span approximately 20.6 $T_f$ in time. These snapshots are stored at regular intervals of $\Delta t_s = 5.76 \cdot 10^{-7}$ s, during which an object traveling at the phase speed of the mode disturbance covers a distance equivalent to ca. 4.5 cells with our mesh resolution. Then, we classify a structure as a cluster of spatially adjacent cells, or pixels, inside which a specific quantity of interest $Q$ exceeds a chosen threshold. More specifically, the structures are identified based on satisfying the condition:
The fluctuating solenoidal $\mathbf{B}'$ field represents the hydrodynamic component associated with the momentum density field. The continuity equation can then be re-expressed as follows:

$$ \frac{\partial \rho'}{\partial t} + \nabla \cdot (\rho \mathbf{u})' = 0, \quad \nabla \cdot (\overline{\rho \mathbf{u}}) = 0. \qquad (13) $$

By inserting equation (12) in (13), a Poisson equation is derived for the fluctuating scalar potential field:

$$ \nabla^2 \psi' = \frac{\partial \rho'}{\partial t}. \qquad (14) $$

The irrotational field $\psi'$ can now be expressed as the linear superposition of an isentropic (acoustic) $\psi'_A$ and isobaric (thermal, or entropic) $\psi'_T$ contribution as follows:

$$ \psi' = \psi'_A + \psi'_T, \quad \nabla^2 \psi'_A = \left(\frac{1}{a^2}\right) \frac{\partial p'}{\partial t}, \quad \nabla^2 \psi'_T = \left(\frac{\partial \rho}{\partial S}\right) \frac{\partial S'}{\partial t}. \qquad (15) $$

The fluctuating momentum density is therefore depicted as a combination of a vortical component $\mathbf{B}'$, an acoustic component $-\nabla \psi'_A$, and a thermal component $-\nabla \psi'_T$. Here, $a = \sqrt{\left(\frac{\partial p}{\partial \rho}\right)_S}$ represents the local speed of sound, and $S$ denotes the entropy of the fluid.

In this study, the Poisson equations of the total irrotational component fluctuation $\psi'$ and the acoustic component $\psi'_A$ are solved. Following this, $\psi'_T$ is determined by substitution. On the outer boundaries, thermal fluctuations are considered to be negligible, while strictly vortical disturbances are assumed (and introduced) within the forcing strip. As a result, identical boundary conditions are enforced for both Poisson problems. The Dirichlet boundary conditions for solving the Poisson problems are formulated by integrating along the boundaries, following the procedure detailed in [38, 39, 22]. The Poisson equations are solved using the open-source library Bcube written in the Julia language, employing finite elements (FE) with third-order Lagrange polynomials. The consistency in order across the different interpolations (SD-based Jaguar and FE-based Bcube) was favored with respect to other lower order, finite-volume based Poisson solvers. The impact of the different interpolating schemes on the resulting $\psi'$ fields has not been studied further. The time derivatives appearing on the right-hand sides of the Poisson equations are computed in Jaguar through an upstream fourth-order finite difference scheme.

## 4. Identification of coherent structures

In the present work, the tracking method considered is based on the work by Lozano-Durán et al. [30]. This section presents a definition of coherent structures and outlines the chosen procedure for their identification and extraction.

It is important to recall that the disturbances of interest to us are of the harmonic type. Hence, their characteristic size can be expressed in terms of their wave length $L_\lambda$. For the time tracking of such structures of size roughly equal to $L_\lambda/2$ in length, we estimate that within one forcing period $T_f$ a sampling frequency of $T_f/n$ is sufficient for time-tracking when $n > 2$. We set $n = 18$ as a conservative choice, since the structures may be smaller depending on the threshold. To be more precise, we conduct our analysis choosing a total of $N_s = 365$ flow field snapshots that span approximately 20.6 $T_f$ in time. These snapshots are stored at regular intervals of $\Delta t_s = 5.76 \cdot 10^{-7}$ s, during which an object traveling at the phase speed of the mode disturbance covers a distance equivalent to ca. 4.5 cells with our mesh resolution. Then, we classify a structure as a cluster of spatially adjacent cells, or pixels, inside which a specific quantity of interest $Q$ exceeds a chosen threshold. More specifically, the structures are identified based on satisfying the condition:
$$Q(x, y, t) > \alpha_p H(x). \qquad (16)$$

Here, $\alpha_p$ represents a threshold derived from percolation analysis [40, 41, 27], and $H(x)$ is a shape function that takes into account the streamwise evolution of the perturbation and is defined as:

$$H(x) = \max_{x,y,t}(Q) \frac{\ln(Q'_w)}{\lambda \max_{x,y,t}(\ln(Q'_w))}, \qquad (17)$$

where $\lambda$ a constant set to 1000 and the subscript $\cdot_w$ refers to the wall values. Given that the dynamics of the second Mack mode are confined to a specific region within the boundary layer, any structure with a barycentre position above twice the boundary layer thickness away from the wall will be discarded. Similarly, any structure with an area $l^2$ below a defined threshold, set to $l = 2.2 \cdot 10^{-3}$ m, will be categorized as noise and consequently not recognized as a valid structure. In comparison, note that the boundary layer thickness span the range of values from 3 to 8 mm. This area criterion resulted from an analysis of the distribution of the coherent structure sizes.

# 5. Preliminary results and discussion

## 5.1. Steady state

The u-velocity and temperature profiles obtained from the DNS are plotted at the end of the computational domain in Figure 4. The normalized wall-normal coordinate is determined using the local boundary-layer thickness, which is based on the total enthalpy criterion $h_0(\delta_{h0}) = C_h h_\infty$ [42], with $C_h = 1.008$. Over the insulated wall, the developed boundary layer is thicker, approximately $\delta_{h0} \simeq 14$ mm, in contrast to $\delta_{h0} \simeq 8$ mm over the cold wall. Unnikrishnan et al. [43] noted that second-mode instabilities tend to exhibit maximum amplitudes near high-gradient regions. Consequently, amplification dynamics are intensified closer to the wall under colder surface conditions. Moreover, as can be seen from the right-most panels of Figure 4, the insulated configuration displays a GIP, while wall cooling suppresses the GIP. Similar behavior has been observed in [44, 43].

## 5.2. Linear stability analysis

A spatial local linear stability analysis (see Appendix A) is performed on the HBL profiles for the frequency $f = 98$ kHz. The resulting phase speeds and spatial growth rates associated with F and S modes are plotted in Figure 5. Horizontal dashed lines in the left plots correspond to the continuous branches of the eigenspectra, including the fast acoustic ($c_r/u_e = 1 + 1/M_e$), entropy/vorticity ($c_r/u_e = 1$), and slow acoustic branches ($c_r/u_e = 1 - 1/M_e$). The dashed vertical lines correspond to positions where specific intermodal interaction phenomena can occur due to phase speed synchronization between two modes. For both the adiabatic and cooled walls, modes F and S appear to originate from the fast and slow acoustic branches near the leading edge, respectively. Subsequently, tracking the trajectory of mode F reveals an intersection with the vorticity/entropy branches. This coincides with the point where the amplification rate of mode S starts moving sharply towards the unstable values (negative $\alpha_i$). Further downstream, the phase velocities of modes F and S synchronize for the adiabatic case. At this juncture, mode S seems to reach its maximum instability. Conversely, under the cold wall condition, this synchronization between the two discrete modes is absent, as observed in [15]. Establishing a cause-and-effect relationship between the intermodal interaction phenomena in the eigenspectra and the stability observations of the modes proves to be difficult. Despite this, insights from a range of analytical and numerical receptivity studies [45, 18, 46] indicate that the synchronization phenomena, characterized by phase velocity and eigenfunction interactions between different
$$Q(x, y, t) > \alpha_p H(x). \qquad (16)$$

Here, $\alpha_p$ represents a threshold derived from percolation analysis [40, 41, 27], and $H(x)$ is a shape function that takes into account the streamwise evolution of the perturbation and is defined as:

$$H(x) = \max_{x,y,t}(Q) \frac{\ln(Q'_w)}{\lambda \max_{x,y,t}(\ln(Q'_w))}, \qquad (17)$$

where $\lambda$ a constant set to 1000 and the subscript $\cdot_w$ refers to the wall values. Given that the dynamics of the second Mack mode are confined to a specific region within the boundary layer, any structure with a barycentre position above twice the boundary layer thickness away from the wall will be discarded. Similarly, any structure with an area $l^2$ below a defined threshold, set to $l = 2.2 \cdot 10^{-3}$ m, will be categorized as noise and consequently not recognized as a valid structure. In comparison, note that the boundary layer thickness span the range of values from 3 to 8 mm. This area criterion resulted from an analysis of the distribution of the coherent structure sizes.

# 5. Preliminary results and discussion

## 5.1. Steady state

The u-velocity and temperature profiles obtained from the DNS are plotted at the end of the computational domain in Figure 4. The normalized wall-normal coordinate is determined using the local boundary-layer thickness, which is based on the total enthalpy criterion $h_0(\delta_{h0}) = C_h h_\infty$ [42], with $C_h = 1.008$. Over the insulated wall, the developed boundary layer is thicker, approximately $\delta_{h0} \simeq 14$ mm, in contrast to $\delta_{h0} \simeq 8$ mm over the cold wall. Unnikrishnan et al. [43] noted that second-mode instabilities tend to exhibit maximum amplitudes near high-gradient regions. Consequently, amplification dynamics are intensified closer to the wall under colder surface conditions. Moreover, as can be seen from the right-most panels of Figure 4, the insulated configuration displays a GIP, while wall cooling suppresses the GIP. Similar behavior has been observed in [44, 43].

## 5.2. Linear stability analysis

A spatial local linear stability analysis (see Appendix A) is performed on the HBL profiles for the frequency $f = 98$ kHz. The resulting phase speeds and spatial growth rates associated with F and S modes are plotted in Figure 5. Horizontal dashed lines in the left plots correspond to the continuous branches of the eigenspectra, including the fast acoustic ($c_r/u_e = 1 + 1/M_e$), entropy/vorticity ($c_r/u_e = 1$), and slow acoustic branches ($c_r/u_e = 1 - 1/M_e$). The dashed vertical lines correspond to positions where specific intermodal interaction phenomena can occur due to phase speed synchronization between two modes. For both the adiabatic and cooled walls, modes F and S appear to originate from the fast and slow acoustic branches near the leading edge, respectively. Subsequently, tracking the trajectory of mode F reveals an intersection with the vorticity/entropy branches. This coincides with the point where the amplification rate of mode S starts moving sharply towards the unstable values (negative $\alpha_i$). Further downstream, the phase velocities of modes F and S synchronize for the adiabatic case. At this juncture, mode S seems to reach its maximum instability. Conversely, under the cold wall condition, this synchronization between the two discrete modes is absent, as observed in [15]. Establishing a cause-and-effect relationship between the intermodal interaction phenomena in the eigenspectra and the stability observations of the modes proves to be difficult. Despite this, insights from a range of analytical and numerical receptivity studies [45, 18, 46] indicate that the synchronization phenomena, characterized by phase velocity and eigenfunction interactions between different
<figure><img src="image_8.png" /><figcaption>Figure 4: Streamwise velocity profile (left), temperature profile (center), and indicator function for the existence of a GIP (right). Where the indicator function crosses the zero value (vertical dashed line), a GIP exists. All profiles taken from $x = 1.5$ m for the cold wall (top) and the insulated case (bottom).</figcaption></figure>

modes, imply energy transfer processes between the modes involved. Therefore, we expect the synchronization patterns (or lack thereof) observed through the LST results presented above to be reflected in the FT-based analysis that will be presented later. Indeed, it is the scope of our work to study in detail these interactions.

The total amplification rate of an unstable wave is quantified by the N-factor, representing the overall exponential growth of an unstable perturbation. It is defined in relation to the final amplitude A at location x and the initial amplitude $A_0$ at the critical point $x_{cr}$. The exact definition is given by:

$$N(x, f) = \ln \left( \frac{A}{A_0} \right) = \int_{x_{cr}}^{x} -\alpha_i(x, f) dx. \qquad (18)$$

The N-factor curves associated to the Mack second mode (mode S here) are presented in Figure 6. The amplification of the second Mack mode is more pronounced on the cold wall compared to the adiabatic case, in line with the fact that cooling destabilizes this mode [47]. Moreover, in the adiabatic case, the maximum amplification of this instability occurs more
<figure><img src="image_9.png" /><figcaption>Figure 4: Streamwise velocity profile (left), temperature profile (center), and indicator function for the existence of a GIP (right). Where the indicator function crosses the zero value (vertical dashed line), a GIP exists. All profiles taken from $x = 1.5$ m for the cold wall (top) and the insulated case (bottom).</figcaption></figure>

modes, imply energy transfer processes between the modes involved. Therefore, we expect the synchronization patterns (or lack thereof) observed through the LST results presented above to be reflected in the FT-based analysis that will be presented later. Indeed, it is the scope of our work to study in detail these interactions.

The total amplification rate of an unstable wave is quantified by the N-factor, representing the overall exponential growth of an unstable perturbation. It is defined in relation to the final amplitude A at location x and the initial amplitude $A_0$ at the critical point $x_{cr}$. The exact definition is given by:

$$N(x, f) = \ln \left( \frac{A}{A_0} \right) = \int_{x_{cr}}^{x} -\alpha_i(x, f) dx. \qquad (18)$$

The N-factor curves associated to the Mack second mode (mode S here) are presented in Figure 6. The amplification of the second Mack mode is more pronounced on the cold wall compared to the adiabatic case, in line with the fact that cooling destabilizes this mode [47]. Moreover, in the adiabatic case, the maximum amplification of this instability occurs more
<figure><img src="image_9.png" /><figcaption>Figure 5: Phase velocities (left) and spatial growth rates (right) of discrete modes F and S at a frequency $f = 98$ kHz computed by local LST for the cold wall (top) and the insulated case (bottom).</figcaption></figure>

upstream. This can be explained by the fact that the adiabatic boundary layer is thicker, and that the frequency of the instability is inversely proportional to the thickness $\delta$ [48]. Hence at equal frequency, the amplification peak is shifted upstream in the adiabatic case.

<figure><img src="image_10.png" /><figcaption>Figure 6: N-factor curves computed by local LST for the cold wall and the insulated case.</figcaption></figure>

### 5.3. Coherent structures tracking based on fluctuating vorticity

In this section, the results of the time-tracking of structures associated with the fluctuations of vorticity over the cold wall are discussed.
<figure><img src="image_10.png" /><figcaption>Figure 5: Phase velocities (left) and spatial growth rates (right) of discrete modes F and S at a frequency $f = 98$ kHz computed by local LST for the cold wall (top) and the insulated case (bottom).</figcaption></figure>

upstream. This can be explained by the fact that the adiabatic boundary layer is thicker, and that the frequency of the instability is inversely proportional to the thickness $\delta$ [48]. Hence at equal frequency, the amplification peak is shifted upstream in the adiabatic case.

<figure><img src="image_11.png" /><figcaption>Figure 6: N-factor curves computed by local LST for the cold wall and the insulated case.</figcaption></figure>

### 5.3. Coherent structures tracking based on fluctuating vorticity

In this section, the results of the time-tracking of structures associated with the fluctuations of vorticity over the cold wall are discussed.
### 5.3.1. Percolation Analysis
In order to determine the threshold parameter $\alpha_p$ in equation 16, we conduct a percolation analysis by systematically varying $\alpha_p$ over a range of values. Of particular interest is the evolution of both the number and the area of structures that are identified at the various thresholds. The methodology, as employed in previous studies [41, 30, 49], is described in great detail in [27]. Within the context of our study, the structures based on the fluctuating vorticity are studied in a first instance as a surrogate for the solenoidal component of the Doak decomposition. We note in passing that we have preliminary knowledge suggesting that the perturbations will organize themselves following a pattern of aligned sign-alternating doublets. Hence, we plan to distinguish between structures with positive and negative vorticity perturbations by treating them as distinct entities during the clustering process. A fluctuating vorticity structure is therefore defined as a clustered set of spatially-connected pixels satisfying the following conditions:

$$ 
\begin{aligned}
\nabla \times \mathbf{u}' &< -\alpha_p H(x), \\
\nabla \times \mathbf{u}' &> \alpha_p H(x).
\end{aligned} 
\qquad (19) $$

The results of the percolation analysis are presented in Figure 7 with the ratio of the area of the largest structure $A_{max}$ divided by the total area $A_{tot}$ occupied by all the structures, and the total number of structures $N_s$ normalized by the largest number of identified structures $N_{s,max}$ over all the values of $\alpha_p$. Logically, at a low value of $\alpha_p$, a relatively small number of structures are detected, covering a broad surface area. By increasing the value of $\alpha_p$, these large structures will break down into several smaller but more numerous structures. This process is referred to as the percolation crisis. If the value of $\alpha_p$ is further increased, the ratio $A_{max}/A_{tot}$ tends to increase slowly while the number of objects continues to decrease. We opt for a threshold value of $\alpha_p = 0.4$, located in the range where the area ratio remains roughly constant and the phenomenon of object splitting is still important.

It is worth noting some discrepancies in the percolation curves obtained in this study compared to those commonly reported in the literature. Indeed, typically at a low value of $\alpha_p$, a large number of pixels are clustered into a single large structure, causing the area curve to start from the value of 1, and the ratio of the number of structures to be very low. Moreover, the area ratio eventually tends to 0 as $\alpha_p$ becomes large, as with a very high threshold, no structures are detected. In our case study, the distinctive feature lies in the separate consideration of structures exhibiting positive and negative vorticity fluctuations. Consequently, in regions where perturbations manifest as sign-alternating doublets, the structures are inherently distinct, irrespective of $\alpha_p$, justifying the observed evolution of the curves in Figure 7.

### 5.3.2. Branch classification
The evolution of the structures over their lifetime can be organized into branches. Additionally, these branches can be classified into several groups based on their topology within a temporal graph that gathers them according to the procedure described by Lozano-Durán and Jiménez [30]. A branch designated as “primary” is one associated with a structure that can unequivocally be identified as an individual entity throughout its entire lifespan. Secondary branches, on the other hand, encompass two classifications: “incoming” for instances where the associated object merges with another to form a new entity, and “outgoing” for scenarios involving the splitting of an object. The fourth category, called “connectors,” consists of branches that originate from the splitting of a structure and conclude by merging into another. Table 2 summarizes classification results for our coherent structures based on fluctuating vorticity. We observe a predominance of primary branches. However, this result is challenging to interpret because the saved temporal sequence does not cover a sufficiently long time interval for the perturbations injected at the initial time to reach the end of the computational domain. Indeed, the duration of the temporal sequence is $0.4\tau_c$, where $\tau_c$ represents the advection time defined as the ratio of the average flow velocity outside
### 5.3.1. Percolation Analysis
In order to determine the threshold parameter $\alpha_p$ in equation 16, we conduct a percolation analysis by systematically varying $\alpha_p$ over a range of values. Of particular interest is the evolution of both the number and the area of structures that are identified at the various thresholds. The methodology, as employed in previous studies [41, 30, 49], is described in great detail in [27]. Within the context of our study, the structures based on the fluctuating vorticity are studied in a first instance as a surrogate for the solenoidal component of the Doak decomposition. We note in passing that we have preliminary knowledge suggesting that the perturbations will organize themselves following a pattern of aligned sign-alternating doublets. Hence, we plan to distinguish between structures with positive and negative vorticity perturbations by treating them as distinct entities during the clustering process. A fluctuating vorticity structure is therefore defined as a clustered set of spatially-connected pixels satisfying the following conditions:

$$ 
\begin{aligned}
\nabla \times \mathbf{u}' &< -\alpha_p H(x), \\
\nabla \times \mathbf{u}' &> \alpha_p H(x).
\end{aligned} 
\qquad (19) $$

The results of the percolation analysis are presented in Figure 7 with the ratio of the area of the largest structure $A_{max}$ divided by the total area $A_{tot}$ occupied by all the structures, and the total number of structures $N_s$ normalized by the largest number of identified structures $N_{s,max}$ over all the values of $\alpha_p$. Logically, at a low value of $\alpha_p$, a relatively small number of structures are detected, covering a broad surface area. By increasing the value of $\alpha_p$, these large structures will break down into several smaller but more numerous structures. This process is referred to as the percolation crisis. If the value of $\alpha_p$ is further increased, the ratio $A_{max}/A_{tot}$ tends to increase slowly while the number of objects continues to decrease. We opt for a threshold value of $\alpha_p = 0.4$, located in the range where the area ratio remains roughly constant and the phenomenon of object splitting is still important.

It is worth noting some discrepancies in the percolation curves obtained in this study compared to those commonly reported in the literature. Indeed, typically at a low value of $\alpha_p$, a large number of pixels are clustered into a single large structure, causing the area curve to start from the value of 1, and the ratio of the number of structures to be very low. Moreover, the area ratio eventually tends to 0 as $\alpha_p$ becomes large, as with a very high threshold, no structures are detected. In our case study, the distinctive feature lies in the separate consideration of structures exhibiting positive and negative vorticity fluctuations. Consequently, in regions where perturbations manifest as sign-alternating doublets, the structures are inherently distinct, irrespective of $\alpha_p$, justifying the observed evolution of the curves in Figure 7.

### 5.3.2. Branch classification
The evolution of the structures over their lifetime can be organized into branches. Additionally, these branches can be classified into several groups based on their topology within a temporal graph that gathers them according to the procedure described by Lozano-Durán and Jiménez [30]. A branch designated as “primary” is one associated with a structure that can unequivocally be identified as an individual entity throughout its entire lifespan. Secondary branches, on the other hand, encompass two classifications: “incoming” for instances where the associated object merges with another to form a new entity, and “outgoing” for scenarios involving the splitting of an object. The fourth category, called “connectors,” consists of branches that originate from the splitting of a structure and conclude by merging into another. Table 2 summarizes classification results for our coherent structures based on fluctuating vorticity. We observe a predominance of primary branches. However, this result is challenging to interpret because the saved temporal sequence does not cover a sufficiently long time interval for the perturbations injected at the initial time to reach the end of the computational domain. Indeed, the duration of the temporal sequence is $0.4\tau_c$, where $\tau_c$ represents the advection time defined as the ratio of the average flow velocity outside
<figure><img src="image_11.png" /><figcaption>Figure 7: Percolation analysis curves: the solid red line represents the ratio between the area of the largest structure, $A_{max}$, and the total volume occupied by structures, $A_{tot}$. The solid blue line depicts the ratio between the number of structures, $N_s$, and the largest number of identified structures, $N_{s,max}$.</figcaption></figure>

of the boundary layer to the distance between the blowing/suction slot and the end of the computational domain. Therefore, it is not possible to determine the lifespan of a branch from its birth to its death when it lives for longer than $0.4\tau_c$. However, it is noteworthy that outgoing branches are so underrepresented compared to the other branch types within the dynamics of the vorticity fluctuation structures. Additionally, the incoming branches appear to be relatively significant in terms of number and occupied area, showing that mergers occur very frequently. Overall, however, we conclude that branch classification requires a longer time series to show definite results and shed light on the merging/splitting processes. Future work will carry out coherent structure tracking based on MPT components over a longer time series.

<table><caption>Table 2: Branch number fraction and branch area fraction displayed for each branch type.</caption><thead><tr><td>Branch type</td><td>$N_b/N_{b,tot}$</td><td>$A_b/A_{b,tot}$</td></tr></thead><tbody><tr><td>Primary</td><td>0.498</td><td>0.447</td></tr><tr><td>Incoming</td><td>0.255</td><td>0.408</td></tr><tr><td>Outgoing</td><td>0.043</td><td>0.030</td></tr><tr><td>Connectors</td><td>0.204</td><td>0.115</td></tr></tbody></table>

### 5.3.3. Vertical evolution

Geometrical properties are derived from the temporal evolution along the identified branches. This section investigates the dynamics of vorticity fluctuation structures along the wall-normal direction by analyzing geometrical parameters, including the minimum and maximum y coordinates along the streamwise direction ($y_{min}$ and $y_{max}$, respectively). Additionally, the study encompasses the structures' heights, represented by $L_y = y_{max} - y_{min}$, and the y-coordinate $y_{cg}$ of the structures' barycentre. Figure 8 illustrates the spatial evolution of these geometric properties of the structures. Dashed horizontal lines indicate, from top to bottom, the boundary layer edge, the critical layer, and the relative sonic line. Additionally, solid vertical lines mark positions where key phenomena in the structures' behavior are observed,
<figure><img src="image_12.png" /><figcaption>Figure 7: Percolation analysis curves: the solid red line represents the ratio between the area of the largest structure, $A_{max}$, and the total volume occupied by structures, $A_{tot}$. The solid blue line depicts the ratio between the number of structures, $N_s$, and the largest number of identified structures, $N_{s,max}$.</figcaption></figure>

of the boundary layer to the distance between the blowing/suction slot and the end of the computational domain. Therefore, it is not possible to determine the lifespan of a branch from its birth to its death when it lives for longer than $0.4\tau_c$. However, it is noteworthy that outgoing branches are so underrepresented compared to the other branch types within the dynamics of the vorticity fluctuation structures. Additionally, the incoming branches appear to be relatively significant in terms of number and occupied area, showing that mergers occur very frequently. Overall, however, we conclude that branch classification requires a longer time series to show definite results and shed light on the merging/splitting processes. Future work will carry out coherent structure tracking based on MPT components over a longer time series.

<table><caption>Table 2: Branch number fraction and branch area fraction displayed for each branch type.</caption><thead><tr><th>Branch type</th><th>N<sub>b</sub>/N<sub>b,tot</sub></th><th>A<sub>b</sub>/A<sub>b,tot</sub></th></tr></thead><tbody><tr><td>Primary</td><td>0.498</td><td>0.447</td></tr><tr><td>Incoming</td><td>0.255</td><td>0.408</td></tr><tr><td>Outgoing</td><td>0.043</td><td>0.030</td></tr><tr><td>Connectors</td><td>0.204</td><td>0.115</td></tr></tbody></table>

### 5.3.3. Vertical evolution

Geometrical properties are derived from the temporal evolution along the identified branches. This section investigates the dynamics of vorticity fluctuation structures along the wall-normal direction by analyzing geometrical parameters, including the minimum and maximum y coordinates along the streamwise direction ($y_{min}$ and $y_{max}$, respectively). Additionally, the study encompasses the structures' heights, represented by $L_y = y_{max} - y_{min}$, and the y-coordinate $y_{cg}$ of the structures' barycentre. Figure 8 illustrates the spatial evolution of these geometric properties of the structures. Dashed horizontal lines indicate, from top to bottom, the boundary layer edge, the critical layer, and the relative sonic line. Additionally, solid vertical lines mark positions where key phenomena in the structures' behavior are observed,
which will be clarified in the discussion that follows. A color code is also employed to categorize structures with similar geometrical properties. Again, those properties will become apparent below.

Observations reveal that slightly downstream from the perturbation injection slit, structures attached to the wall coexist with other structures situated slightly above. Figure 9 provides a snapshot view of the vorticity fluctuation contours zoomed in on this area, providing qualitative insights into the behavior of the structures. The presence of elongated structures of vorticity fluctuations with opposite signs overlaid is noticeable, indicating the presence of a shear layer associated with the imposition of a blowing/suction condition at the wall. After $x = 0.32$ m, a point marked by the first vertical line, the structures appear to organize into different patterns. Two groups of structures are distinguishable in Figure 8: a first group of objects (in green) attached to the wall and extending up to the relative sonic line, and a second group of objects (in blue) between the sonic line and the boundary layer edge. It is also evident from the $y_{cg}$ locations (see Figure 8) that this second group of objects organizes around the critical layer. This spatial arrangement is clearly depicted in Figure 9, which provides a zoomed-in view of this area. The presence of sign-alternating doublets is observed, with two regions of opposite vorticity perturbation: one trapped between the wall and the relative sonic line, and the other around the critical layer. It is noted that the behavior of vorticity fluctuation structures varies according to the sign of the fluctuation. Structures exhibiting negative fluctuations display a more significant size near the wall compared to their counterparts with positive fluctuations. Conversely, this trend is reversed for the alternate set of structures situated above the relative sonic line. From the position $x = 0.83$ m (marked by the second vertical line), corresponding to the point at which the S mode of linear stability analysis becomes unstable (see Figure 5), a change in the behavior of these structures becomes observable. The structures aligned with the critical layer become thicker over time, mainly because they extend downward while remaining aligned with the critical layer (see Figure 8). This phenomenon reaches a critical point at $x = 1.05$ m, marked by the third vertical line, corresponding to the point where the S mode is most unstable according to the LST (see Figure 5). At this point, the upper structures appear to merge with objects attached to the wall, and the new structures obtained then extend from the wall to the boundary layer edge (in yellow). In Figure 9, it is observed that this merging occurs between regions of vorticity fluctuation with the same sign but located at two different heights. These observations regarding the behavior of vorticity fluctuation structures seem to indicate the generation of a shear layer that results from applying a perturbation to the wall in the form of a blowing/suction source. This layer eventually gives rise to sign-alternating doublets exhibiting trapped wave behavior between the wall and the relative sonic line. Additionally, a zone of vorticity perturbation organizes around the critical layer. Finally, we observe that the positions marking changes in the stability of the S mode from the LST play a crucial role in the dynamics of vorticity fluctuation structures.

5.3.4. *Advection velocities* The streamwise advection velocity of the structures is computed based on the temporal displacement of the structures' streamwise barycenter position, $x_{cg}$. In Figure 10, we present a comparative analysis of the advection velocities of vorticity fluctuation structures and the phase speeds of the discrete modes F and S obtained from LST analysis. Interestingly, in the proximity of the excitation point, the structures seem to be advected at the bulk flow velocity $u_e$. Shortly after, their advection velocity follows a trend similar to that of mode F, including the presence of multiple peaks. These could be indicative of the superposition of multiple modes associated with the blowing/suction forcing, which excites mode F, mode S, acoustic waves, and entropy/vorticity waves simultaneously just downstream of the actuator [50]. What can be stated based on our current results is that this region of varying convection speed also corresponds to significant topological changes in the structure's shape as can be seen
which will be clarified in the discussion that follows. A color code is also employed to categorize structures with similar geometrical properties. Again, those properties will become apparent below.

Observations reveal that slightly downstream from the perturbation injection slit, structures attached to the wall coexist with other structures situated slightly above. Figure 9 provides a snapshot view of the vorticity fluctuation contours zoomed in on this area, providing qualitative insights into the behavior of the structures. The presence of elongated structures of vorticity fluctuations with opposite signs overlaid is noticeable, indicating the presence of a shear layer associated with the imposition of a blowing/suction condition at the wall. After $x = 0.32$ m, a point marked by the first vertical line, the structures appear to organize into different patterns. Two groups of structures are distinguishable in Figure 8: a first group of objects (in green) attached to the wall and extending up to the relative sonic line, and a second group of objects (in blue) between the sonic line and the boundary layer edge. It is also evident from the $y_{cg}$ locations (see Figure 8) that this second group of objects organizes around the critical layer. This spatial arrangement is clearly depicted in Figure 9, which provides a zoomed-in view of this area. The presence of sign-alternating doublets is observed, with two regions of opposite vorticity perturbation: one trapped between the wall and the relative sonic line, and the other around the critical layer. It is noted that the behavior of vorticity fluctuation structures varies according to the sign of the fluctuation. Structures exhibiting negative fluctuations display a more significant size near the wall compared to their counterparts with positive fluctuations. Conversely, this trend is reversed for the alternate set of structures situated above the relative sonic line. From the position $x = 0.83$ m (marked by the second vertical line), corresponding to the point at which the S mode of linear stability analysis becomes unstable (see Figure 5), a change in the behavior of these structures becomes observable. The structures aligned with the critical layer become thicker over time, mainly because they extend downward while remaining aligned with the critical layer (see Figure 8). This phenomenon reaches a critical point at $x = 1.05$ m, marked by the third vertical line, corresponding to the point where the S mode is most unstable according to the LST (see Figure 5). At this point, the upper structures appear to merge with objects attached to the wall, and the new structures obtained then extend from the wall to the boundary layer edge (in yellow). In Figure 9, it is observed that this merging occurs between regions of vorticity fluctuation with the same sign but located at two different heights. These observations regarding the behavior of vorticity fluctuation structures seem to indicate the generation of a shear layer that results from applying a perturbation to the wall in the form of a blowing/suction source. This layer eventually gives rise to sign-alternating doublets exhibiting trapped wave behavior between the wall and the relative sonic line. Additionally, a zone of vorticity perturbation organizes around the critical layer. Finally, we observe that the positions marking changes in the stability of the S mode from the LST play a crucial role in the dynamics of vorticity fluctuation structures.

5.3.4. *Advection velocities* The streamwise advection velocity of the structures is computed based on the temporal displacement of the structures' streamwise barycenter position, $x_{cg}$. In Figure 10, we present a comparative analysis of the advection velocities of vorticity fluctuation structures and the phase speeds of the discrete modes F and S obtained from LST analysis. Interestingly, in the proximity of the excitation point, the structures seem to be advected at the bulk flow velocity $u_e$. Shortly after, their advection velocity follows a trend similar to that of mode F, including the presence of multiple peaks. These could be indicative of the superposition of multiple modes associated with the blowing/suction forcing, which excites mode F, mode S, acoustic waves, and entropy/vorticity waves simultaneously just downstream of the actuator [50]. What can be stated based on our current results is that this region of varying convection speed also corresponds to significant topological changes in the structure's shape as can be seen
<figure><img src="image_12.png" /><figcaption>Figure 8: Spatial evolution of the geometrical features related to the vertical extent and position of vortical coherent structures: maximum wall-normal coordinate (top left), minimum wall-normal coordinate (bottom left), height (top right), wall-normal barycenter position (bottom right). The vertical solid lines represent the delimitation of region of interest I, II, III and IV. The dashed horizontal lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively. Structures within each region are represented by distinct color: ● for region I, ● for regions II and III attached to the wall, ● for regions II and III not attached to the wall, and • for region IV. Additionally, ⊕ and ⊙ denote structures associated with positive and negative fluctuations, respectively, where the distinction is relevant.</figcaption></figure>

on Figure 8.

Further downstream, the structures are advected with a velocity close to the phase velocity of mode S. This transition appears to occur around $x = 0.83$ m, corresponding to the point where the S mode becomes unstable according to linear stability analysis (see Figure 5), and also coincides with a change in the vertical evolution behavior of the structures (Figure 8). Such behavior appears to be explained by the stability of different modes present in the boundary layer. Near the actuator, the modes of the fast acoustic branch and the vorticity/entropy branch are excited in addition to the F and S modes, as mentioned earlier. Further downstream, the F mode is relatively less stable (growth rate close to 0), along with certain hydrodynamic and entropic modes. The perturbations propagating into the boundary layer are then dominated by the latter. Finally, the S mode becomes unstable and eventually prevails. A comparison of the wall-normal perturbation amplitudes obtained from DNS and the eigenfunctions of the discrete modes F and S obtained by LST at two different locations, $x = 0.55$ m and $x = 1.35$ m, is presented in Figure 11. The consistency between the wall-normal perturbation profiles and the eigenfunctions of mode F at $x = 0.55$ m suggests that the fluctuations are dominated by this mode in this region. At $x = 1.35$ m, the trend reverses, and the perturbation profiles align closely with the eigenfunctions of mode S. These results support the analysis conducted
<figure><img src="image_13.png" /></figure>

Figure 8: Spatial evolution of the geometrical features related to the vertical extent and position of vortical coherent structures: maximum wall-normal coordinate (top left), minimum wall-normal coordinate (bottom left), height (top right), wall-normal barycenter position (bottom right). The vertical solid lines represent the delimitation of region of interest I, II, III and IV. The dashed horizontal lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively. Structures within each region are represented by distinct color: ● for region I, ● for regions II and III attached to the wall, ● for regions II and III not attached to the wall, and • for region IV. Additionally, ⊕ and ⊙ denote structures associated with positive and negative fluctuations, respectively, where the distinction is relevant.

on Figure 8.

Further downstream, the structures are advected with a velocity close to the phase velocity of mode S. This transition appears to occur around $x = 0.83$ m, corresponding to the point where the S mode becomes unstable according to linear stability analysis (see Figure 5), and also coincides with a change in the vertical evolution behavior of the structures (Figure 8). Such behavior appears to be explained by the stability of different modes present in the boundary layer. Near the actuator, the modes of the fast acoustic branch and the vorticity/entropy branch are excited in addition to the F and S modes, as mentioned earlier. Further downstream, the F mode is relatively less stable (growth rate close to 0), along with certain hydrodynamic and entropic modes. The perturbations propagating into the boundary layer are then dominated by the latter. Finally, the S mode becomes unstable and eventually prevails. A comparison of the wall-normal perturbation amplitudes obtained from DNS and the eigenfunctions of the discrete modes F and S obtained by LST at two different locations, $x = 0.55$ m and $x = 1.35$ m, is presented in Figure 11. The consistency between the wall-normal perturbation profiles and the eigenfunctions of mode F at $x = 0.55$ m suggests that the fluctuations are dominated by this mode in this region. At $x = 1.35$ m, the trend reverses, and the perturbation profiles align closely with the eigenfunctions of mode S. These results support the analysis conducted
<figure><img src="image_13.png" /><figcaption>Figure 9: Instantaneous magnified view of vorticity fluctuation in three zones of interest I, II and IV. The black dashed lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively.</figcaption></figure>

on the comparison of the phase velocities of the two discrete modes from LST and the advection velocity of the fluctuating vorticity structures.

<figure><img src="image_14.png" /><figcaption>Figure 10: Phase speed velocity of the modes present in the cooled-wall boundary layer obtained through LST, and mean advection velocity of the vorticity structures found through coherent structure tracking.</figcaption></figure>

5.4. *Fluid-thermodynamic decomposition*

We apply Doak's MPT decomposition to instantaneous flow fields of fluctuating quantities extracted from the DNS (cooled and adiabatic walls). Figure 12 presents snapshots of the
<figure><img src="image_14.png" /><figcaption>Figure 9: Instantaneous magnified view of vorticity fluctuation in three zones of interest I, II and IV. The black dashed lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively.</figcaption></figure>

on the comparison of the phase velocities of the two discrete modes from LST and the advection velocity of the fluctuating vorticity structures.

<figure><img src="image_15.png" /><figcaption>Figure 10: Phase speed velocity of the modes present in the cooled-wall boundary layer obtained through LST, and mean advection velocity of the vorticity structures found through coherent structure tracking.</figcaption></figure>

5.4. *Fluid-thermodynamic decomposition*

We apply Doak's MPT decomposition to instantaneous flow fields of fluctuating quantities extracted from the DNS (cooled and adiabatic walls). Figure 12 presents snapshots of the
<figure><img src="image_15.png" /><figcaption>Figure 11: Wall-normal perturbation amplitudes from DNS and F and S mode eigenfunctions obtained by LST at $x = 0.55$ m (left) and $x = 1.35$ m (right). The black dashed lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively.</figcaption></figure>

magnitudes of the instantaneous vortical (||$B$||), acoustic (||$\nabla\psi_A$||), and thermal (||$\nabla\psi_T$||) fluctuations over the cold plate surface, shown for the entire computational domain. It is clear from the contours that within the region where the second Mack mode experiences maximum amplification (approximately $x = 1.35$ m as indicated by the LST), all three FT components undergo amplification within the boundary layer, while outside of the boundary layer the vortical and the acoustic components are the most intense. The three components maintain a relatively low intensity until reaching the amplification region of the second Mack mode. The vortical and acoustic components feature radiation outward from the boundary layer, guided by the oblique shocklet induced by the blowing/suction source term. Simultaneously, this shocklet appears to generate vortical fluctuations. In contrast, the thermal component is concentrated mainly within the boundary layer, and displays no radiation tendency away from the boundary layer. This observation aligns well with findings from previous studies [22, 26]. Similar observations were made for the insulated flat plate case (not shown for brevity).

Figure 13 allows us to be more quantitative, illustrating the magnitudes of the FT components at the wall surface for both the cold wall and the insulated wall cases. The first thing to note is that the magnitudes are significantly lower over the adiabatic wall. This trend follows from the fact that perturbations are able to grow significantly more in the cooled wall case, as shown by the N-factor curves in Figure 6. Recalling from equation 18 that the amplitude grows as $e^N$, differences in FT-component magnitudes of one order of magnitude can be expected and are indeed observed.

The region of growth leading to the maximum magnitudes in these components is clearly aligned with the location where the second Mack mode predicted by LST is most unstable in both cases (see Figure 5). Moreover, the acoustic component is the most intense of the three for both wall temperatures, followed by the thermal component, while the vortical component experiences the mildest amplification. In the case of the cold wall, the acoustic component appears to be much more predominant with respect to the other two components in the second Mack mode amplification region. This observation is consistent with the destabilizing influence of wall cooling and with the acoustic nature of this mode. The agreement observed between our findings and the prior research conducted by Unnikrishnan *et al.* [22] and Long *et al.* [26] validates our implementation of the Doak decomposition on flow fields of hypersonic boundary layers.
<figure><img src="image_16.png" /><figcaption>Figure 11: Wall-normal perturbation amplitudes from DNS and F and S mode eigenfunctions obtained by LST at $x = 0.55$ m (left) and $x = 1.35$ m (right). The black dashed lines correspond to the boundary layer edge, the critical layer, and the relative sonic line, from top to bottom, respectively.</figcaption></figure>

magnitudes of the instantaneous vortical (||$B$||), acoustic (||$\nabla\psi_A$||), and thermal (||$\nabla\psi_T$||) fluctuations over the cold plate surface, shown for the entire computational domain. It is clear from the contours that within the region where the second Mack mode experiences maximum amplification (approximately $x = 1.35$ m as indicated by the LST), all three FT components undergo amplification within the boundary layer, while outside of the boundary layer the vortical and the acoustic components are the most intense. The three components maintain a relatively low intensity until reaching the amplification region of the second Mack mode. The vortical and acoustic components feature radiation outward from the boundary layer, guided by the oblique shocklet induced by the blowing/suction source term. Simultaneously, this shocklet appears to generate vortical fluctuations. In contrast, the thermal component is concentrated mainly within the boundary layer, and displays no radiation tendency away from the boundary layer. This observation aligns well with findings from previous studies [22, 26]. Similar observations were made for the insulated flat plate case (not shown for brevity).

Figure 13 allows us to be more quantitative, illustrating the magnitudes of the FT components at the wall surface for both the cold wall and the insulated wall cases. The first thing to note is that the magnitudes are significantly lower over the adiabatic wall. This trend follows from the fact that perturbations are able to grow significantly more in the cooled wall case, as shown by the N-factor curves in Figure 6. Recalling from equation 18 that the amplitude grows as $e^N$, differences in FT-component magnitudes of one order of magnitude can be expected and are indeed observed.

The region of growth leading to the maximum magnitudes in these components is clearly aligned with the location where the second Mack mode predicted by LST is most unstable in both cases (see Figure 5). Moreover, the acoustic component is the most intense of the three for both wall temperatures, followed by the thermal component, while the vortical component experiences the mildest amplification. In the case of the cold wall, the acoustic component appears to be much more predominant with respect to the other two components in the second Mack mode amplification region. This observation is consistent with the destabilizing influence of wall cooling and with the acoustic nature of this mode. The agreement observed between our findings and the prior research conducted by Unnikrishnan *et al.* [22] and Long *et al.* [26] validates our implementation of the Doak decomposition on flow fields of hypersonic boundary layers.
<figure><img src="image_16.png" /><figcaption>Figure 12: Instantaneous view of magnitudes of vortical (top), acoustic (center) and thermal (bottom) perturbations over the cold plate surface.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 13: Magnitudes of the fluid-thermodynamic components at $y = 0$ in the case of the cold wall (left) and the adiabatic wall (right).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 12: Instantaneous view of magnitudes of vortical (top), acoustic (center) and thermal (bottom) perturbations over the cold plate surface.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 13: Magnitudes of the fluid-thermodynamic components at $y = 0$ in the case of the cold wall (left) and the adiabatic wall (right).</figcaption></figure>
## 6. Conclusion

In this study, we have delved into the dynamics of high-speed boundary layer flows over insulated and cooled walls, exploring the link between stability analysis and the properties of coherent structures of vorticity found in the boundary layer. This study is an initial exploration on how coherent structure tracking methods and the FT component decomposition of the flow field can uncover the mechanisms underlying the second Mack mode instability.

A local linear stability analysis lead to the expected distinct behavior in the branching of discrete modes for adiabatic and cooled wall conditions, with an absence of synchronization in the cooled wall case. Coherent structure tracking based on fluctuating vorticity provided a qualitative examination of the spatial and temporal evolution of structures within the boundary layer. The spatial evolution of coherent structures aligns with changes in the linear stability of the discrete modes, underscoring the interplay between coherent structures and the underlying instability mechanisms in high-speed boundary layers. A noteworthy result stems from the evolution of the vertical dimension and location of the structures, revealing the coexistence of vorticity structures at different heights along the wall that feature alternating sign doublets during mode F predominance. These structures organize around the critical layer, and above the relative sonic line. Structure mergers at a critical point aligned with the maximum instability of the S mode has also been revealed, the dynamical implications of which need to be further investigated.

A classification of the coherent structure dynamics into branches within graphs was attempted. It revealed that for this time-periodic flow, a longer time sequence is necessary in order for the objects that appear during the simulation to have the time to be washed out of the domain.

The MPT-based FT decomposition introduced by Doak has been implemented as a post-processing tool to our boundary layer data, which is the preliminary step towards time-tracking of coherent structures based on this quantity. Individual snapshots confirm the insights into the different energetic modes within the boundary layer that have been observed previously. Here, we have linked the quantitative trends in the acoustic and thermal components of the FT decomposition with the amplification region of the second Mack mode instability. The acoustic component, particularly pronounced in the cold wall case, emphasizes the destabilizing influence of wall cooling. Future work will involve carrying out the time-tracking used in this study for vortical structures but applied to FT-based objects instead, and on the basis of a data set spanning a longer time series.

## Acknowledgments

This work was supported in part by the European Research Council under the Caust grant ERC-AdG-101018287. The authors express their gratitude to the developers of Bcube at ONERA/DMPE, namely Maxime Bouyges, Ghislain Blanchard, and Lokman Bennani, for their valuable assistance in configuring and managing the solution of the Poisson problems.

## Appendix A. Linear Stability Theory

The Linear Stability Theory (LST) equations are derived from the governing Navier-Stokes equations. Initially, LST assumes that the instantaneous flow variables $q = (p, u, v, w, T)^T$ can be decomposed into mean and fluctuating components, represented as $\bar{q}$ and $q'$, respectively, where $q = \bar{q} + q'$. In cases of natural transition, the disturbances are assumed to be very small, which allows the linearization of the Navier-Stokes equations, considering that quadratic and higher-order terms are negligible. The parallel or quasi-parallel flow assumption is applied, i.e. the mean flow is considered as predominantly dependent on the wall-normal direction $y$. Under these conditions, normal mode solutions are sought in the form:
## 6. Conclusion

In this study, we have delved into the dynamics of high-speed boundary layer flows over insulated and cooled walls, exploring the link between stability analysis and the properties of coherent structures of vorticity found in the boundary layer. This study is an initial exploration on how coherent structure tracking methods and the FT component decomposition of the flow field can uncover the mechanisms underlying the second Mack mode instability.

A local linear stability analysis lead to the expected distinct behavior in the branching of discrete modes for adiabatic and cooled wall conditions, with an absence of synchronization in the cooled wall case. Coherent structure tracking based on fluctuating vorticity provided a qualitative examination of the spatial and temporal evolution of structures within the boundary layer. The spatial evolution of coherent structures aligns with changes in the linear stability of the discrete modes, underscoring the interplay between coherent structures and the underlying instability mechanisms in high-speed boundary layers. A noteworthy result stems from the evolution of the vertical dimension and location of the structures, revealing the coexistence of vorticity structures at different heights along the wall that feature alternating sign doublets during mode F predominance. These structures organize around the critical layer, and above the relative sonic line. Structure mergers at a critical point aligned with the maximum instability of the S mode has also been revealed, the dynamical implications of which need to be further investigated.

A classification of the coherent structure dynamics into branches within graphs was attempted. It revealed that for this time-periodic flow, a longer time sequence is necessary in order for the objects that appear during the simulation to have the time to be washed out of the domain.

The MPT-based FT decomposition introduced by Doak has been implemented as a post-processing tool to our boundary layer data, which is the preliminary step towards time-tracking of coherent structures based on this quantity. Individual snapshots confirm the insights into the different energetic modes within the boundary layer that have been observed previously. Here, we have linked the quantitative trends in the acoustic and thermal components of the FT decomposition with the amplification region of the second Mack mode instability. The acoustic component, particularly pronounced in the cold wall case, emphasizes the destabilizing influence of wall cooling. Future work will involve carrying out the time-tracking used in this study for vortical structures but applied to FT-based objects instead, and on the basis of a data set spanning a longer time series.

## Acknowledgments

This work was supported in part by the European Research Council under the Caust grant ERC-AdG-101018287. The authors express their gratitude to the developers of Bcube at ONERA/DMPE, namely Maxime Bouyges, Ghislain Blanchard, and Lokman Bennani, for their valuable assistance in configuring and managing the solution of the Poisson problems.

## Appendix A. Linear Stability Theory

The Linear Stability Theory (LST) equations are derived from the governing Navier-Stokes equations. Initially, LST assumes that the instantaneous flow variables $q = (p, u, v, w, T)^T$ can be decomposed into mean and fluctuating components, represented as $\bar{q}$ and $q'$, respectively, where $q = \bar{q} + q'$. In cases of natural transition, the disturbances are assumed to be very small, which allows the linearization of the Navier-Stokes equations, considering that quadratic and higher-order terms are negligible. The parallel or quasi-parallel flow assumption is applied, i.e. the mean flow is considered as predominantly dependent on the wall-normal direction $y$. Under these conditions, normal mode solutions are sought in the form:
$$ \mathbf{q}'(x, y, t) = \hat{\mathbf{q}}(y) \exp [i(\alpha x + \beta z - \omega t)] + \text{c.c.} \quad (A.1) $$

Here, $\alpha$ and $\beta$ represent the $x$ and $z$ components of the wavenumber vector $\mathbf{k}$, while $\omega$ is the frequency, and $\hat{\mathbf{q}}(y)$ represents eigenfunctions that describe the mode structure within the boundary layer. By introducing the normal modes into the linearized Navier-Stokes equations, a one-dimensional differential eigenvalue problem for the modal functions is obtained:

$$ \left( \mathbf{A} \frac{d}{dy^2} + \mathbf{B} \frac{d}{dy} + \mathbf{C} \right) \hat{\mathbf{q}} = 0, \quad (A.2) $$

Here, $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are complex square matrices of size 5. In the context of a boundary layer flow, the perturbation field can be described by the following boundary conditions:

$$ \text{At } y = 0, \quad \hat{u} = \hat{v} = \hat{w} = 0, \quad (A.3) $$

$$ \hat{T} = 0 \quad \text{for an isothermal wall,} \quad (A.4) $$

$$ \frac{d\hat{T}}{dy} = 0 \quad \text{for an adiabatic wall,} \quad (A.5) $$

$$ \text{At } y \to \infty, \quad \hat{u} = \hat{v} = \hat{w} = \hat{T} \to 0 \quad (A.6) $$

For spatial stability analysis, the frequency $\omega$ is real, whereas $\alpha$ is a complex eigenvalue, with $\alpha = \alpha_r + i\alpha_i$. When $\alpha_i < 0$, the flow is unstable, exhibiting a spatial growth rate $\alpha_i$ and a streamwise wavenumber $\alpha_r$. The local phase velocity is defined as $c_r = \omega/\alpha_r$. Additionally, in the context of the Mack second mode, which primarily exhibits two-dimensional dynamics, $\beta$ is assumed to be 0 in this study. The relationship for the eigenvalues, often referred to as the dispersion relation, can be expressed as:

$$ \alpha = \Omega(\omega). \quad (A.7) $$

The evaluation of the dispersion relation for a boundary layer profile is a fundamental aspect of stability theory. The eigenvalues, along with their corresponding eigenfunctions, provide a complete description of the normal modes for a given boundary layer profile. In this study, these computations are carried out using the in-house LST solver, MAMOUT [36], which numerically solves the dispersion relation for the HBL profiles extracted from the precursor solution (with the CEDRE finite-volume solver). This is achieved through a high-order discretization scheme, as detailed in [51]. It should be mentioned that a convergence analysis with respect to the base flow resolution was conducted to assess the accuracy of our stability results.

## References

[1] Saric W S, Reed H L and Kerschen E J 2002 Boundary-layer receptivity to freestream disturbances Annu. Rev. Fluid Mech. 34 291–319

[2] Morkovin M V 1994 Transition in open flow systems a reassessment Bull. Am. Phys. Soc. 39 1882

[3] Goldstein M E 1983 The evolution of Tollmien-Schlichting waves near a leading edge J. Fluid Mech. 127 59–81

[4] Mack L M 1969 Boundary-layer stability theory (Jet Propulsion Laboratory)

[5] Mack L M 1984 Boundary-layer linear stability theory AGARD Rep. 709 1–3

[6] Özgen S and Kırcah S A 2008 Linear stability analysis in compressible, flat-plate boundary-layers, Theor. Comput. Fluid Dyn. 22 1–20

[7] Mack L M 1990 On the inviscid acoustic-mode instability of supersonic shear flows: Part 1: Two-dimensional waves, Theor. Comput. Fluid Dyn. 2 97–123

[8] Malmuth N, Fedorov A, Shalaev V, Cole J, Hites M, Williams D, and Khokhlov A 1998 Problems in high speed flow prediction relevant to control Proc. Int. Conf. 2nd AIAA, Theoretical Fluid Mechanics Meeting (Albuquerque) (AIAA Paper 2695)
$$ \mathbf{q}'(x, y, t) = \hat{\mathbf{q}}(y) \exp [i(\alpha x + \beta z - \omega t)] + c.c. \quad (A.1) $$

Here, $\alpha$ and $\beta$ represent the $x$ and $z$ components of the wavenumber vector $\mathbf{k}$, while $\omega$ is the frequency, and $\hat{\mathbf{q}}(y)$ represents eigenfunctions that describe the mode structure within the boundary layer. By introducing the normal modes into the linearized Navier-Stokes equations, a one-dimensional differential eigenvalue problem for the modal functions is obtained:

$$ \left( \mathbf{A} \frac{d}{dy^2} + \mathbf{B} \frac{d}{dy} + \mathbf{C} \right) \hat{\mathbf{q}} = 0, \quad (A.2) $$

Here, $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are complex square matrices of size 5. In the context of a boundary layer flow, the perturbation field can be described by the following boundary conditions:

$$ \text{At } y = 0, \quad \hat{u} = \hat{v} = \hat{w} = 0, \quad (A.3) $$

$$ \hat{T} = 0 \quad \text{for an isothermal wall,} \quad (A.4) $$

$$ \frac{d\hat{T}}{dy} = 0 \quad \text{for an adiabatic wall,} \quad (A.5) $$

$$ \text{At } y \to \infty, \quad \hat{u} = \hat{v} = \hat{w} = \hat{T} \to 0 \quad (A.6) $$

For spatial stability analysis, the frequency $\omega$ is real, whereas $\alpha$ is a complex eigenvalue, with $\alpha = \alpha_r + i\alpha_i$. When $\alpha_i < 0$, the flow is unstable, exhibiting a spatial growth rate $\alpha_i$ and a streamwise wavenumber $\alpha_r$. The local phase velocity is defined as $c_r = \omega/\alpha_r$. Additionally, in the context of the Mack second mode, which primarily exhibits two-dimensional dynamics, $\beta$ is assumed to be 0 in this study. The relationship for the eigenvalues, often referred to as the dispersion relation, can be expressed as:

$$ \alpha = \Omega(\omega). \quad (A.7) $$

The evaluation of the dispersion relation for a boundary layer profile is a fundamental aspect of stability theory. The eigenvalues, along with their corresponding eigenfunctions, provide a complete description of the normal modes for a given boundary layer profile. In this study, these computations are carried out using the in-house LST solver, MAMOUT [36], which numerically solves the dispersion relation for the HBL profiles extracted from the precursor solution (with the CEDRE finite-volume solver). This is achieved through a high-order discretization scheme, as detailed in [51]. It should be mentioned that a convergence analysis with respect to the base flow resolution was conducted to assess the accuracy of our stability results.

## References

[1] Saric W S, Reed H L and Kerschen E J 2002 Boundary-layer receptivity to freestream disturbances Annu. Rev. Fluid Mech. 34 291–319

[2] Morkovin M V 1994 Transition in open flow systems a reassessment Bull. Am. Phys. Soc. 39 1882

[3] Goldstein M E 1983 The evolution of Tollmien-Schlichting waves near a leading edge J. Fluid Mech. 127 59–81

[4] Mack L M 1969 Boundary-layer stability theory (Jet Propulsion Laboratory)

[5] Mack L M 1984 Boundary-layer linear stability theory AGARD Rep. 709 1–3

[6] Özgen S and Kırcah S A 2008 Linear stability analysis in compressible, flat-plate boundary-layers, Theor. Comput. Fluid Dyn. 22 1–20

[7] Mack L M 1990 On the inviscid acoustic-mode instability of supersonic shear flows: Part 1: Two-dimensional waves, Theor. Comput. Fluid Dyn. 2 97–123

[8] Malmuth N, Fedorov A, Shalaev V, Cole J, Hites M, Williams D, and Khokhlov A 1998 Problems in high speed flow prediction relevant to control Proc. Int. Conf. 2nd AIAA, Theoretical Fluid Mechanics Meeting (Albuquerque) (AIAA Paper 2695)
[9] Fedorov A 2011 Transition and stability of high-speed boundary layers *Annu. Rev. Fluid Mech.* **43** 79–95

[10] Mack L M 1975 *A numerical method for the prediction of high-speed boundary-layer transition using linear theory Proc. Int. Conf. Analyses Requiring Advanced Computers (Hampton) Part 1*

[11] Demetriades A 1974 *Hypersonic viscous flow over a slender cone. III-Laminar instability and transition Proc. Int. Conf. 7th Fluid and Plasma Dynamics Conference (Palo Alto) (AIAA Paper 535)*

[12] Egorov I V, Fedorov A V, and Soudakov V G 2006 Direct numerical simulation of disturbances generated by periodic suction-blowing in a hypersonic boundary layer *Theor. Comput. Fluid Dyn.* **20** 41–54

[13] Knisely C P and Zhong X 2019 Sound radiation by supersonic unstable modes in hypersonic blunt cone boundary layers. II. Direct numerical simulation *Phys. Fluids* **31**

[14] Fedorov A V and Khokhlov A P 2002 Receptivity of hypersonic boundary layer to wall disturbances *Theor. Comput. Fluid Dyn.* **15** 231–254

[15] Fedorov A and Tumin A 2011 High-speed boundary-layer instability: old terminology and a new framework *AIAA J.* **49** 1647–1657

[16] Fedorov A and Tumin A 2003 Initial-value problem for hypersonic boundary-layer flows *AIAA J.* **41** 379–389

[17] Ma Y and Zhong X 2003 Receptivity of a supersonic boundary layer over a flat plate. Part 1. Wave structures and interactions *J. Fluid Mech.* **488** 31–78

[18] Ma Y and Zhong X 2003 Receptivity of a supersonic boundary layer over a flat plate. Part 2. Receptivity to free-stream sound *J. Fluid Mech.* **488** 79–121

[19] Kuehl J J 2018 Thermoacoustic interpretation of second-mode instability *AIAA J.* **56** 3585–3592

[20] Tian X and Wen C 2021 Growth mechanisms of second-mode instability in hypersonic boundary layers *J. Fluid Mech.* **908** R4

[21] Chen Y, Guo P and Wen C 2023 A unified explanation of energy growth sources for unstable modes in flat-plate boundary layers *J. Fluid Mech.* **972** A5

[22] Unnikrishnan S and Gaitonde D V 2019 Interactions between vortical, acoustic and thermal components during hypersonic transition *J. Fluid Mech.* **868** 611–647

[23] Chu B-T 1965 On the energy transfer to small disturbances in fluid flow (Part I) *Acta Mech.* **1** 215–234

[24] Doak P E 1989 Momentum potential theory of energy flux carried by momentum fluctuations *J. Sound Vib.* **131** 67–90

[25] Jahanbakhshi R and Zaki T A 2021 Optimal heat flux for delaying transition to turbulence in a high-speed boundary layer *J. Fluid Mech.* **916** A46

[26] Long T, Dong Y, Zhao R and Wen C 2021 Mechanism of stabilization of porous coatings on unstable supersonic mode in hypersonic boundary layers *Phys. Fluids* **33**

[27] Atzori M, Vinuesa R, Lozano-Durán A and Schlatter P 2018 Characterization of turbulent coherent structures in square duct flow *J. Phys.: Conf. Ser.* **1001** 012008

[28] Tumuklu O, Theofilis V, and Levin D A 2018 On the unsteadiness of shock-laminar boundary layer interactions of hypersonic flows over a double cone *Phys. Fluids* **30**

[29] Chuvakhov P V and Fedorov A V 2016 Spontaneous radiation of sound by instability of a highly cooled hypersonic boundary layer *J. Fluid Mech.* **805** 188–206

[30] Lozano-Durán A and Jiménez J 2014 Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades *J. Fluid Mech.* **759** 432–471

[31] Refloch A, Courbet B, Murrone A, Villedieu P, Laurent C, Gilbank P, Troyes J, Tessé L, Chaineray G, Dargaud J B, and others 2011 CEDRE software *Aerospace Lab 2* 1–10

[32] Cassagne A, Boussuge J-F, Villedieu N, Puigt G, d'Ast I and Genot A 2015 *JAGUAR: a new CFD code dedicated to massively parallel high-order LES computations on complex geometry Proc. Int. Conf. 50th 3AF International Conference on Applied Aerodynamics (Toulouse)*

[33] Fièvet R, Deniau H, Piot E 2020 Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods *J. Comput. Phys.* **408** 109276

[34] Veilleux A, Puigt G, Deniau H, Daviller G 2022 A stable Spectral Difference approach for computations with triangular and hybrid grids up to the 6th order of accuracy *J. Comput. Phys.* **449** 110774

[35] Veilleux A, Puigt G, Deniau H, Daviller G 2022 Stable Spectral Difference Approach Using Raviart-Thomas Elements for 3D Computations on Tetrahedral Grids *J. Sci. Comput.* **91** 7

[36] Brazier J-P and Perraud J 2022 *Laminar-Turbulent Transition Modelling in High-Speed Flows Proc. Int. Conf. 2nd International Conference on High-Speed Vehicle Science Technology (Bruges)*

[37] Poinsot T J and Lele S K 1992 Boundary conditions for direct simulations of compressible viscous flows *J. Comput. Phys.* **101** 104–129

[38] Daviller G, Comte P, and Jordan P 2009 Flow decomposition for the study of source mechanisms (AIAA Paper 3305)

[39] Jordan P, Daviller G, and Comte P 2013 Doak's momentum potential theory of energy flux used to study a solenoidal wavepacket *J. Sound Vib.* **332** 3924–3936
[9] Fedorov A 2011 Transition and stability of high-speed boundary layers *Annu. Rev. Fluid Mech.* **43** 79–95

[10] Mack L M 1975 *A numerical method for the prediction of high-speed boundary-layer transition using linear theory Proc. Int. Conf. Analyses Requiring Advanced Computers (Hampton) Part 1*

[11] Demetriades A 1974 *Hypersonic viscous flow over a slender cone. III-Laminar instability and transition Proc. Int. Conf. 7th Fluid and Plasma Dynamics Conference (Palo Alto) (AIAA Paper 535)*

[12] Egorov I V, Fedorov A V, and Soudakov V G 2006 Direct numerical simulation of disturbances generated by periodic suction-blowing in a hypersonic boundary layer *Theor. Comput. Fluid Dyn.* **20** 41–54

[13] Knisely C P and Zhong X 2019 Sound radiation by supersonic unstable modes in hypersonic blunt cone boundary layers. II. Direct numerical simulation *Phys. Fluids* **31**

[14] Fedorov A V and Khokhlov A P 2002 Receptivity of hypersonic boundary layer to wall disturbances *Theor. Comput. Fluid Dyn.* **15** 231–254

[15] Fedorov A and Tumin A 2011 High-speed boundary-layer instability: old terminology and a new framework *AIAA J.* **49** 1647–1657

[16] Fedorov A and Tumin A 2003 Initial-value problem for hypersonic boundary-layer flows *AIAA J.* **41** 379–389

[17] Ma Y and Zhong X 2003 Receptivity of a supersonic boundary layer over a flat plate. Part 1. Wave structures and interactions *J. Fluid Mech.* **488** 31–78

[18] Ma Y and Zhong X 2003 Receptivity of a supersonic boundary layer over a flat plate. Part 2. Receptivity to free-stream sound *J. Fluid Mech.* **488** 79–121

[19] Kuehl J J 2018 Thermoacoustic interpretation of second-mode instability *AIAA J.* **56** 3585–3592

[20] Tian X and Wen C 2021 Growth mechanisms of second-mode instability in hypersonic boundary layers *J. Fluid Mech.* **908** R4

[21] Chen Y, Guo P and Wen C 2023 A unified explanation of energy growth sources for unstable modes in flat-plate boundary layers *J. Fluid Mech.* **972** A5

[22] Unnikrishnan S and Gaitonde D V 2019 Interactions between vortical, acoustic and thermal components during hypersonic transition *J. Fluid Mech.* **868** 611–647

[23] Chu B-T 1965 On the energy transfer to small disturbances in fluid flow (Part I) *Acta Mech.* **1** 215–234

[24] Doak P E 1989 Momentum potential theory of energy flux carried by momentum fluctuations *J. Sound Vib.* **131** 67–90

[25] Jahanbakhshi R and Zaki T A 2021 Optimal heat flux for delaying transition to turbulence in a high-speed boundary layer *J. Fluid Mech.* **916** A46

[26] Long T, Dong Y, Zhao R and Wen C 2021 Mechanism of stabilization of porous coatings on unstable supersonic mode in hypersonic boundary layers *Phys. Fluids* **33**

[27] Atzori M, Vinuesa R, Lozano-Durán A and Schlatter P 2018 Characterization of turbulent coherent structures in square duct flow *J. Phys.: Conf. Ser.* **1001** 012008

[28] Tumuklu O, Theofilis V, and Levin D A 2018 On the unsteadiness of shock-laminar boundary layer interactions of hypersonic flows over a double cone *Phys. Fluids* **30**

[29] Chuvakhov P V and Fedorov A V 2016 Spontaneous radiation of sound by instability of a highly cooled hypersonic boundary layer *J. Fluid Mech.* **805** 188–206

[30] Lozano-Durán A and Jiménez J 2014 Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades *J. Fluid Mech.* **759** 432–471

[31] Refloch A, Courbet B, Murrone A, Villedieu P, Laurent C, Gilbank P, Troyes J, Tessé L, Chaineray G, Dargaud J B, and others 2011 CEDRE software *Aerospace Lab 2* 1–10

[32] Cassagne A, Boussuge J-F, Villedieu N, Puigt G, d'Ast I and Genot A 2015 *JAGUAR: a new CFD code dedicated to massively parallel high-order LES computations on complex geometry Proc. Int. Conf. 50th 3AF International Conference on Applied Aerodynamics (Toulouse)*

[33] Fièvet R, Deniau H, Piot E 2020 Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods *J. Comput. Phys.* **408** 109276

[34] Veilleux A, Puigt G, Deniau H, Daviller G 2022 A stable Spectral Difference approach for computations with triangular and hybrid grids up to the 6th order of accuracy *J. Comput. Phys.* **449** 110774

[35] Veilleux A, Puigt G, Deniau H, Daviller G 2022 Stable Spectral Difference Approach Using Raviart-Thomas Elements for 3D Computations on Tetrahedral Grids *J. Sci. Comput.* **91** 7

[36] Brazier J-P and Perraud J 2022 *Laminar-Turbulent Transition Modelling in High-Speed Flows Proc. Int. Conf. 2nd International Conference on High-Speed Vehicle Science Technology (Bruges)*

[37] Poinsot T J and Lele S K 1992 Boundary conditions for direct simulations of compressible viscous flows *J. Comput. Phys.* **101** 104–129

[38] Daviller G, Comte P, and Jordan P 2009 Flow decomposition for the study of source mechanisms (AIAA Paper 3305)

[39] Jordan P, Daviller G, and Comte P 2013 Doak's momentum potential theory of energy flux used to study a solenoidal wavepacket *J. Sound Vib.* **332** 3924–3936
[40] Moisy F and Jiménez J 2004 Geometry and clustering of intense structures in isotropic turbulence *J. Fluid Mech.* **513** 111–133

[41] Del Alamo JC, Jiménez J, Zandonade P and Moser RD 2006 Self-similar vortex clusters in the turbulent logarithmic region *J. Fluid Mech.* **561** 329–358

[42] Hao Z, Yan C, Zhou L and Qin Y 2017 Development of a boundary layer parameters identification method for transition prediction with complex grids *J. Aerosp. Eng.* **231** 2068–2084

[43] Unnikrishnan S and Gaitonde DV 2021 Instabilities and transition in cooled wall hypersonic boundary layers *J. Fluid Mech.* **915** A26

[44] Zhang S, Liu J, and Luo J 2016 Effect of wall-cooling on Mack-mode instability in high speed flat-plate boundary layers *Appl. Math. Mech.*, **37**, 1219–1230

[45] Fedorov AV and Khokhlov AP 2001 Prehistory of instability in a hypersonic boundary layer *Theor. Comput. Fluid Dyn.* **14** 359–375

[46] Ma Y and Zhong X 2005 Receptivity of a supersonic boundary layer over a flat plate. Part 3. Effects of different types of free-stream disturbances *J. Fluid Mech.* **532** 63–109

[47] Mack L M 1975 Linear stability theory and the problem of supersonic boundary-layer transition AIAA J. **13** 278–289

[48] Stetson K, Thompson E, Donaldson J and Siler L 1989 *Laminar boundary layer stability experiments on a cone at Mach 8. V-Tests with a cooled model, Proc. Int. Conf. 20th Fluid Dynamics, Plasma Dynamics and Lasers Conference (Buffalo)* (AIAA Paper 6.1989-1895)

[49] Osawa K and Jiménez J 2018 Intense structures of different momentum fluxes in turbulent channels *Phys. Rev. Fluids* **3** 084603

[50] Wang X, Zhong X and Ma Y 2011 Response of a hypersonic boundary layer to wall blowing-suction AIAA J. **49** 1336–1353

[51] Malik M R 1990 Numerical methods for hypersonic boundary layer stability *J. Comput. Phys.* **86** 376–413
[40] Moisy F and Jiménez J 2004 Geometry and clustering of intense structures in isotropic turbulence *J. Fluid Mech.* **513** 111–133

[41] Del Alamo JC, Jiménez J, Zandonade P and Moser RD 2006 Self-similar vortex clusters in the turbulent logarithmic region *J. Fluid Mech.* **561** 329–358

[42] Hao Z, Yan C, Zhou L and Qin Y 2017 Development of a boundary layer parameters identification method for transition prediction with complex grids *J. Aerosp. Eng.* **231** 2068–2084

[43] Unnikrishnan S and Gaitonde DV 2021 Instabilities and transition in cooled wall hypersonic boundary layers *J. Fluid Mech.* **915** A26

[44] Zhang S, Liu J, and Luo J 2016 Effect of wall-cooling on Mack-mode instability in high speed flat-plate boundary layers *Appl. Math. Mech.*, **37**, 1219–1230

[45] Fedorov AV and Khokhlov AP 2001 Prehistory of instability in a hypersonic boundary layer *Theor. Comput. Fluid Dyn.* **14** 359–375

[46] Ma Y and Zhong X 2005 Receptivity of a supersonic boundary layer over a flat plate. Part 3. Effects of different types of free-stream disturbances *J. Fluid Mech.* **532** 63–109

[47] Mack L M 1975 Linear stability theory and the problem of supersonic boundary-layer transition AIAA J. **13** 278–289

[48] Stetson K, Thompson E, Donaldson J and Siler L 1989 *Laminar boundary layer stability experiments on a cone at Mach 8. V-Tests with a cooled model, Proc. Int. Conf. 20th Fluid Dynamics, Plasma Dynamics and Lasers Conference (Buffalo)* (AIAA Paper 6.1989-1895)

[49] Osawa K and Jiménez J 2018 Intense structures of different momentum fluxes in turbulent channels *Phys. Rev. Fluids* **3** 084603

[50] Wang X, Zhong X and Ma Y 2011 Response of a hypersonic boundary layer to wall blowing-suction AIAA J. **49** 1336–1353

[51] Malik M R 1990 Numerical methods for hypersonic boundary layer stability *J. Comput. Phys.* **86** 376–413