# Linear stability analysis of compressible boundary layer over an insulated
wall: Existence of multiple new unstable modes for Mach number beyond 3

Neha Chaturvedi1<EMAIL_ADDRESS><EMAIL_ADDRESS>Swagata Bhaumik2<EMAIL_ADDRESS>Rituparn Somvanshi3<EMAIL_ADDRESS>1,2,3Department Of
Mechanical Engineering, IIT (ISM) Dhanbad, Dhanbad - 826004, India

###### Abstract

We investigate the linear stability of parallel two-dimensional ($2D$)
compressible boundary layer flow over a smooth, adiabatic flat plate. We
consider both two- and three-dimensional ($3D$) disturbances, which yield
$6^{th}$\- and $8^{th}$-order systems for wall-normal variation of spectral
amplitudes, respectively. We perform spatial stability analysis of the flow
using the compound matrix method (CMM) to remove the stiffness of the problem,
unlike the conventional techniques using Gram-Schmidt ortho-normalization or
discretizing the governing equations by appropriate finite difference schemes.
This method has yet to be employed for linear stability analysis of
compressible boundary layers. We consider the Mach number $M$ of the flow to
vary from low subsonic cases ($M=0.1$) to supersonic cases ($M=6$). First, we
validate the compound matrix method by comparing the $M=0.1$ case results with
the incompressible boundary layer governed by the Orr-Sommerfeld equation. The
results depict an excellent match for primary and secondary modes. The flow
becomes increasingly stable with an increase in $M$ up to $M=2$. We further
note subsonic cases to be more stable when considering $3D$ disturbances. On
the contrary, this conclusion holds only for $3D$ disturbances having spanwise
wavenumber $\beta$ greater than a specific optimum value corresponding to
supersonic cases up to $M=2$. We show that $2D$ disturbances are more
spatially unstable than $3D$ perturbations with no spanwise variation (i.e.,
$\beta=0$) up to $M=1.4$, and the opposite scenario happens for $1.4<M<2$.
Only one primary spatially unstable mode is found for boundary layer cases up
to $M=2$. We locate a series of unstable modes for $M>3$, and the number of
such modes is much more than two, as reported in Mack (AGARD Report No. 709,
Part 3, 1984) from viscous calculations. Mack reports only two unstable modes
up to $M=4.5$ that subsequently fuse with an increase in $M$. Our results show
the number and the frequency extent of the unstable zone for these modes
increase significantly with an increase in Mach number, e.g., we find $6$,
$8$, and $10$ spatially unstable modes up to streamwise wavenumber of
$\alpha_{r}=0.4$ for $M=3$, $4$ and $6$, respectively. The calculation of
group velocity for these modes shows that these propagate downstream at a
higher speed of $0.78$ to $0.89$ times the free-stream speed than those
corresponding to incompressible, subsonic, and low supersonic ($M<2$) cases.
While the shape of the neutral curves for the second unstable mode for $M=4$
and $6$ is similar to the fused neutral curve shown in Mack for $M=4.8$, the
characteristics of spatially unstable higher-order modes, to the best of our
knowledge, are not been reported elsewhere so far regarding the viscous
stability of supersonic boundary layer.

###### keywords:

Boundary Layer, Linear stability analysis, Compressible flow, Dispersion
relation, Normal mode analysis, Compound Matrix Method

##  1 Introduction

The stability of high-speed boundary layer and associated transition
prediction is of great practical importance in designing subsonic and
supersonic aircraft, hypersonic re-entry vehicles, gas turbine blades, and
vanes [1, 2]. This information is essential to predict the above objects’ skin
friction, overall drag, and heat transfer characteristics. Significant
theoretical [3, 4, 5], experimental [6, 7], and numerical [8, 9, 10, 11, 12]
advances have been made to illustrate the transition to turbulence for high
speed flows over the last several decades. A comprehensive account of various
routes in this regard is given in Reshotko [13] and Lee [2].

The transition of boundary layer flows from laminar to the turbulent stage is
conventionally attributed to the growth of instability waves describing the
evolution of small disturbances over the base flow [14, 15]. Initial
pioneering efforts to determine the stability of the incompressible boundary
layer are due to Tollmein [14, 16], Schlichting [17], and Heisenberg [18], who
obtained the analytical solution of the classical Orr-Sommerfeld equation
(OSE) that governs the linear evolution of small perturbations in viscous
incompressible flow framework, for a base flow that does not vary along
streamwise direction (parallel flow approximation). These initial results
demonstrated that complex interactions of viscous terms with the inertial
terms of the governing equation can result in the amplification of
monochromatic waves for flows that are stable according to the inviscid
instability theories of Kelvin and Rayleigh [19, 20]. The growth of these
viscous tuned monochromatic waves (later termed as Tollmien-Schlichting or
$TS$-waves) in space/time is subsequently thought to be responsible for flow
transition from laminar to turbulent state, ever since Schubauer & Skramstad
[21] reported corresponding experimental detections under controlled
environment.

The boundary layer transition process is more complicated for compressible
flows or flows with heat transfer than the incompressible flows [10, 11].
Spatial or temporal linear stability analysis of transonic and supersonic
flows has been reported in Mack [22, 3] and Reed [4]. These studies are
obtained by normal mode analysis in the form of temporally or spatially
unstable disturbances. Mack [3] noted the existence of two unstable modes from
$M=3$ to approximately $4.5$ for the viscous stability calculations of
compressible flows over an insulated wall. The second mode corresponds to high
wavenumber disturbances; therefore, its neutral curve is observed to lie on
top of the first mode. Two neutral curves are noted to merge with an increase
in Mach number to $4.8$ and beyond, denoting effectively the existence of one
unstable zone whose spectral range is much wider than those corresponding to
subsonic or low supersonic cases [23]. The existence of two separate unstable
modes observed at intermediate supersonic Mach numbers is also noted for high-
speed flows over hot and cold isothermal wall cases [24], which subsequently,
with an increase in flow Mach number, fuses to form a much wider single
unstable zone. Experimental investigations on high-speed boundary layer
transition confirm the existence and dominance of high-frequency (or high
wavenumber) second mode [25, 26, 27, 28]. Further details are also given in
Fedorov & Tumin [29] and Fedorov et. al. [30].

Conventionally, the stability calculations are performed by two approaches
once the governing ordinary differential equation (ODE) for the spectral
perturbation amplitudes under parallel flow approximation is obtained. The
problem is cast as essentially an eigenvalue problem in the first approach,
corresponding to the matrix obtained by discretizing the governing equations
along the wall-normal directions by suitable numerical finite-difference
schemes or using Chebyshev polynomials [31, 32, 8]. This method is more
appropriate for determining the temporal instability of the boundary layer as
the streamwise and spanwise wavenumbers appear nonlinearly in these equations.
Moreover, most of the eigenvalues provided by this method are numerical and
spurious, whose number depends upon the number of discretization points and
equations involved. Screening out the few real physical eigenmodes selectively
from this pool of spurious ones is challenging. An alternative second approach
for stability calculations of the viscous boundary layer is illustrated in
Mack[33, 3], which does not yield numerous spurious modes and is suitable for
analysis of spatial, temporal, or even spatiotemporal instability of the
boundary layer. Following this method, integration of spectral amplitude
equations by suitable numerical schemes like Runge-Kutta or Runge-Kutta-
Fehlberg methods [24, 34] is carried out from the known initial conditions at
the free-stream. Straightforward numerical integration is impossible as the
resulting coupled ODEs are inherently stiff, amplifying the spurious modes.
Therefore, Gram-Schmidt ortho-normalization is applied at every step to
circumvent this issue[33, 24, 34]. Subsequently, the eigenvalues are obtained
by satisfying the appropriate dispersion relation at the wall.

Ng & Reid [35], Yiantsios & Higgins [36] and Allen & Bridges [37] have
proposed the compound matrix method (CMM) to remove the associated
characteristic stiffness of such ODEs so that one can perform straightforward
numerical integration without the solution diverging due to the appearance of
growing non-physical modes. Compared to the method proposed by Mack[33], here,
no orthonormalization is required; instead, the integration is carried over
the compound variables, which are essentially the minors of the modes
representing the physically realizable solutions. This method has been
successfully applied for the stability calculations of incompressible wall-
bounded [15, 38] and mixed-convection boundary layer [39] problems. CMM has
also found applications in solving stiff governing equations related to the
stability of plane Poiseuille flow containing multiple intermediate interfaces
[40], nonlinear traveling waves [41], and pre-stressed elastic tube under
axial compression [42]. This method has so far not been used for stability
calculations of the compressible boundary layers for which the constitution of
corresponding compound matrix equations are more complex than incompressible
hydrodynamic or mixed-convection boundary layer, as illustrated later in Sec.
3. It is performed here.

Results show that while the instability characteristic matches existing
results for subsonic and low supersonic cases, we report the existence of
multiple unstable modes for $M>3$ with increasing wavenumber. The number of
such unstable modes increases with an increase in Mach number beyond $3$. In
contrast, Mack [3] or subsequent similar investigations [24, 34] report at
most two distinct unstable modes at intermediate supersonic Mach number cases.
Our calculations show that the extent of the corresponding unstable spectral
range is much broader than subsonic or low supersonic cases, which is enhanced
with an increase in $M$. Several such unstable modes are characterized and
illustrated in Sec. 5.4. It is the novel element in the present manuscript.
The existence of multiple (much more than two) unstable viscous modes has not
been reported before for high-speed flows.

The paper is organized as follows. Section 2 describes linearized disturbance
equations, equations for respective spectral amplitudes and the far-field
variation of disturbances corresponding to compressible boundary layers.
Sections 3 and 4 illustrate the dispersion function for determining the
eigenmodes of compressible shear-layer over the insulated wall and wall-normal
variation of the mean flow as obtained from corresponding self-similar
equations[43], respectively. Section 5 provides the results of the instability
analysis as obtained from CMM for the compressible boundary layer over the
adiabatic wall. We start with validating the CMM methodology by comparing the
results for $M=0.1$ with corresponding incompressible results in Sec. 5.1.
Subsequently, we explore the effects of variation of Mach number and spanwise
wavenumber ($3D$-effect) on the flow instability in Secs. 5.2 and 5.3,
respectively. Results for high supersonic Mach numbers ($M=3$, $4$, and $6$)
are provided in Sec. 5.4, which describes the characteristics of multiple
unstable modes for these flow cases. The summary and conclusion are provided
in Sec. 6.

##  2 linearized disturbance equations for compressible boundary layer

The dimensional form of the governing unsteady Navier-Stokes equations (NSE)
for the compressible flow, in vectorial notations are given as

$\displaystyle\frac{\partial\tilde{\rho}}{\partial\tilde{t}}+\tilde{\nabla}\cdot\left(\tilde{\rho}\tilde{\bf
v}\right)=0$ (1) $\displaystyle\frac{\partial}{\partial
t}\left(\tilde{\rho}\tilde{\bf
v}\right)+\tilde{\nabla}\cdot\left(\tilde{\rho}\tilde{\bf v}\tilde{\bf
v}\right)=-\tilde{\nabla}\tilde{p}+\tilde{\nabla}\cdot\tilde{\bf\tau}$ (2)
$\displaystyle\frac{\partial}{\partial
t}\left(\tilde{\rho}\tilde{e}_{t}\right)+\tilde{\nabla}\cdot\left(\tilde{\rho}\tilde{\bf
v}\tilde{h}_{t}\right)=\tilde{\nabla}\cdot\left(\tilde{\bf\tau}\cdot\tilde{\bf
v}\right)-\tilde{\nabla}\cdot\tilde{\bf q}$ (3)

where $\tilde{\rho}$ represents fluid density; $\tilde{T}$ represents fluid
temperature; $p$ represents thermodynamic pressure; $\tilde{\bf v}$ represents
fluid velocity; $\tilde{e}_{i}$ and $\tilde{e}_{t}=\frac{1}{2}\left(\tilde{\bf
v}\cdot\tilde{\bf v}\right)+\tilde{e}_{i}$ define the specific internal and
total energy, respectively; $h_{i}=\tilde{e}_{i}+\tilde{p}/\tilde{\rho}$ and
$\tilde{h}_{t}=\tilde{e}_{t}+\tilde{p}/\tilde{\rho}$ represent specific
internal and total enthalpy, respectively; $\tilde{\bf q}$ represents the
heat-flux and $\tilde{\bf\tau}$ represents the viscous stress tensor. For an
isotropic and Newtonian fluid obeying Fourier law of heat conduction, the
symmetric stress tensor and the heat-flux vector are given as
$\tilde{\bf\tau}=-\tilde{\lambda}D{\bf I}+2\tilde{\mu}\tilde{\bf\epsilon}$ and
$\tilde{\bf q}=-\tilde{\kappa}\tilde{\nabla}\tilde{T}$, respectively [44].
Here, $\tilde{\bf\epsilon}=\frac{1}{2}\left(\tilde{\nabla}\tilde{\bf
v}+\tilde{\nabla}\tilde{\bf v}^{T}\right)$, and
$D=\tilde{\nabla}\cdot\tilde{\bf v}$ are the symmetric strain-rate tensor and
the volumetric dilatation rate, respectively. The variables, $\tilde{\mu}$,
$\tilde{\lambda}$ and $\tilde{\kappa}$ denote the dynamic and the second
coefficient of viscosity and the thermal conductivity of the fluid,
respectively. According to Stokes’ hypothesis, [44]
$\tilde{\lambda}=2\tilde{\mu}/3$, which indicates that the bulk viscosity
coefficient of the flow $k=0$. Here, we assume a calorically perfect gas, and
therefore, the equation of state is given as
$\tilde{p}=R\;\tilde{\rho}\;\tilde{T}$ so that $\tilde{e}_{i}=c_{v}\tilde{T}$
and $h_{i}=c_{p}\tilde{T}$ where, $R$, $c_{v}=R/(\gamma-1)$, $c_{p}=\gamma
R/(\gamma-1)$, and $\gamma=c_{p}/c_{v}$ are the universal gas constant,
specific heat at constant volume and constant pressure and ratio of specific
heats, respectively. We also assume here that $\tilde{\mu}$,
$\tilde{\lambda}$, and $\tilde{\kappa}$ are functions of temperature
$\tilde{T}$ alone.

We use the free-stream values for velocity $\tilde{U}_{\infty}$, temperature
$\tilde{T}_{\infty}$, density $\tilde{\rho}_{\infty}$, pressure
$\tilde{p}_{\infty}=R\tilde{\rho}_{\infty}\tilde{T}_{\infty}$, viscosity
$\tilde{\mu}_{\infty}$, heat-conductivity $\tilde{\kappa}_{\infty}$ as the
corresponding reference values to non-dimensionalize the above equations. The
non-dimensional Reynolds number, Mach number and Prandtl number of the flow is
given as $Re=\tilde{\rho}_{\infty}\tilde{U}_{\infty}L/\tilde{\mu}_{\infty}$,
$M=\tilde{U}_{\infty}/\sqrt{\gamma R\tilde{T}_{\infty}}$ and
$Pr=c_{p}\mu_{\infty}/k_{\infty}$, respectively, where $L$ is the reference
length-scale.

Figure 1: The schematic of the $2D$ parallel flow-approximation showing a
tentative variation of the velocity and temperature profiles of a wall-bounded
shear layer. Here, $\delta_{*}$ represents the displacement thickness of the
shear-layer, which is treated as constant here under the parallel flow
approximation.

Next, to consider the disturbance evolution, we split each flow variable
$\phi$ into a mean and fluctuating component as
$\phi=\bar{\phi}+\epsilon\hat{\phi}$, where $\epsilon<<1$. Therefore, we
assume the fluctuating quantities are much smaller than the mean flow
variables. The mean quantities are functions of spatial coordinates only,
while the fluctuating quantities depend on both space and time. By
substituting the above decomposition into Eqs. (1-3), we obtain the linearized
equations for the disturbance quantities. We denote $x$, $y$, and $z$ as the
non-dimensional streamwise, wall-normal, and spanwise coordinates,
respectively. We consider a $2D$ parallel mean flow (as shown in Fig. 1), and
therefore, the mean quantities, i.e., $\bar{\varrho}$, $\bar{P}$, $\bar{T}$
and $\bar{U}$ are only functions of the wall-normal coordinate $y$ only. Here,
$\bar{U}$ is the streamwise component of the mean flow. As the boundary layer
is assumed to be $2D$ and parallel, the wall-normal and spanwise component of
the mean flow is assumed to be zero, i.e., $\bar{V}=\bar{W}=0$. We consider
the local displacement thickness $\delta_{*}$ as the reference length-scale,
i.e., $L=\delta_{*}$. Let, $\hat{\textbf{v}}=(\hat{u},\hat{v},\hat{w})$
represents perturbation velocity components, and $\hat{p}$, $\hat{\theta}$ and
$\hat{\rho}$ denote disturbance pressure, temperature, and density,
respectively. The non-dimensional linearized disturbance equations under $2D$
parallel flow approximation are given as [3, 34]

$\displaystyle\frac{D\hat{\rho}}{Dt}+\hat{v}\frac{d\bar{\varrho}}{dy}+\bar{\varrho}\hat{D}$
$\displaystyle=$ $\displaystyle 0$ (4)
$\displaystyle\bar{\varrho}\frac{D\hat{u}}{Dt}+\bar{\varrho}\hat{v}\frac{d\bar{U}}{dy}+\frac{1}{\gamma
M^{2}}\frac{\partial\hat{p}}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{1}{Re}\left(\frac{\partial\hat{\tau}_{xx}}{\partial
x}+\frac{\partial\hat{\tau}_{xy}}{\partial
y}+\frac{\partial\hat{\tau}_{xz}}{\partial z}\right)$ (5)
$\displaystyle\bar{\varrho}\frac{D\hat{v}}{Dt}+\frac{1}{\gamma
M^{2}}\frac{\partial\hat{p}}{\partial y}$ $\displaystyle=$
$\displaystyle\frac{1}{Re}\left(\frac{\partial\hat{\tau}_{yx}}{\partial
x}+\frac{\partial\hat{\tau}_{yy}}{\partial
y}+\frac{\partial\hat{\tau}_{yz}}{\partial y}\right)$ (6)
$\displaystyle\bar{\varrho}\frac{D\hat{w}}{Dt}+\frac{1}{\gamma
M^{2}}\frac{\partial\hat{p}}{\partial z}$ $\displaystyle=$
$\displaystyle\frac{1}{Re}\left(\frac{\partial\hat{\tau}_{zx}}{\partial
x}+\frac{\partial\hat{\tau}_{zy}}{\partial
y}+\frac{\partial\hat{\tau}_{zz}}{\partial z}\right)$ (7)
$\displaystyle\bar{\varrho}\bar{c}_{p}\frac{D\hat{\theta}}{Dt}-\frac{\gamma-1}{\gamma}\frac{D\hat{p}}{Dt}+\bar{\varrho}\bar{c}_{p}\hat{v}\frac{d\bar{T}}{dy}$
$\displaystyle=$
$\displaystyle\frac{1}{Re}\frac{\gamma-1}{M^{2}}\left[\hat{\tau}_{xy}\frac{d\bar{U}}{dy}+\bar{\mu}\frac{d\bar{U}}{dy}\left(\frac{\partial\hat{u}}{\partial
y}+\frac{\partial\hat{v}}{\partial x}\right)\right]$ (8)
$\displaystyle-\frac{1}{Re}\frac{1}{Pr}\left(\frac{\partial\hat{q}_{x}}{\partial
x}+\frac{\partial\hat{q}_{y}}{\partial y}+\frac{\partial\hat{q}_{z}}{\partial
z}\right)$

where, the operator $\frac{D}{Dt}=\left(\frac{\partial}{\partial
t}+\bar{U}(y)\frac{\partial}{\partial x}\right)$. The perturbation stress
tensor components are given as
$\hat{\tau}_{xx}=2\bar{\mu}\frac{\partial\hat{u}}{\partial
x}+\bar{\lambda}\hat{D}$,
$\hat{\tau}_{yx}=\hat{\tau}_{xy}=\bar{\mu}\left(\frac{\partial\hat{u}}{\partial
y}+\frac{\partial\hat{v}}{\partial
x}\right)+\frac{d\bar{\mu}}{d\bar{T}}\frac{d\bar{U}}{dy}\hat{\theta}$,
$\hat{\tau}_{xz}=\hat{\tau}_{zx}=\bar{\mu}\left(\frac{\partial\hat{u}}{\partial
z}+\frac{\partial\hat{w}}{\partial x}\right)$,
$\hat{\tau}_{yy}=2\bar{\mu}\frac{\partial\hat{v}}{\partial
y}+\bar{\lambda}\hat{D}$,
$\hat{\tau}_{yz}=\hat{\tau}_{zy}=\bar{\mu}\left(\frac{\partial\hat{w}}{\partial
y}+\frac{\partial\hat{v}}{\partial z}\right)$, and
$\hat{\tau}_{zz}=2\bar{\mu}\frac{\partial\hat{w}}{\partial
z}+\bar{\lambda}\hat{D}$ where $\hat{D}=\left(\frac{\partial\hat{u}}{\partial
x}+\frac{\partial\hat{v}}{\partial y}+\frac{\partial\hat{w}}{\partial
z}\right)$. Components of perturbation heat-flux vector are given as
$\hat{q}_{x}=-\bar{\kappa}\frac{\partial\hat{\theta}}{\partial x}$,
$\hat{q}_{y}=-\left(\bar{\kappa}\frac{\partial\hat{\theta}}{\partial
y}+\frac{d\bar{\kappa}}{d\bar{T}}\frac{d\bar{T}}{dy}\hat{\theta}\right)$ and
$\hat{q}_{z}=-\bar{\kappa}\frac{\partial\hat{\theta}}{\partial z}$. The
relationship between $\hat{p}$, $\hat{\theta}$ and $\hat{\rho}$ are derived
from the perturbation equation of state given as
$\hat{p}=\left(\hat{\rho}\bar{T}+\bar{\varrho}\hat{\theta}\right)$.

###  2.1 Variation of spectral amplitudes of the linearized disturbances

Equations (4-8) are subsequently expressed in the spectral domain by using the
Fourier-Laplace transform of the perturbation quantities as

$\displaystyle\left(\hat{u},\hat{v},\hat{w},\hat{p},\hat{\rho},\hat{\theta}\right)^{T}=\int\int\left(\chi(y),\varphi(y),\Upsilon(y),\Pi(y),\zeta(y),\Theta(y)\right)^{T}e^{i\left(\alpha
x+\beta z-\omega t\right)}d\alpha d\beta$ (9)

where $\alpha$ and $\beta$ represent streamwise and spanwise wavenumber,
respectively, while $\omega$ denote the circular frequency of the perturbation
components. The variables $\chi$, $\varphi$, $\Upsilon$, $\Pi$, $\zeta$, and
$\Theta$ indicate the spectral amplitudes of disturbance quantities $\hat{u}$,
$\hat{v}$, $\hat{w}$, $\hat{p}$, $\hat{\rho}$, and $\hat{\theta}$,
respectively. We define two auxiliary variables
$\psi=\left(\alpha\chi+\beta\Upsilon\right)$ and
$\Omega=\left(\beta\chi-\alpha\Upsilon\right)$. The variable $\psi$ represents
the spectral amplitude of $\hat{D_{1}}=\left(\frac{\partial\hat{u}}{\partial
x}+\frac{\partial\hat{w}}{\partial z}\right)$, while $\Omega$ denotes the
spectral amplitude of the perturbation wall-normal component of vorticity
$\hat{\xi}_{y}=\left(\frac{\partial\hat{u}}{\partial
z}-\frac{\partial\hat{w}}{\partial z}\right)$. We note that in the context of
the stability of incompressible viscous flows, the Squire equation is
essentially the wall-normal variation of the perturbation wall-normal
vorticity component $\hat{\xi}_{y}$ [45].

Substituting the Fourier-Laplace transform of perturbations given by Eqs (9)
into the linearized disturbance evolution equations (4-8), we obtain the wall-
normal variation of respective spectral amplitudes as [3, 34]

$\displaystyle
iQ\zeta+\varphi\left(\frac{d\bar{\varrho}}{dy}\right)+\bar{\varrho}\left(\varphi^{\prime}+i\psi\right)$
$\displaystyle=$ $\displaystyle 0$ (10) $\displaystyle
i\bar{\varrho}Q\chi+\varphi\left(\bar{\varrho}\frac{d\bar{U}}{dy}\right)+\frac{1}{\gamma
M^{2}}i\alpha\Pi$ $\displaystyle=$
$\displaystyle\frac{\bar{\mu}}{Re}\left[\chi^{\prime\prime}-\Delta^{2}\chi+i\alpha\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\left(\varphi^{\prime}+i\psi\right)\right]$
(11)
$\displaystyle+\frac{1}{Re}\biggl{[}C_{1}\left(i\alpha\varphi+\chi^{\prime}\right)+C_{2}\Theta+C_{3}\Theta^{\prime}\biggl{]}$
$\displaystyle i\bar{\varrho}Q\varphi+\frac{1}{\gamma M^{2}}\Pi^{\prime}$
$\displaystyle=$
$\displaystyle\frac{\bar{\mu}}{Re}\left[\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)\varphi^{\prime\prime}-\Delta^{2}\varphi+i\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\psi^{\prime}\right]+$
(12)
$\displaystyle\frac{1}{Re}\biggl{[}i\frac{\bar{\lambda}}{\bar{\mu}}C_{1}\psi+\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)C_{1}\varphi^{\prime}+i\alpha
C_{3}\Theta\biggl{]}$ $\displaystyle i\bar{\varrho}Q\Upsilon+\frac{1}{\gamma
M^{2}}i\beta\Pi$ $\displaystyle=$
$\displaystyle\frac{\bar{\mu}}{Re}\left[\Upsilon^{\prime\prime}-\Delta^{2}\Upsilon+i\beta\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\left(\varphi^{\prime}+i\psi\right)\right]$
(13)
$\displaystyle+\frac{C_{1}}{Re}\biggl{[}i\beta\varphi+\Upsilon^{\prime}\biggl{]}$
$\displaystyle
iQ\left(\bar{\varrho}c_{p}\Theta-\frac{\gamma-1}{\gamma}\Pi\right)+\left(\bar{\varrho}c_{p}\frac{d\bar{T}}{dy}\right)\varphi$
$\displaystyle=$
$\displaystyle\frac{1}{Re}\frac{1}{Pr}\biggl{[}\bar{\kappa}\left(\Theta^{\prime\prime}-\Delta^{2}\Theta\right)+C_{9}\Theta+2C_{10}\Theta^{\prime}\biggl{]}+$
(14)
$\displaystyle\frac{(\gamma-1)M^{2}}{Re}\left[C_{8}\Theta+2\bar{\mu}\frac{d\bar{U}}{dy}\left(\chi^{\prime}+i\alpha\varphi\right)\right]$
$\displaystyle\Pi=\left(\bar{T}\zeta+\bar{\varrho}\Theta\right)$ (15)

where, $(\cdot)^{\prime}=d(\cdot)/dy$, $\Delta=\sqrt{\alpha^{2}+\beta^{2}}$,
and $Q=\left(\alpha\bar{U}-\omega\right)$. Other variables appearing in Eqs.
(10-14) are given as $C_{1}=\frac{d{\mu}}{d\bar{T}}\frac{d\bar{T}}{dy}$,
$C_{2}=\left(\frac{d^{2}{\mu}}{d\bar{T}^{2}}\frac{d\bar{T}}{dy}\frac{d\bar{U}}{dy}+\frac{d{\mu}}{d\bar{T}}\frac{d^{2}\bar{U}}{dy^{2}}\right)$,
$C_{3}=\left(\frac{d{\mu}}{d\bar{T}}\frac{d\bar{U}}{dy}\right)$,
$C_{8}=\left(\frac{d{\mu}}{d\bar{T}}\left(\frac{d\bar{U}}{dy}\right)^{2}\right)$,
$C_{9}=\left(\frac{d^{2}{\kappa}}{d\bar{T}^{2}}\left(\frac{d\bar{T}}{dy}\right)^{2}+\frac{d{\kappa}}{d\bar{T}}\frac{d^{2}\bar{T}}{dy^{2}}\right)$,
and $C_{10}=\left(\frac{d{\kappa}}{d\bar{T}}\frac{d\bar{T}}{dy}\right)$.
Substituting Eq. (15) in Eqs. (10-14), and further simplifying, one gets
equations in terms of $\psi$, $\varphi$, $\Theta$ and $\Omega$ as

$\displaystyle\left[\frac{\bar{\mu}}{Re}\right]\psi^{\prime\prime}+\left[\frac{C_{1}}{Re}\right]\psi^{\prime}+\left[i\frac{\bar{\mu}}{Re}\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\Delta^{2}+\frac{F_{1}}{\gamma
M^{2}}\Delta^{2}\right]\varphi^{\prime}+\left[\frac{\alpha
C_{3}}{Re}\right]\Theta^{\prime}$ $\displaystyle+\left[i\frac{F_{1}}{\gamma
M^{2}}\Delta^{2}-\frac{\bar{\mu}}{Re}\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)\Delta^{2}-i\bar{\varrho}Q\right]\psi+\left[\frac{F_{1}}{\gamma
M^{2}}\frac{1}{\bar{\varrho}}\frac{d\bar{\varrho}}{dy}\Delta^{2}-\alpha\bar{\varrho}\frac{d\bar{U}}{dy}+i\frac{C_{1}}{Re}\Delta^{2}\right]\varphi$
$\displaystyle+\left[\frac{\alpha C_{2}}{Re}-i\frac{\bar{\varrho}}{\gamma
M^{2}}\Delta^{2}\right]\Theta=0$ (16)
$\displaystyle\left[\frac{\bar{\mu}}{Re}\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)Q^{2}-i\frac{\bar{\varrho}\bar{T}}{\gamma
M^{2}}Q\right]\varphi^{\prime\prime}+\left[i\frac{\bar{\mu}}{Re}\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)Q^{2}+\frac{\bar{\varrho}\bar{T}}{\gamma
M^{2}}Q\right]\psi^{\prime}$
$\displaystyle+\left[\frac{C_{1}}{Re}\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)Q^{2}-i\frac{G_{2}}{\gamma
M^{2}}\right]\varphi^{\prime}+\left[-\frac{\bar{\varrho}}{\gamma
M^{2}}Q^{2}\right]\Theta^{\prime}+\left[i\frac{C_{1}}{Re}\left(\frac{\bar{\lambda}}{\bar{\mu}}\right)Q^{2}+\frac{G_{3}}{\gamma
M^{2}}\right]\psi$
$\displaystyle+\left[-\frac{\bar{\mu}}{Re}Q^{2}\Delta^{2}-i\bar{\varrho}Q^{3}-i\frac{G_{4}}{\gamma
M^{2}}\right]\varphi+\left[i\frac{\alpha C_{3}}{Re}-\frac{1}{\gamma
M^{2}}\frac{d\bar{\varrho}}{dy}\right]\Theta=0$ (17)
$\displaystyle\left[\frac{\bar{\kappa}}{RePr}\right]\Theta^{\prime\prime}+\left[2\alpha\frac{\bar{\mu}}{Re}\left(\gamma-1\right)\frac{M^{2}}{\Delta^{2}}\frac{d\bar{U}}{dy}\right]\psi^{\prime}+\left[-\frac{\gamma-1}{\gamma}F_{1}Q\right]\varphi^{\prime}$
$\displaystyle+\left[2\frac{C_{10}}{RePr}\right]\Theta^{\prime}+\left[2\beta\frac{\bar{\mu}}{Re}\left(\gamma-1\right)\frac{M^{2}}{\Delta^{2}}\frac{d\bar{U}}{dy}\right]\Omega^{\prime}+\left[-i\frac{\gamma-1}{\gamma}F_{1}Q\right]\psi$
$\displaystyle+\left[2i\alpha\frac{\bar{\mu}}{Re}\left(\gamma-1\right)M^{2}\frac{d\bar{U}}{dy}-\frac{\gamma-1}{\gamma}F_{1}\frac{Q}{\bar{\varrho}}\frac{d\bar{\varrho}}{dy}-\bar{\varrho}c_{p}\bar{T}\right]\varphi$
$\displaystyle+\left[-\frac{\bar{\kappa\Delta^{2}}}{RePr}+\frac{C_{8}}{Re}\left(\gamma-1\right)M^{2}+\frac{C_{9}}{RePr}-i\bar{\varrho}Q\left(c_{p}-\frac{\gamma-1}{\gamma}\right)\right]\Theta=0$
(18)
$\displaystyle\left[\frac{\bar{\mu}}{Re}\right]\Omega^{\prime\prime}+\left[\frac{\beta
C_{3}}{Re}\right]\Theta^{\prime}+\left[\frac{C_{1}}{Re}\right]\Omega^{\prime}+\left[-\beta\bar{\varrho}\frac{d\bar{U}}{dy}\right]\varphi+\left[\frac{\beta
C_{2}}{Re}\right]\Theta+\left[-\frac{\bar{\mu}}{Re}\Delta^{2}-i\bar{\varrho}Q\right]\Omega=0$
(19)

where $F_{1}=\bar{\varrho}\bar{T}/Q$,
$F_{2}=\frac{1}{\bar{\varrho}}\frac{d}{dy}\left(\frac{\bar{\varrho}^{2}\bar{T}}{Q}\right)$,
$F_{3}=\frac{d}{dy}\left(\frac{\bar{T}}{Q}\frac{d\bar{\varrho}}{dy}\right)$,
$G_{2}=F_{2}Q^{2}$, $G_{3}=\frac{dF_{1}}{dy}Q^{2}$, and $G_{4}=F_{3}Q^{2}$.

For $3D$ disturbances ($\beta\neq 0$ and $\hat{w}\neq 0$), the above set of
ODEs constitute an $8^{th}$-order system. However, when the disturbances are
$2D$ in nature ($\beta=0$ and $\hat{w}=0$), it reduces to $6^{th}$-order
system of ODEs [3]. For $2D$ disturbances, there will not be any equation for
$\Omega$ as it is trivially zero. This feature of Eqs. (16-19), which governs
the wall-normal variation of the disturbance amplitudes in a parallel
compressible boundary layer is different from its incompressible counterpart,
i.e., the Orr-Sommerfeld equation (OSE) [15]. The OSE is a $4^{th}$-order ODE
irrespective of whether the disturbances are $2D$ or $3D$.

For the stability calculations of the boundary layer on insulated wall, we
need to satisfy perturbation no-slip, zero-normal velocity and zero-normal
disturbance heat-flux conditions at the wall, i.e.,
$\hat{u}=\hat{v}=\hat{w}=\frac{\partial\hat{\theta}}{\partial y}=0$ at $y=0$.
Hence, in terms of spectral amplitudes, we need to specify the following
conditions at $y=0$,

$\displaystyle\psi(0)=\varphi(0)=\Omega(0)=0\;\;\text{and}\;\;\Theta^{\prime}(0)=0$
(20)

These homogeneous wall conditions also need to be supplemented by decaying
conditions on all perturbation variables in the free-stream, i.e.,
$\left(\psi(y),\chi(y),\varphi(y),\Upsilon(y),\Pi(y),\varrho(y),\Theta(y)\right)^{T}\rightarrow
0$ as $y\rightarrow\infty$. Next, we discuss the exponential decay rate of
various modes in the free-stream.

###  2.2 Far-field variation of disturbances

At the far-field ($y\rightarrow\infty$) all the disturbances should decay,
i.e., $\left(\hat{\textbf{v}},\hat{\rho},\hat{\theta}\right)\rightarrow 0$
whereas $(\bar{U},\bar{\varrho},\bar{T},\bar{\mu},\bar{k})\rightarrow 1$. The
Eqs. (16-19), therefore, degenerate into a system of ODEs with constant
coefficients. The exponential behavior of the modes is obtained from the
corresponding characteristic polynomial equation. We can reduce the algebraic
complexity by substituting $\varphi$, the spectral amplitude of the
perturbation $v$-velocity, with $\Pi$, the spectral amplitude of disturbance
pressure. Let, $\psi_{\infty}$, $\varphi_{\infty}$, $\Pi_{\infty}$,
$\Omega_{\infty}$, and $\Theta_{\infty}$ denote the values of $\psi$,
$\varphi$, $\Pi$, $\Theta$, and $\Omega$ at the free-stream (i.e., at
$y\rightarrow\infty$). The degenerated system of ODEs at the free-stream in
terms of $\psi_{\infty}$, $\Pi_{\infty}$, $\Theta_{\infty}$, and
$\Omega_{\infty}$ are given as

$\displaystyle\psi^{\prime\prime}_{\infty}$ $\displaystyle=$
$\displaystyle\left[\Delta^{2}+iReQ_{\infty}\right]\psi_{\infty}+b_{12}\Pi_{\infty}+b_{13}\Theta_{\infty}$
(21) $\displaystyle\Pi^{\prime\prime}_{\infty}$ $\displaystyle=$
$\displaystyle b_{22}\Pi_{\infty}+b_{23}\Theta_{\infty}$ (22)
$\displaystyle\Theta^{\prime\prime}_{\infty}$ $\displaystyle=$ $\displaystyle
b_{32}\Pi_{\infty}+b_{33}\Theta_{\infty}$ (23)
$\displaystyle\Omega^{\prime\prime}_{\infty}$ $\displaystyle=$
$\displaystyle\left[\Delta^{2}+iReQ_{\infty}\right]\Omega_{\infty}$ (24)

where, $b_{12}=\left[i\frac{Re}{\gamma
M^{2}}-Q_{\infty}\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\right]\Delta^{2}$,
$b_{13}=\left[Q_{\infty}\Delta^{2}\left(1+\frac{\bar{\lambda}}{\bar{\mu}}\right)\right]$,
$b_{22}=\tilde{b}_{22}/b_{2}$, $b_{23}=\tilde{b}_{23}/b_{2}$,
$b_{2}=\left[iQ_{\infty}\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)+\frac{Re}{\gamma
M^{2}}\right]$,
$\tilde{b}_{22}=\left[\Delta^{2}b_{2}-ReQ^{2}_{\infty}\left(1-\frac{\gamma-1}{\gamma}Pr\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)\right)\right]$,
$\tilde{b}_{23}=ReQ^{2}_{\infty}\left[1-Pr\left(2+\frac{\bar{\lambda}}{\bar{\mu}}\right)\right]$,
$b_{32}=-\left[i\frac{\gamma-1}{\gamma}RePrQ_{\infty}\right]$,
$b_{33}=\left[\Delta^{2}+iRePrQ_{\infty}\right]$, and
$Q_{\infty}=\left(\alpha-\omega\right)$. Equations (21-24) show that at the
free-stream $\Omega_{\infty}$ is decoupled from $\varphi_{\infty}$,
$\Pi_{\infty}$ and $\Theta_{\infty}$. Considering Eqs. (21-24), we note that
modes vary in the free-stream as $\phi_{j}e^{-\left(\Lambda_{j}y\right)}$
where $j=1,\cdots,8$ and the exponents $\Lambda_{j}$ are given as

$\displaystyle\Lambda_{1,2}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\left(\Delta^{2}+iReQ_{\infty}\right)}$ (25)
$\displaystyle\Lambda_{3,4}$ $\displaystyle=$
$\displaystyle\pm\frac{1}{\sqrt{2}}\sqrt{\mathcal{T}_{\infty}+\sqrt{\mathcal{T}^{2}_{\infty}-4\mathcal{D}_{\infty}}}$
(26) $\displaystyle\Lambda_{5,6}$ $\displaystyle=$
$\displaystyle\pm\frac{1}{\sqrt{2}}\sqrt{\mathcal{T}_{\infty}-\sqrt{\mathcal{T}^{2}_{\infty}-4\mathcal{D}_{\infty}}}$
(27) $\displaystyle\Lambda_{7,8}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\left(\Delta^{2}+iReQ_{\infty}\right)}$ (28)

where $\mathcal{T}_{\infty}=\left(b_{22}+b_{33}\right)$ and
$\mathcal{D}_{\infty}=\left(b_{22}b_{33}-b_{23}b_{32}\right)$. Similar
expressions of the behavior of the modes in the free-stream for compressible
boundary layers are also given in Mack [3] and Özgen & Kırcalı [34],
respectively. The free-stream variation of modes $j=1,2$ and $j=7,8$ are
identical to that corresponding to the viscous mode obtained from the Orr-
Sommerfeld equation[38, 15]. Mack[3] postulated that as $Re\rightarrow\infty$,
modes-$3$ and $4$ become independent of $M$ and depends only on $Re$ and $Pr$,
while modes-$5$ and $6$ do not depend on $Re$ and $Pr$. Modes-$3$ and $4$,
therefore, have been termed as the viscous-temperature mode, while modes-$5$
and $6$ as the inviscid mode in Mack[3]. This aspect is illustrated in the
next section. As the physically realizable disturbances decay in the free-
stream, we retain only the modes that decay in the free-stream. Hence, the
physical eigenfunctions are linear combinations of modes $j=1,3,5$, and $7$
for $3D$ disturbances and $j=1,3$, and $5$ for $2D$ disturbances. In general,
for $3D$ disturbances, therefore, we can denote
$\Phi=\left(\psi,\varphi,\Theta,\Omega\right)^{T}$ as

$\displaystyle\Phi=c_{1}\Phi_{1}+c_{3}\Phi_{3}+c_{5}\Phi_{5}+c_{7}\Phi_{7}$
(29)

whereas for $2D$ disturbances $\Phi$ can be expressed as

$\displaystyle\Phi=c_{1}\Phi_{1}+c_{3}\Phi_{3}+c_{5}\Phi_{5}$ (30)

where $c_{1}$, $c_{3}$, $c_{5}$, and $c_{7}$ are arbitrary constants. Solving
Eq. (21-24), the spectral amplitudes for the $j^{th}$ mode at the free-stream
(where $j=1,3$ and $5$) are given as

$\displaystyle\begin{Bmatrix}\psi_{j_{\infty}}\\\ \varphi_{j_{\infty}}\\\
\Theta_{j_{\infty}}\\\
\Omega_{j_{\infty}}\end{Bmatrix}=\begin{Bmatrix}b_{12}\left(b_{33}+\Lambda_{j}\right)-b_{13}b_{32}\\\
i\left[Q_{\infty}\left(b_{32}+b_{33}+\Lambda_{j}^{2}\right)\left(b_{11}+\Lambda_{j}^{2}\right)-b_{12}\left(b_{33}+\Lambda_{j}^{2}\right)+b_{13}b_{32}\right]\\\
{b_{32}\left(b_{11}+\Lambda_{j}^{2}\right)}\\\ 0\\\
\end{Bmatrix}e^{\left(-\Lambda_{j}y\right)}$ (31)

where, $b_{11}=\left(\Delta^{2}+iReQ_{\infty}\right)$. The free-stream
variation for the Squire mode corresponding to $j=7$ (relevant only for $3D$
perturbations), is given as
$\left[\psi_{7_{\infty}},\varphi_{7_{\infty}},\Theta_{7_{\infty}},\Omega_{7_{\infty}}\right]^{T}=\left[0,0,0,1\right]^{T}\;e^{\left(-\Lambda_{7}y\right)}$.

##  3 Formulation of the Compound Matrix Method (CMM) for calculation of
eigenvalues

Figure 2: $\Lambda_{1_{r}}$, $\Lambda_{3_{r}}$, and $\Lambda_{5_{r}}$ plotted
in the $(\alpha_{r},\alpha_{i})$-plane, where $\alpha=\alpha_{r}+i\alpha_{i}$,
for indicated values of $Re$ and $M$ when $\omega_{r}=0.1$ and $\beta=0$.

From the previous discussions, we note that physically realizable spectral
amplitude of $3D$ disturbances that decay in the free-stream can be written as
a linear combination of four independent modes $\Phi_{1}$, $\Phi_{3}$,
$\Phi_{5}$, and $\Phi_{7}$, respectively. These modes decay in the free-stream
as $e^{-\Lambda_{1_{r}}y}$, $e^{-\Lambda_{3_{r}}y}$, $e^{-\Lambda_{5_{r}}y}$,
and $e^{-\Lambda_{7_{r}}y}$, respectively where the $\Lambda_{j_{r}}$ denotes
the real part of $\Lambda_{j}$. In Fig. 2, we plot contours of
$\Lambda_{1_{r}}$, $\Lambda_{3_{r}}$, and $\Lambda_{5_{r}}$ in the
$(\alpha_{r},\alpha_{i})$-plane, ($\alpha=\alpha_{r}+i\alpha_{i}$) for
indicated values of $Re$ and $M$ corresponding to $\omega_{r}=0.1$,
$\omega_{i}=0$, and $\beta=0$, where $\omega=\omega_{r}+i\omega_{i}$.

Equation (25) suggests that $\Lambda_{1_{r}}$ is independent of the Mach
number, and the corresponding contours are shown in Fig. 2(a). We show the
contours of $\Lambda_{3_{r}}$ for $M=0.6$ and $4$ for $Re=1000$ and
$\omega_{r}=0.1$ in Fig. 2(b). We note that contours of $\Lambda_{3_{r}}$ are
virtually indistinguishable for $M=0.6$ and $4$. So, following Mack[3], we
identify this mode as the viscous-temperature or entropic mode. Theoretically,
for an ideal pure entropic mode, these exponents may be given as
$\Lambda^{th}_{3,4}\approx\pm\sqrt{\left(\Delta^{2}+iRePrQ_{\infty}\right)}$.
In Fig. 2(b), we also show contours of $\Lambda^{th}_{3_{r}}$ for $M=0.6$ and
$4$. We note that contour lines of actual $\Lambda_{3_{r}}$ and
$\Lambda^{th}_{3_{r}}$ are almost identical, except at higher wavenumber
regions, where slight deviation is noted. This observation justifies the
postulation of Mack[3] in terming these modes as viscous-temperature or
entropic modes.

Similarly, in Figs. 2(c,d), we show the contours of $\Lambda_{5_{r}}$ for two
different values of $Re$ for $M=0.6$ and $4$, respectively. Figure 2(c) shows
that for $M=0.6$, $\Lambda_{5_{r}}$ is practically independent of $\alpha_{i}$
when $\alpha_{r}>0.05$ as the contour lines are almost vertical. In contrast,
for $M=4$, contours of $\Lambda_{5_{r}}$ is almost horizontal for
$\alpha_{r}>0.15$ and therefore, $\Lambda_{5_{r}}$ is independent of
$\alpha_{r}$ in this range. We also observe that as the Mach number is
increased from $0.6$ to $4$, the values of $\Lambda_{5_{r}}$ reduce by one
order of magnitude. For $M=4$, slight difference is noted between the contours
for $\Lambda_{5_{r}}$ corresponding to $Re=1000$ and $3000$ for
$\alpha_{r}>0.25$. Except this small difference, both Figs. 2(c,d) show that
contours of $\Lambda_{5_{r}}$ are virtually independent of $Re$. We denote
$\Lambda^{th}_{5,6}=\pm\sqrt{\left(\Delta^{2}-Q_{\infty}^{2}M^{2}\right)}$,
which is the exponent for the free-stream variation of the modes from inviscid
instability analysis of compressible shear-layers[3]. We also show the contour
lines of $\Lambda^{th}_{5_{r}}$ in Figs. 2(c,d). We note that, while for
$M=0.6$, contour lines of actual $\Lambda_{5_{r}}$ and $\Lambda^{th}_{5_{r}}$
are identical for both $Re$ cases, for $M=4$, slight deviation is noted at
higher $\alpha_{r}$ regions. Almost perfect match of $\Lambda_{5_{r}}$ and
$\Lambda^{th}_{5_{r}}$ in Figs. 2(c,d) justifies the postulation of Mack[3] as
these to be independent of $Re$ and $Pr$ when $Re\rightarrow\infty$.
Therefore, we term this mode as the inviscid acoustic mode.

As the modes decay at different rates, the straightforward integration of Eqs.
(16-19) from free-stream to wall is not possible. It would lead to the
generation of spurious modes, which exponentially grow along the direction of
integration, causing the numerical solutions to blow up. Ng and Reid [35] and
Allen and Bridges [37] proposed the compound matrix method (CMM) to circumvent
this stiffness problem. In CMM, one solves a set of auxiliary equations
derived from the original Eqs. (16-19) in terms of compound variables. The
compound variables are well-defined combinations of the fundamental modes
$\Phi_{j}$ as defined in Eq. (29) such that these auxiliary variables grow or
decay exponentially at comparable rates. The application of the CMM removes
the stiffness of the original perturbation equations, as also illustrated in
Barker et. al. [46] and Yiantsios & Higgins [36].

Following this method, the $6^{th}$ order system for $2D$ disturbance would
yield ${}^{6}C_{3}=20$ auxiliary compound matrix equations, while this number
would be ${}^{8}C_{4}=70$ for $3D$ disturbances which follow the $8^{th}$
order stiff ODE. Using CMM, the original boundary value problem (BVP) is
converted into an initial value problem (IVP), where the initial conditions
are defined in the free-stream. The eigenvalues representing complex
streamwise wavenumber $\alpha=\alpha_{r}+i\alpha_{i}$, for a particular
combination of $Re$, $M$, $Pr$, $\omega$ and $\beta$ for spatial stability
analysis, can be found by integrating these auxiliary equations from free
stream to wall, subject to specified initial free-stream conditions, and
satisfying the dispersion relation obtained from boundary condition at the
wall as illustrated next.

The stability Eqs. (16-19) can be recast as a system of first-order ODEs as

$\displaystyle\left\\{\textbf{X}^{\prime}\right\\}=\left[E\right]\left\\{\textbf{X}\right\\}$
(32)

where the elements of the state vector
$\left\\{\textbf{X}\right\\}=\left[X_{1},X_{2},X_{3},X_{4},X_{5},X_{6},X_{7},X_{8}\right]^{T}$
are defined as $X_{1}=\psi$, $X_{2}=\psi^{\prime}$, $X_{3}=\varphi$,
$X_{4}=\varphi^{\prime}$, $X_{5}=\Theta$, $X_{6}=\Theta^{\prime}$,
$X_{7}=\Omega$ and $X_{8}=\Omega^{\prime}$. The matrix $\left[E\right]$ for
$3D$ disturbances is a $8\times 8$ matrix. Similarly for the $2D$ case,
$\left[E\right]$ is a $6\times 6$ matrix. Following Eq. (29), we can express X
as as
$\textbf{X}=c_{1}\textbf{X}_{1}+c_{3}\textbf{X}_{3}+c_{5}\textbf{X}_{5}+c_{7}\textbf{X}_{7}$,
where $\textbf{X}_{1}$, $\textbf{X}_{3}$, $\textbf{X}_{5}$ and
$\textbf{X}_{7}$ are the linearly independent physical modes that decay in the
free-stream.

Following the methodology proposed in Allen and Bridges [37] for a fourth-
order-system, we project the solutions of the above eighth-order system on a
subspace of $C^{8}$ into $\bigwedge^{4}(C^{8})$ with the help of
$\textbf{X}_{1}$, $\textbf{X}_{3}$, $\textbf{X}_{5}$, and $\textbf{X}_{7}$.
The problem is, therefore, reduced to linking these four-dimensional subspaces
of $C^{8}$ satisfying Eqs. (16-19) with a corresponding point in the vector
space $\bigwedge^{4}(C^{8})$. Mathematically, any subspace spanned by four
linearly independent vectors $e_{1}$, $e_{2}$, $e_{3}$, and $e_{4}$ can be
represented notationally as a point $e_{1}\wedge e_{2}\wedge e_{3}\wedge
e_{4}$, in $\bigwedge^{4}(C^{8})$. The boundary conditions at $y=0$ define a
specific four-dimensional subspace of $C^{8}$, as illustrated later.

Let us define the new basis variables as
$e_{j}=[X_{j,1},X_{j,3},X_{j,5},X_{j,7}]^{T}$, where $j=1,\cdots,8$ in
$C^{8}$. Therefore, all elements of $e_{j}\wedge e_{k}\wedge e_{l}\wedge
e_{m}$ form the basis for $\bigwedge^{4}(C^{8})$ with the dimension
${}^{8}C_{4}=70$. The complete solution matrix, with these basis vectors, is
obtained as

$\displaystyle\begin{bmatrix}X_{1,1}&X_{1,3}&X_{1,5}&X_{1,7}\\\
X_{2,1}&X_{2,3}&X_{2,5}&X_{2,7}\\\ X_{3,1}&X_{3,3}&X_{3,5}&X_{3,7}\\\
X_{4,1}&X_{4,3}&X_{4,5}&X_{4,7}\\\ X_{5,1}&X_{5,3}&X_{5,5}&X_{5,7}\\\
X_{6,1}&X_{6,3}&X_{6,5}&X_{6,7}\\\ X_{7,1}&X_{7,3}&X_{7,5}&X_{7,7}\\\
X_{8,1}&X_{8,3}&X_{8,5}&X_{8,7}\end{bmatrix}=\begin{bmatrix}\psi_{1}&\psi_{3}&\psi_{5}&\psi_{7}\\\
\psi^{\prime}_{1}&\psi^{\prime}_{3}&\psi^{\prime}_{5}&\psi^{\prime}_{7}\\\
\varphi_{1}&\varphi_{3}&\varphi_{5}&\varphi_{7}\\\
\varphi^{\prime}_{1}&\varphi^{\prime}_{3}&\varphi^{\prime}_{5}&\varphi^{\prime}_{7}\\\
\Theta_{1}&\Theta_{3}&\Theta_{5}&\Theta_{7}\\\
\Theta^{\prime}_{1}&\Theta^{\prime}_{3}&\Theta^{\prime}_{5}&\Theta^{\prime}_{7}\\\
\Omega_{1}&\Omega_{3}&\Omega_{5}&\Omega_{7}\\\
\Omega^{\prime}_{1}&\Omega^{\prime}_{3}&\Omega^{\prime}_{5}&\Omega^{\prime}_{7}\end{bmatrix}$
(33)

Here, $X_{j,k}$ represents the $j^{th}$ element of the solution vector X
corresponding to $k^{th}$ mode, where $j=1,\cdots,8$ and $k=1,3,5,7$.
Therefore, the seventy-two compound variables are constructed as the $4\times
4$ minors of the solution matrix given by Eq. (34). This can be symbolically
denoted as

$\displaystyle
Y_{n}=\mathcal{Y}_{j,k,l,m}=\begin{vmatrix}X_{j,1}&X_{j,3}&X_{j,5}&X_{j,7}\\\
X_{k,1}&X_{k,3}&X_{k,5}&X_{k,7}\\\ X_{l,1}&X_{l,3}&X_{l,5}&X_{l,7}\\\
X_{m,1}&X_{m,3}&X_{m,5}&X_{m,7}\end{vmatrix}$ (34)

where $1\leq j<k<l<m\leq 8$ and $n=n(j,k,l,m)$. The details of the functional
relationship between the indices are given in Appendix-$I$. Denoting
$\left\\{\textbf{Z}\right\\}=\left[Y_{1},\cdots,Y_{70}\right]^{T}$, we get
that $Z$ satisfies the set of linear coupled ODEs given by

$\displaystyle\left\\{\textbf{Z}^{\prime}\right\\}=\left[F\right]\left\\{\textbf{Z}\right\\}$
(35)

Here, $[F]$ is $70\times 70$ matrix for $3D$ disturbance. We can determine the
elements of the $[F]$ matrix from the elements of the matrix $[E]$ given in
Eq. (32). The detailed methodology is illustrated in Appendix-$II$. One
readily notes from Eq. (34) that at the free-stream
$Y_{n_{\infty}}\simeq\mathcal{C}_{n}\exp\left(-\left(\Lambda_{1}+\Lambda_{3}+\Lambda_{5}+\Lambda_{7}\right)y\right)$,
where $\mathcal{C}_{n}$ is some constant. Therefore, all the modes in the
free-stream decay at the identical exponential rate
$e^{\left(-\left(\Lambda_{1}+\Lambda_{3}+\Lambda_{5}+\Lambda_{7}\right)y\right)}$,
removing the stiffness of the original problem. Thus, we can solve Eq. (35) by
using any standard integration procedure and no special treatment like
repeated orthonormalization of the solution as adopted in Mack[33] and Özgen &
Kırcalı [34] is required.

At the free-stream ($y\rightarrow\infty$) exact analytic nature of the modes
$\Phi_{1}$, $\Phi_{3}$, $\Phi_{5}$ and $\Phi_{7}$ are known from Eq. (31).
This equation makes it possible to specify initial conditions for the compound
variables $Y_{n}$ at the free-stream. Using the definitions of $Y_{n}$ given
in Eq. (34) where $n=1,\cdots,70$, and the analytical behavior of $\Phi_{j}$
(here, $j=1,3,5,7$) at the free-stream given by Eq. (31), one obtains the
initial conditions for $Y_{n}$ at $y\rightarrow\infty$ as

$\displaystyle
Y_{n_{\infty}}=\mathcal{Y}_{{j,k,l,m}_{\infty}}=\begin{vmatrix}X_{{j,1}_{\infty}}&X_{{j,3}_{\infty}}&X_{{j,5}_{\infty}}&X_{{j,7}_{\infty}}\\\
X_{{k,1}_{\infty}}&X_{{k,3}_{\infty}}&X_{{k,5}_{\infty}}&X_{{k,7}_{\infty}}\\\
X_{{l,1}_{\infty}}&X_{{l,3}_{\infty}}&X_{{l,5}_{\infty}}&X_{{l,7}_{\infty}}\\\
X_{{m,1}_{\infty}}&X_{{m,3}_{\infty}}&X_{{m,5}_{\infty}}&X_{{m,7}_{\infty}}\end{vmatrix}e^{\left(-\left(\Lambda_{1}+\Lambda_{3}+\Lambda_{5}+\Lambda_{7}\right)y_{\infty}\right)}$
(36)

The system of equation given by Eq. (35) now can be integrated as an initial
value problem (IVP) from $y=y_{\infty}$ to wall i.e., $y=0$ by any standard
ODE solving technique like forth-order Runge-Kutta method $RK_{4}$. While
using $Y_{n_{\infty}}$ as the initial condition for Eq. (35) at
$y=y_{\infty}$, it is preferred to scale the corresponding values by
$e^{\left(-\left(\Lambda_{1}+\Lambda_{3}+\Lambda_{5}+\Lambda_{7}\right)y_{\infty}\right)}$.

To find out the eigenvalues of Eqs. (16-19), one needs to satisfy the
homogeneous boundary conditions at $y=0$ given by Eq. (20). From Eq. (20), we
can write at $y=0$

$\displaystyle\psi(0)$ $\displaystyle=$ $\displaystyle
c_{1}\psi_{1}(0)+c_{3}\psi_{3}(0)+c_{5}\psi_{5}(0)+c_{7}\psi_{7}(0)$ (37)
$\displaystyle=$ $\displaystyle
c_{1}X_{1,1}(0)+c_{3}X_{1,3}(0)+c_{5}X_{1,5}(0)+c_{7}X_{1,7}(0)=0$
$\displaystyle\varphi(0)$ $\displaystyle=$ $\displaystyle
c_{1}\varphi_{1}(0)+c_{3}\varphi_{3}(0)+c_{5}\varphi_{5}(0)+c_{7}\varphi_{7}(0)$
(38) $\displaystyle=$ $\displaystyle
c_{1}X_{3,1}(0)+c_{3}X_{3,3}(0)+c_{5}X_{3,5}(0)+c_{7}X_{3,7}(0)=0$
$\displaystyle\Theta^{\prime}(0)$ $\displaystyle=$ $\displaystyle
c_{1}\Theta^{\prime}_{1}(0)+c_{3}\Theta^{\prime}_{3}(0)+c_{5}\Theta^{\prime}_{5}(0)+c_{7}\Theta^{\prime}_{7}(0)$
(39) $\displaystyle=$ $\displaystyle
c_{1}X_{6,1}(0)+c_{3}X_{6,3}(0)+c_{5}X_{6,5}(0)+c_{7}X_{6,7}(0)=0$
$\displaystyle\Omega(0)$ $\displaystyle=$ $\displaystyle
c_{1}\Omega_{1}(0)+c_{3}\Omega_{3}(0)+c_{5}\Omega_{5}(0)+c_{7}\Omega_{7}(0)$
(40) $\displaystyle=$ $\displaystyle
c_{1}X_{7,1}(0)+c_{3}X_{7,3}(0)+c_{5}X_{7,5}(0)+c_{7}X_{7,7}(0)=0$

To have a non-trivial solution, therefore, we need to satisfy the
characteristic determinant of the linear system of equations to be identically
zero at the wall, i.e.,

$\displaystyle
Y_{23}(0)=\mathcal{Y}_{1,3,6,7}(0)=\begin{vmatrix}\psi_{1}(0)&\psi_{3}(0)&\psi_{5}(0)&\psi_{7}(0)\\\
\varphi_{1}(0)&\varphi_{3}(0)&\varphi_{5}(0)&\varphi_{7}(0)\\\
\Theta^{\prime}_{1}(0)&\Theta^{\prime}_{3}(0)&\Theta^{\prime}_{5}(0)&\Theta^{\prime}_{7}(0)\\\
\Omega_{1}(0)&\Omega_{3}(0)&\Omega_{5}(0)&\Omega_{7}(0)\end{vmatrix}=\begin{vmatrix}X_{1,1}(0)&X_{1,3}(0)&X_{1,5}(0)&X_{1,7}(0)\\\
X_{3,1}(0)&X_{3,3}(0)&X_{3,5}(0)&X_{3,7}(0)\\\
X_{6,1}(0)&X_{6,3}(0)&X_{6,5}(0)&X_{6,7}(0)\\\
X_{7,1}(0)&X_{7,3}(0)&X_{7,5}(0)&X_{7,7}(0)\end{vmatrix}=0\hskip 28.45274pt$
(41)

The above equation provides the dispersion relation as
$D_{r}+iD_{i}=Y_{23}(0)=\mathcal{Y}_{1,3,6,7}(0)=0$. Integrating the resultant
auxiliary system of Eqs. (35) subject to the initial conditions given
corresponding to the normalized $Y_{n_{\infty}}$ at $y=y_{max}$, and
satisfying the above dispersion relation, we can detect the eigenvalues for a
given set of parameters.

For $2D$ disturbances, Eqs. (16-19) constitute a $6^{th}$-order system, and
hence, there would be three physical linear independent modes $\Phi_{1}$,
$\Phi_{3}$, and $\Phi_{5}$ that decay in the free-stream. Consequently, the
application of CMM would yield ${}^{6}C_{3}=20$ compound matrix equations. The
rest of the procedure is similar to what is described above.

##  4 Variation of mean wall-bounded shear layer over adiabatic flat plate

Figure 3: (a) $\bar{U}(y)$, and (b) $\bar{T}(y)$, plotted as a function of
$\hat{\eta}=\tilde{y}\sqrt{Re_{x}}/\tilde{x}$, where
$Re_{x}=\tilde{\rho}_{\infty}\tilde{U}_{\infty}\tilde{x}/\tilde{\mu}_{\infty}$
is the Reynolds number based on the streamwise coordinate $\tilde{x}$. (c) The
displacement thickness parameter $c_{\delta}$ plotted as a function of the
free-stream Mach number $M$ where,
$c_{\delta}=\delta_{*}\sqrt{Re_{x}}/\tilde{x}$.

For a $2D$ compressible wall-bounded shear layer, the self-similar equations
can be obtained from the corresponding boundary-layer equations[47] following
Illingworth [48] or Howarth–Dorodnitsyn [43] transformation from the physical
$(\tilde{x},\tilde{y})$-plane to transformed $(\xi,\eta)$-coordinate system,
where
$\xi=\int^{\tilde{x}}_{0}\tilde{\rho}_{e}\tilde{\mu}_{e}\tilde{U}_{e}d\tilde{x}$,
and
$\eta=\frac{\tilde{U}_{e}}{\sqrt{\xi}}\int^{\tilde{y}}_{0}\tilde{\rho}d\tilde{y}$.
Here, $\tilde{\rho}_{e}$, $\tilde{U}_{e}$ and $\tilde{\mu}_{e}$ indicate
density, streamwise velocity, and dynamic viscosity at the edge of the shear-
layer. Here, we are not considering any external pressure gradient to be
applied, and hence the edge conditions are not functions of streamwise
coordinate and are treated to be constants, i.e.,
$\tilde{\rho}_{e}=\tilde{\rho}_{\infty}$, $\tilde{U}_{e}=\tilde{U}_{\infty}$,
and $\tilde{\mu}_{e}=\tilde{\mu}_{\infty}$. The self-similar equations are
given as [43]

$\displaystyle\left(\bar{\varrho}\bar{\mu}f_{\eta\eta}\right)_{\eta}+ff_{\eta\eta}$
$\displaystyle=$ $\displaystyle 0$ (42)
$\displaystyle\left(\frac{1}{Pr}\bar{\varrho}\bar{\kappa}g_{\eta}\right)_{\eta}+\bar{c}_{p}fg_{\eta}+\left(\gamma-1\right)M^{2}\bar{\varrho}\bar{\mu}f_{\eta\eta}^{2}$
$\displaystyle=$ $\displaystyle 0$ (43)

where $\left(\cdot\right)_{\eta}=d\left(\cdot\right)/d\eta$. The self-similar
variables $f(\eta)$ and $g(\eta)$ are defined such that $f_{\eta}=\bar{U}$ and
$g=\bar{T}$. For no-slip adiabatic wall, one needs to satisfy
$f(0)=f_{\eta}(0)=g_{\eta}(0)=0$ while at the free-stream
$\left(f_{\eta},g\right)\rightarrow 1$ as $\eta\rightarrow\infty$. Following
the boundary layer approximation,
$\frac{\partial\tilde{p}}{\partial\tilde{y}}=0$ in the shear layer, and hence,
considering the ideal equation of state is obeyed, $\bar{\varrho}\bar{T}=1$.
Here, we consider the dynamic viscosity $\tilde{\mu}$, heat conductivity
$\tilde{\kappa}$ and specific heat at constant pressure $\tilde{c}_{p}$ are
functions of temperature only[49] as

$\displaystyle\bar{\mu}$ $\displaystyle=$
$\displaystyle\bar{T}^{3/2}\left(\frac{1+110/T_{ref}}{\bar{T}+110/T_{ref}}\right)$
(44) $\displaystyle\bar{\kappa}$ $\displaystyle=$
$\displaystyle\bar{T}^{1/2}\left(\frac{1+S_{3}10^{-S_{4}}}{1+\left(S_{3}/\bar{T}\right)10^{-S_{4}/\bar{T}}}\right)$
(45) $\displaystyle\bar{c}_{p}$ $\displaystyle=$
$\displaystyle\left(\frac{1+\frac{\gamma-1}{\gamma}{\left(\theta_{2}^{2}e^{-\theta_{2}}\right)}/{\left(1-e^{\theta_{2}}\right)^{2}}}{1+\frac{\gamma-1}{\gamma}{\left(\theta_{1}^{2}e^{-\theta_{1}}\right)}/{\left(1-e^{\theta_{1}}\right)^{2}}}\right)$
(46)

where $S_{3}=\left(245.4/T_{ref}\right)$, $S_{4}=12/T_{ref}$,
$\theta_{1}=3055/T_{ref}$, $\theta_{2}=\theta_{1}/\bar{T}$ and $T_{ref}=303K$.
The variation of dynamics viscosity denoted in Eq. (44) is the well-known
Sutherland’s law [49, 48] for air. In Fig. 3(a,b), we plot the self-similar
velocity $\bar{U}$, and temperature $\bar{T}$ as a function of
$\hat{\eta}=\tilde{y}\sqrt{Re_{x}}/\tilde{x}$ for the adiabatic plate
condition. Here,
$Re_{x}=\tilde{\rho}_{\infty}\tilde{U}_{\infty}\tilde{x}/\tilde{\mu}_{\infty}$
is the Reynolds number based on local streamwise coordinate $\tilde{x}$. We
note that the increase in Mach number increases the hydrodynamic and thermal
boundary layer thickness. As noted in Fig. 3(b), this also increases the wall
temperature. We use the displacement thickness $\delta_{*}$ as the reference
length scale in the stability calculations to determine the Reynolds number
$Re$. For any self-similar wall-bounded laminar shear layer,
$\delta_{*}=c_{\delta}\tilde{x}/\sqrt{Re_{x}}$, where the factor $c_{\delta}$
is a function of Mach number $M$, $\gamma$ and Prandtl number $Pr$. For an
incompressible, isothermal, and laminar flow, $c_{\delta}=1.72$ [48]. In Fig.
3(c), we show the variation of $c_{\delta}$ as a function of $M$. As the Mach
number $M$ increases, $c_{\delta}$ monotonically increases.

##  5 Results and discussion

Traditionally, two approaches have been adopted to perform stability analysis
of boundary layers, namely (1) temporal stability analysis, which considers
complex $\omega$ with real $\alpha$ and $\beta$, and (2) spatial stability
theory, which assumes complex $\alpha$ (and $\beta$) treating $\omega$ as
real. A more elaborate and realistic analysis can be performed using
spatiotemporal stability analysis where $\alpha$, $\beta$, and $\omega$ are
all treated as complex. Detailed description in this regard may be found in
Drazin and Reid [15], Tollmien [14], Schmid et. al. [45], and Sengupta [38].
Here, we only focus on the spatial stability analysis and treat $\omega$ and
$\beta$ as real and $\alpha=\alpha_{r}+i\alpha_{i}$ as complex. Therefore,
$\omega_{i}=0$, where $\omega=\omega_{r}+i\omega_{i}$. As illustrated before,
spatial instability occurs when $\alpha_{i}<0$.

###  5.1 Comparison of spatial stability between incompressible and
compressible boundary layers at $M=0.1$

Figure 4: Comparison between the spatial stability analysis of 2D
incompressible zero-pressure gradient Blasius and compressible boundary layer
on the insulated wall at $M=0.1$. (a,b) Spatial eigenvalues for $Re=1500$ and
$\omega_{r}=0.1$, (c) comparison of the corresponding Neutral curves, and
(d,e) group velocity corresponding to mode-$1$ in the $(Re,\omega_{r})$-plane.
(a) CMM applied to OSE for 2D incompressible ZPG boundary layer
---
Mode No. | $\alpha_{r}$ | $\alpha_{i}$ | $Y_{1r}(0)$ | $Y_{1i}(0)$
$1$ | $0.29373724$ | $-7.03994013\times 10^{-3}$ | $1.6941\times 10^{-19}$ | $6.8821\times 10^{-20}$
$2$ | $0.17675906$ | $0.12104521$ | $-2.9461\times 10^{-16}$ | $2.0643\times 10^{-17}$
$3$ | $0.13301162$ | $0.21791196$ | $-9.1898\times 10^{-14}$ | $2.4733\times 10^{-14}$
(b) CMM applied to 2D compressible boundary layer over adiabatic plate for
$M=0.1$
Mode No. | $\alpha_{r}$ | $\alpha_{i}$ | $Y_{7r}(0)$ | $Y_{7i}(0)$
$1$ | $0.29324967$ | $-5.9321327\times 10^{-3}$ | $-0.1728\times 10^{-8}$ | $-0.2485\times 10^{-8}$
$2$ | $0.17673123$ | $0.12151321$ | $-0.5690\times 10^{-5}$ | $-0.1086\times 10^{-4}$
$3$ | $0.21522309$ | $0.21439414$ | $-0.2256\times 10^{-4}$ | $0.2402\times 10^{-4}$
Table 1: The spatial eigenvalues obtained after Newton-Raphson polishing
tabulated for $Re=1500$ and $\omega_{r}=0.1$ for 2D zero-pressure gradient
incompressible and compressible boundary layer at $M=0.1$ over an adiabatic
wall. The corresponding values of the real and imaginary part of the
dispersion relations, after achieving convergence, are also enlisted.

We first compare the spatial stability analysis of incompressible $2D$ zero-
pressure gradient (ZPG) Blasius boundary layer and compressible boundary layer
on the insulated wall at $M=0.1$ for the validation of the developed approach
using CMM. We choose a Reynolds number of $Re=1500$ and $\omega_{r}=0.1$. The
stability analysis of the incompressible Blasius boundary layer is carried out
by solving the Orr-Sommerfeld equation (OSE)[15]

$\displaystyle\varphi^{\prime\prime\prime\prime}-2\Delta^{2}\varphi^{\prime\prime}+\Delta^{4}\varphi=iRe\left[Q\left(\varphi^{\prime\prime}-\Delta^{2}\varphi\right)-\alpha\frac{d^{2}\bar{U}}{dy^{2}}\varphi\right]$
(47)

Equation (47) shows that OSE is a forth-order ODE. When CMM is applied to the
OSE, the corresponding dispersion relation is obtained as the no-slip and the
zero-normal boundary condition at the wall. These conditions lead to the
dispersion relation given as
$Y_{1}=\left(\varphi_{1}\varphi_{3}^{\prime}-\varphi_{1}^{\prime}\varphi_{3}\right)=0$
at the wall, where $\varphi_{1}$ and $\varphi_{3}$ are the independent modes
of the OSE, which decays in the free-stream[15, 38]. In Fig. 4(a), we plot the
contours of $Y_{1r}(0)=0$ and $Y_{1i}(0)=0$ in the
$(\alpha_{r},\alpha_{i})$-plane as obtained by integrating the CMM equation
corresponding to OSE from free-stream to the wall. The intersection of the
contour lines corresponding to $Y_{1r}(0)=0$ an $Y_{1i}(0)=0$ indicates a
spatial eigenvalue for OSE. We note the existence of three intersection
points, which are noted as modes-$1$, $2$, and $3$, respectively. The obtained
values, noted in Fig. 4(a) for the three modes, can be further polished by the
Newton-Raphson method using CMM. These values are noted in Table 1(a). The
value of $\alpha_{i}$ for mode-$1$ is noted to be negative, indicating this
mode displays spatial instability, and disturbances corresponding to this mode
grow in space as $e^{-\alpha_{i}x}$. The other two modes are spatially stable.

Figure 4(b) shows the contours lines of $Y_{7r}(0)=0$ and $Y_{7i}(0)=0$ when
CMM is applied to a compressible boundary layer on the adiabatic plate for
$M=0.1$. The dispersion relation is given as $Y_{7}(0)=0$ as discussed in Sec.
3 and therefore, the intersection points of $Y_{7r}(0)=0$ and $Y_{7i}(0)=0$
indicate eigenmodes. Figure 4(b) shows that modes-$1$, $2$, and $3$ are
clearly noted. These values are subsequently polished using the Newton-Raphson
method as tabulated in Table 1(b). The convergence is achieved when
$|\delta\alpha|=|\alpha_{n+1}-\alpha_{n}|<10^{-16}$ and $|\delta
Y_{7}(0)|=|Y_{7,n+1}(0)-Y_{7,n}(0)|<10^{-16}$, where the subscript $n$
indicates the $n^{th}$ iteration of the Newton-Raphson method. Figure 4(b)
also shows that there are significant numerical disturbances when
$\alpha_{r}<0.14$. Even if several crossings are noted on the top part of the
domain, these are spurious as no convergence in the Newton-Raphson method
corresponding to these values is achieved. We note that modes-$1$ and $2$ are
similar for incompressible and compressible boundary layers. Mode-$3$, for the
compressible boundary layer, has a higher value of $\alpha_{r}$ than the
incompressible one. Like the incompressible one, the mode-$1$ for the
compressible boundary layer displays spatial instability, and modes-$2$ and
$3$ are spatially stable.

We find that for both 2D incompressible and $M=0.1$ compressible boundary
layer on the insulated wall, modes-$2$ and $3$ are always spatially stable
when $\omega_{r}$ and $Re$ is varied over a wide range, i.e., $\alpha_{i}>0$
always for these modes. These modes are also noted to disappear when
$\omega_{r}$ and $Re$ are below certain corresponding values, as also reported
in Sengupta, Nair, & Rana [50] for ZPG incompressible boundary layer. Next, we
track mode-$1$ in the $(Re,\omega_{r})$-plane by performing Newton-Raphson
polishing while using a neighboring polished mode as the initial guess. In
Fig. 4(c), we compare $\alpha_{i}=0$ contour for mode-$1$ for both cases. This
curve is called the Neutral curve [15]. Inside this curve, $\alpha_{i}<0$, and
the flow is spatially unstable, while outside this curve, $\alpha_{i}>0$, and
the flow is spatially stable. We define the critical Reynolds number $Re_{cr}$
as the lowermost limiting Reynolds number below which the mode-$1$ is always
stable. We note identical $Re_{cr}$ ($Re_{cr}\simeq 520$) and lower branch of
the neutral curve for both 2D cases. The value of $Re_{cr}\simeq 520$ for
incompressible boundary layers is reported in various literature [15, 16, 45,
38]. The upper branch of the neutral curve shows a slight difference due to
the effect of flow compressibility. We note that the $M=0.1$ compressible
boundary layer is slightly more stable than the incompressible one. In Fig.
4(c,d), we further compare the group velocity $V_{g}$ of disturbances
corresponding to mode-$1$. The group velocity is computed as
$V_{g}=\partial\omega_{r}/\partial\alpha_{r}$ at fixed $Re$, which indicates
the speed at which the wave energy travels. Figure 4(c,d) shows that the group
velocity contours for compressible and incompressible shear layers are almost
identical.

###  5.2 Effects of flow compressibility on spatial stability

Figure 5: The Neutral curve plotted for indicated Mach number cases in (a)
$(\sqrt{Re_{x}},\omega_{r})$\- and (b) $(Re,\omega_{r})$-plane for $2D$
disturbance. Group velocity contours shown for (c) $M=0.6$ and (d) $M=2$ in
the $(Re,\omega_{r})$-plane. The corresponding neutral curve in frames (c,d)
is represented as a dashed line.

We further explore the effects of flow compressibility on the spatial
stability of the compressible wall-bounded shear layer on an adiabatic wall.
Here, we consider the Mach number cases ranging from $0.1$ to $2$ and only
their corresponding $2D$ stability characteristics. As noted in Fig. 4(a),
multiple spatial eigenmodes exist for each case. However, only one mode
displays spatial instability, and all other modes are spatially stable, i.e.,
$\alpha_{i}>0$. In Figs. 5(a,b), we show the neutral curve for the only one
spatially unstable mode plotted in (a) $(\sqrt{Re_{x}},\omega_{r})$\- and (b)
$(Re,\omega_{r})$-planes, where $Re_{x}$ is the Reynolds number based on local
streamwise coordinate and $Re=c_{\delta}\sqrt{Re_{x}}$ as described in Sec. 4.
Here, we provide two different measures as the value of $c_{\delta}$ depends
on the $M$ (see Fig. 3(c)), indicating $Re$ not to be a good measure of the
extent of the unstable zone when different Mach number cases are compared.
Figures 5(a,b) show that with an increase in Mach number, the extent of the
spatially unstable zone reduces. This reduction is marginal for subsonic cases
(especially when viewed at the $(\sqrt{Re_{x}},\omega_{r})$-plane) and
subsequently becomes more severe as the flow becomes supersonic. An
explanation in this regard was given in Mack[3], which concluded that the
viscous instability weakens with an increase in flow Mach number, while the
dominance of inviscid instability extends to lower Reynolds numbers. We show
the group velocity contours for the spatially unstable modes in Figs. 5(c,d)
corresponding to $M=0.6$ and $2$, respectively, which shows that at higher
Mach numbers, the disturbances travel relatively faster. The contour lines
become increasingly parallel while the deviation in the values becomes
significantly lesser and is closer to $0.5$, i.e., the speed of a pure
acoustic disturbance. For $M=2$, these indicate the lessening of viscous
dispersive effects, and the disturbances increasingly travel following
inviscid acoustic modes as postulated in Mack[3].

###  5.3 Effects of spanwise wavenumber $\beta$ on spatial stability

Figure 6: The Neutral curve plotted in $(Re,\omega_{r})$-plane corresponding
to $2D$ compressible boundary layer on insulated plate corresponding to
indicated spanwise wavenumber cases for (a) $M=1.8$ and (b) $M=2$. Group
velocity contours shown for $M=2$ case in $(Re,\omega_{r})$-plane when (c)
$\beta=0.1$ and (d) $\beta=0.5$. The corresponding neutral curve in frames
(c,d) is represented as a dashed line.

Next, we investigate the effects of spanwise wavenumber $\beta$ on the spatial
stability of the compressible boundary layer over an adiabatic wall. We
consider $M=1.8$ and $2$ with the spanwise wavenumber $\beta$ ranging from $0$
to $0.5$. As the problem is essentially $3D$, the CMM is applied to the
$8^{th}$-order system (even for the case corresponding to $\beta=0$), which
results in $72$ compound matrix equations.

Multiple spatial eigenmodes exist for each of these cases, while only one mode
displays spatial instability, and the rest are always spatially stable. In
Figs. 6(a,b), we show the neutral curve for the sole spatially unstable mode
in $(Re,\omega_{r})$-plane for $M=1.8$ and $2$, respectively. We observe that
the value of the critical Reynolds number $Re_{cr}$ and the extent of the
spatially unstable zone depends on $\beta$. For $M=1.8$, the neutral curve for
$\beta=0$ and $\beta=0.1$ is almost indistinguishable. The subsequent increase
in $\beta$ decreases $Re_{cr}$ and increases the extent of the unstable zone
as noted for $\beta=0.3$ in Fig. 6(a). However, a subsequent increase in
$\beta$ increases $Re_{cr}$ and decreases the extent of the spatially unstable
zone. Beyond $\beta=0.4$, the decrease in the extent of the unstable zone is
more drastic, as noted by comparing the neutral curves for $\beta=0.4$,
$0.45$, and $0.5$ in Fig. 6(a). Therefore, an optimum value of $\beta$ close
to $0.3$ exists for which $Re_{cr}$ is minimum, and the extent of the
spatially unstable zone is maximum for $M=1.8$. Similar conclusions are also
noted for $M=2$ in Fig. 6(b). For $M=2$, the decrease in $Re_{cr}$ and
increase in the extent in the unstable zone from $\beta=0$ to $0.3$ is more
severe than $M=1.8$. Subsequently, we note an increase in $Re_{cr}$ and a
decrease in the extent of the unstable zone, which is less severe than the
$M=1.8$ case, as the extent of the unstable zone for $\beta=0.5$ is still
larger than $\beta=0$. The group velocity contours for $\beta=0.1$ and $0.5$
corresponding to $M=2$ are depicted in Fig. 6(c,d), which shows that
disturbances are less dispersive for $\beta=0.1$ than $0.5$ for lower
frequency components ($\omega_{r}\lesssim 0.15$). These components,
collectively as a group, are also noted to travel slower for $\beta=0.5$ than
$\beta=0.1$.

Regarding the stability calculation of the compressible boundary layers, there
is ambiguity regarding the appropriate order of the equation. For the
incompressible flows, whether the disturbances are $2D$ or $3D$ ($\beta\neq 0$
or $\hat{w}\neq 0$), the stability equations are given by the appropriate form
of the $4^{th}$-order OSE given by Eq. (47). For the compressible boundary
layers, such simplifications are not straight-forward, as also noted in
Mack[3], Dunn & Lin [51] and Lees & Reshotko [52]. Even when $\beta=0$, the
$3D$ compressible stability equations are not identical to the $2D$ version.
Considering Eqs. (16-19) we note that $\Omega$ contributes only to the
equation for $\Theta$ through the dissipation term
$\mathcal{T}_{1}=\left[2\beta\frac{\bar{\mu}}{Re}\left(\gamma-1\right)\frac{M^{2}}{\Delta^{2}}\frac{d\bar{U}}{dy}\right]$.
Mack[3] proposed to simply set this term to zero while using the resultant
$6^{th}$-order system for eigenvalue calculations so that the accuracy is
improved specifically at higher supersonic Mach numbers. We note from Eqs.
(16-19) that for the $2D$ case, we only need to satisfy Eqs. (16-18) with
$\mathcal{T}_{1}=0$ corresponding to
$\psi(0)=\varphi(0)=\Theta^{\prime}(0)=0$, while for the $3D$ case with
$\beta=0$, we also need to satisfy an additional decoupled equation

$\displaystyle\left[\frac{\bar{\mu}}{Re}\right]\Omega^{\prime\prime}+\left[\frac{C_{1}}{Re}\right]\Omega^{\prime}+\left[-\frac{\bar{\mu}}{Re}\Delta^{2}-i\bar{\varrho}Q\right]\Omega=0$
(48)

along with $\Omega(0)=0$. Therefore, it can be deduced that an eigenvalue for
the $2D$ case will represent the eigenvalue for the $3D$ case with $\beta=0$
corresponding to the degenerative trivial condition of $\hat{w}=0$ or
$\Omega=0$. An eigenvalue for the $3D$ case with $\beta=0$ should not
necessarily represent an eigenvalue for the respective $2D$ case, as for the
former, non-trivial Squire mode is not imposed at free-stream while using CMM.
These points are illustrated next using Fig. 7 and Table 2.

Figure 7: The neutral curve plotted in $(Re,\omega_{r})$-plane corresponding to $2D$ ($6^{th}$-order system) and $3D$ disturbance case with $\beta=0$ ($8^{th}$-order system) for (a) $M=0.6$, (b) $M=1.4$, (c) $M=1.8$ and (d) $M=2$. The corresponding $Re_{cr}$ is also indicated in the frames. (a) $Re=2500$, $\omega_{r}=0.06$ |  |
---|---|---
$M$ | Formulation | $\alpha_{r}$ | $\alpha_{i}$ | $Y_{jr}(0)$ | $Y_{ji}(0)$
$0.6$ | $2D$ | $0.18893069$ | $-9.41920470\times 10^{-3}$ | $3.6520\times 10^{-9}$ | $1.5799\times 10^{-9}$
$3D$ with $\beta=0$ | $0.18952319$ | $-6.22944070\times 10^{-3}$ | $1.2104\times 10^{-9}$ | $1.4429\times 10^{-9}$
$1.4$ | $2D$ | $0.14093782$ | $-3.36921910\times 10^{-3}$ | $1.4945\times 10^{-8}$ | $-1.5668\times 10^{-8}$
$3D$ with $\beta=0$ | $0.14325045$ | $-2.94359420\times 10^{-3}$ | $1.5272\times 10^{-10}$ | $-3.8199\times 10^{-10}$
$1.8$ | $2D$ | $0.11736789$ | $-6.51037790\times 10^{-4}$ | $-2.9220\times 10^{-9}$ | $-6.0412\times 10^{-10}$
$3D$ with $\beta=0$ | $0.12099984$ | $-8.21889200\times 10^{-4}$ | $1.4344\times 10^{-9}$ | $1.8298\times 10^{-9}$
$2.0$ | $2D$ | $0.10879819$ | $-7.24156630\times 10^{-5}$ | $-1.4542\times 10^{-9}$ | $-2.3321\times 10^{-9}$
$3D$ with $\beta=0$ | $0.11232695$ | $-2.36114590\times 10^{-4}$ | $1.2846\times 10^{-9}$ | $-1.4493\times 10^{-10}$
(b) $Re=2500$, $\omega_{r}=0.1$ |  |
$M$ | Formulation | $\alpha_{r}$ | $\alpha_{i}$ | $Y_{jr}(0)$ | $Y_{ji}(0)$
$0.6$ | $2D$ | $0.28802121$ | $4.43784240\times 10^{-3}$ | $2.1877\times 10^{-9}$ | $2.7341\times 10^{-8}$
$3D$ with $\beta=0$ | $0.29012486$ | $1.39401840\times 10^{-2}$ | $-4.9300\times 10^{-10}$ | $-9.2410\times 10^{-10}$
$1.4$ | $2D$ | $0.22044748$ | $1.0555690\times 10^{-3}$ | $-2.4435\times 10^{-9}$ | $-4.6050\times 10^{-9}$
$3D$ with $\beta=0$ | $0.22088128$ | $3.2584926\times 10^{-3}$ | $1.1763\times 10^{-9}$ | $-9.2065\times 10^{-10}$
$1.8$ | $2D$ | $0.18988059$ | $7.01090500\times 10^{-4}$ | $-2.3585\times 10^{-8}$ | $-9.6791\times 10^{-9}$
$3D$ with $\beta=0$ | $0.19109581$ | $1.28319800\times 10^{-3}$ | $3.6364\times 10^{-11}$ | $4.5576\times 10^{-10}$
$2.0$ | $2D$ | $0.17793112$ | $6.38880180\times 10^{-4}$ | $7.0879\times 10^{-9}$ | $-8.7970\times 10^{-9}$
$3D$ with $\beta=0$ | $0.17952602$ | $8.81986290\times 10^{-4}$ | $1.5717\times 10^{-9}$ | $-1.1064\times 10^{-9}$
Table 2: The spatial eigenvalues corresponding to the $2D$ and $3D$ case with $\beta=0$ tabulated for $M=0.6$, $1.4$, $1.8$ and $2$ when $Re=2500$ and (a) $\omega_{r}=0.06$ and (b) $\omega_{r}=0.1$. The corresponding values of the real and imaginary part of the wall dispersion function $Y_{j}(0)$ are also enlisted, where $j=7$ and $23$ for $2D$ and $3D$ cases, respectively. Formulation | $M=0.6$ | $M=1.8$ | $M=2$
---|---|---|---
$Re_{cr}$ | $\sqrt{Re_{x_{cr}}}$ | $Re_{cr}$ | $\sqrt{Re_{x_{cr}}}$ | $Re_{cr}$ | $\sqrt{Re_{x_{cr}}}$
$2D$ | $566.35$ | $303.51$ | $1360.70$ | $458.15$ | $2043.58$ | $630.73$
$3D$ with $\beta=0$ | $588.12$ | $315.17$ | $1089.19$ | $366.73$ | $1829.06$ | $564.52$
$\beta=0.1$ | $617.32$ | $330.82$ | $1089.05$ | $366.68$ | $1256.57$ | $387.83$
$\beta=0.3$ | $1214.11$ | $650.65$ | $990.23$ | $333.41$ | $1022.91$ | $315.71$
$\beta=0.4$ | $-$ | $-$ | $1223.46$ | $411.93$ | $1144.14$ | $353.12$
$\beta=0.45$ | $-$ | $-$ | $1533.75$ | $516.41$ | $1329.31$ | $410.28$
$\beta=0.5$ | $-$ | $-$ | $2392.37$ | $805.51$ | $1742.46$ | $537.79$
Table 3: $Re_{cr}$ based on displacement thickness and local streamwise
coordinate is tabulated for $2D$ compressible boundary layer corresponding to
indicated disturbance cases.

In Fig. 7(a-d), we show the neutral curve in $(Re,\omega_{r})$-plane
corresponding to $2D$ ($6^{th}$-order system) and $3D$ case with $\beta=0$
($8^{th}$-order system) for $M=0.6$, $1.4$, $1.8$, and $2$, respectively. In
Table 2, we show the values of the eigenvalue $\alpha$ for $Re=2500$ and
$\omega_{r}=0.06$ and $0.1$, respectively, for the above cases. The first
point, i.e., $Re=2500$ and $\omega_{r}=0.06$ is inside the neutral curve,
while the second point, i.e., $Re=2500$ and $\omega_{r}=0.1$ is in the
spatially stable region for all the cases. The converged value of the
dispersion function $Y_{j}(0)$ is also noted in the table, where $j=7$ and
$23$ for $2D$ and $3D$ cases, respectively. Figure 7 shows that for the
subsonic to low supersonic Mach numbers (up to at least $M=1.4$), the unstable
zone corresponding to the $2D$ case is larger along with a slightly lower
value of $Re_{cr}$. A tentative reason may be provided following Mack[3] that
the viscous modes dominate for these cases, and as non-zero perturbation
spanwise velocity $\hat{w}$ needs to extract energy from the mean flow for
sustenance, the stability of the flow is enhanced. However, for $M=1.8$ and
$2$ cases, the $3D$ case with $\beta=0$ is more spatially unstable than the
corresponding $2D$ case. This is also noted for the corresponding $\alpha$’s
tabulated in table 2 for $Re=2500$ and $\omega_{r}=0.06$. A possible reason
for such behavior may be the dominance of apparently inviscid mode for higher
Mach numbers as pointed in Mack[3], and Lees & Reshotko [52]. These points
need further detailed investigation, possibly by comparing the wall-normal
eigenstructure.

We have also investigated the effect of $\beta$ for a subsonic Mach number of
$M=0.6$ (results not presented here). For the subsonic case, stability of the
flow increases with an increase in spanwise wavenumber both in terms of
reduction in $Re_{cr}$ and the extent of the spatially unstable zone in
$(Re,\omega_{r})$-plane. In Table 3, we tabulate the critical Reynolds number
$Re_{cr}$ based on displacement thickness and local streamwise coordinate for
$M=0.6,1.8$ and $2$. It is often concluded that a $2D$ mean flow is more
stable for $3D$ disturbances than pure $2D$ disturbances (see i.e., Drazin &
Reid [15]). While this conclusion is valid for incompressible boundary
layers[53] and also for subsonic compressible adiabatic boundary layer cases,
as noted here for $M=0.6$, Figs. 6 and 7 and Tables 2 and 3 prove that, it
strictly does not hold for supersonic boundary layer cases.

###  5.4 Spatial stability corresponding to $2D$ disturbances for $M>3$

Next, we investigate the spatial instability of compressible boundary layers
for $M>3$. Mack[3, 23] reported the appearance of the second unstable mode at
finite Reynolds number for $M>3$, which is described to be inviscid in nature
due to the presence of generalized inflection point
$y_{s}=(U^{\prime}/T^{\prime})^{\prime}$ at certain height. The neutral curve
for the second mode is at the top of that corresponding to the first viscous
mode in Mack [3] for $3<M<4.5$ insulated wall-bounded boundary layer. As $M$
increases, the amplification rate of the second mode increases while the
corresponding neutral curve tends to merge with that for the viscous first
mode. The merger between the first and the second mode happens according to
Mack[3] beyond $M=4.6$ for the adiabatic wall-bounded boundary layer. The
existence of the high-frequency second mode is also reported in Kendall[25],
Stetson[28], Fedorov & Tumin [29] and other references contained therein. The
above observations motivate us to relook at the $2D$ spatial instability of
the supersonic boundary layers beyond $M>3$ by CMM.

Figure 8: Spatial $2D$ modes in $(\alpha_{r},\alpha_{i})$-plane for $Re=1000$ and $\omega_{r}=0.1$ and selected neutral curves shown for (a,b) $M=4$, (c,d) $M=3$ and (e,f) $M=6$ for $2D$ compressible boundary layer over an insulated plate. | Mode No. | $\alpha_{r}$ | $\alpha_{i}$ | $Y_{7r}(0)$ | $Y_{7i}(0)$
---|---|---|---|---|---
(a) $Re=1000$, | $1$ | $0.13099846$ | $2.61582640\times 10^{-3}$ | $-6.4153\times 10^{-9}$ | $-6.9270\times 10^{-9}$
and $\omega_{r}=0.1$ | $2$ | $0.15160470$ | $2.94952140\times 10^{-3}$ | $-9.7098\times 10^{-9}$ | $-2.2834\times 10^{-9}$
$3$ | $0.18859601$ | $2.67518660\times 10^{-3}$ | $-2.1216\times 10^{-9}$ | $-1.9402\times 10^{-8}$
$4$ | $0.22759297$ | $2.65627840\times 10^{-3}$ | $8.5090\times 10^{-9}$ | $2.6440\times 10^{-9}$
$5$ | $0.26758295$ | $2.79169690\times 10^{-3}$ | $-1.2911\times 10^{-8}$ | $9.8854\times 10^{-10}$
$6$ | $0.30832326$ | $3.03706200\times 10^{-3}$ | $-1.1347\times 10^{-8}$ | $2.6369\times 10^{-8}$
$7$ | $0.34962437$ | $3.35736340\times 10^{-3}$ | $-2.7888\times 10^{-8}$ | $-3.6595\times 10^{-9}$
$8$ | $0.39132136$ | $3.73667760\times 10^{-3}$ | $-3.7469\times 10^{-8}$ | $2.0684\times 10^{-9}$
(b) $Re=15000$, | $1$ | $0.19498196$ | $1.93222380\times 10^{-4}$ | $-1.0033\times 10^{-8}$ | $7.4700\times 10^{-9}$
and $\omega_{r}=0.15$ | $2$ | $0.21168292$ | $-2.76303440\times 10^{-4}$ | $5.4583\times 10^{-9}$ | $-3.6451\times 10^{-9}$
$3$ | $0.24550970$ | $-7.74741060\times 10^{-4}$ | $-4.8893\times 10^{-10}$ | $9.6411\times 10^{-10}$
$4$ | $0.28347549$ | $-1.00149330\times 10^{-3}$ | $-5.3637\times 10^{-8}$ | $7.0336\times 10^{-9}$
$5$ | $0.32309651$ | $-1.08311510\times 10^{-3}$ | $1.1205\times 10^{-9}$ | $7.0153\times 10^{-9}$
$6$ | $0.36358333$ | $-1.06486810\times 10^{-3}$ | $-5.1550\times 10^{-9}$ | $-4.2671\times 10^{-9}$
$7$ | $0.40459979$ | $-9.47628640\times 10^{-4}$ | $-5.1089\times 10^{-9}$ | $-1.1243\times 10^{-8}$
$8$ | $0.44597900$ | $-7.45557130\times 10^{-4}$ | $-2.8856\times 10^{-8}$ | $2.5628\times 10^{-8}$
(c) $Re=29000$, | $1$ | $0.19497371$ | $-2.6368667\times 10^{-5}$ | $-2.0988\times 10^{-9}$ | $-4.6073\times 10^{-10}$
and $\omega_{r}=0.15$ | $2$ | $0.21163477$ | $-3.9839061\times 10^{-4}$ | $-1.6398\times 10^{-8}$ | $-2.3914\times 10^{-8}$
$3$ | $0.24537614$ | $-8.3571899\times 10^{-4}$ | $-6.9107\times 10^{-9}$ | $2.6590\times 10^{-9}$
$4$ | $0.28324479$ | $-9.1820088\times 10^{-4}$ | $1.2776\times 10^{-8}$ | $1.6192\times 10^{-8}$
$5$ | $0.32281691$ | $-7.9247856\times 10^{-4}$ | $-1.9835\times 10^{-9}$ | $-2.8773\times 10^{-9}$
$6$ | $0.36329904$ | $-5.8810069\times 10^{-4}$ | $6.9281\times 10^{-10}$ | $-3.2294\times 10^{-9}$
$7$ | $0.40427428$ | $-3.4727529\times 10^{-4}$ | $9.0748\times 10^{-9}$ | $2.0358\times 10^{-9}$
$8$ | $0.44549695$ | $-5.2784686\times 10^{-5}$ | $-1.6727\times 10^{-8}$ | $-2.0506\times 10^{-8}$
Table 4: The spatial eigenvalues corresponding to the indicated modes for
$M=4$ obtained after Newton-Raphson polishing tabulated for (a)
$(Re,\omega_{r})=(1000,0.1)$, (b) $(Re,\omega_{r})=(15000,0.15)$, and (c)
$(Re,\omega_{r})=(29000,0.15)$ for the $2D$ compressible boundary layer over
insulated wall. The corresponding values of the real and imaginary parts of
$Y_{7}(0)$ are also enlisted after convergence is achieved.

Figure 8(a) shows the location of spatial modes for $M=4$, $Re=1000$ and
$\omega_{r}=0.1$ in $(\alpha_{r},\alpha_{i})$-plane. The mode-$1$ is viscous
in nature, while on the right of that, we note a train of modes with
increasing $\alpha_{r}$. Here, we number these modes as $2-8$. In the figure,
the range of $\alpha_{r}$ is taken up to $0.4$. As we increase the range of
$\alpha_{r}$, more and more such modes are noted to exist. Modes $1-8$ for
$Re=1000$ and $\omega_{r}=0.1$ are also tabulated in Table 4(a). Figure 8(b)
shows the neutral curves for modes $1-6$. The shape of the neutral curves for
modes-$7$ and $8$ are similar to modes $3$ and $8$ (not shown here). The
$Re_{cr}$ for mode-$1$ is considerably higher than that of other modes. While
$Re_{cr}$ is noted to be approximately $26051$ and $6500$ for modes $1$ and
$2$, respectively, $Re_{cr}\simeq 3500$ corresponding to modes $3-8$. From the
neutral curve for mode-$2$ onwards, we also note that these modes are
spatially unstable over a broader range of frequencies than mode-$1$,
subsonic, and low supersonic Mach number ($M\leq 2$) cases. For the last
cases, disturbances corresponding to $\omega_{r}>0.15$ are noted to be
spatially stable. From mode-$2$ onwards, the maximum limit of $\omega_{r}$ up
to which spatial instability may occur increases from $0.4$ to approximately
$0.7$ for $M=4$. We also note the upper branch of the neutral curve for modes
$4$, $5$, and $6$ have a dropping tendency at high Reynolds numbers. Whether
these curves drop continuously and subsequently close at certain higher $Re$
will be investigated in the future. We also enlist spatial eigenvalues
corresponding to modes $1-8$ obtained after Newton-Raphson polishing in Table
4 for $Re=15000$ and $29000$, respectively, when $\omega_{r}=0.15$ for
reference. The first point is inside the spatially unstable zone corresponding
to modes $2-8$ but outside the neutral curve for mode-$1$. The second point is
within the unstable zone corresponding to modes $1-8$. We note by comparing
the values of $\alpha_{i}$ from Tables 1-4 that values of $\alpha_{i}$’s
corresponding to spatial instability for $M=4$ are roughly one order of
magnitude lower than subsonic or low supersonic ($M<2$) cases.

Figure 9: Contours of $\alpha_{r}$ and group velocity plotted in
$(Re,\omega_{r})$-plane for Modes-$1$, $3$ and $6$ corresponding to $M=4$. The
corresponding neutral curve is also shown for reference.

The number of such modes increases with the increase in Mach number for $M>3$.
This can be concluded by comparing Figs. 8(c,e) with Fig. 8(a). We note only
$6$ modes up to $\alpha_{r}=0.4$ for $M=3$, whereas $8$ and $10$ modes are
noted within this range for $M=4$ and $6$, respectively. All these modes
depict spatial instability beyond respective critical $Re$. The $Re_{cr}$ for
mode-$1$ corresponding to $M=4$ and $6$ is $26,051$ and $25305$, respectively
while that for $M=3$ is approximately $6.52\times 10^{5}$ (results not shown
here). The selected neutral curves for $M=3$ and $6$ are shown in Figs. 8(d)
and 8(f), respectively, which show that the maximum limit for $\omega_{r}$ up
to which spatially unstable disturbances are noted, increases from $M=3$ to
$6$. Therefore, in summary, for $M>3$, the number of unstable modes increases
drastically with the Mach number. This also coincides with the increase in the
maximum limit for $\omega_{r}$ for which spatially unstable perturbations are
observed.

It is pertinent to note that the shape of the neutral curves for mode-$2$
corresponding to $M=4$ and $6$ resembles the fused neutral curve shown for
$M=4.8$ in Mack[3]. Except for this similarity, the nature and number of
spatially unstable zones obtained here by CMM are different than that reported
in Mack[3] and other subsequent references to date for the compressible
boundary layer. To the best of our knowledge, the existence of such higher-
order modes from viscous calculations for supersonic compressible boundary
layers has yet to be reported in the literature so far. This is one of the
novel aspects of the present investigation.

Next, we illustrate the nature of the spatial modes for higher supersonic Mach
number cases. Figure 9 shows contours of $\alpha_{r}$ and group velocity
$V_{g}=\partial\omega_{r}/\partial\alpha_{r}$ in $(Re,\omega_{r})$-plane for
modes-$1$, $3$ and $6$ corresponding to $M=4$. We note that $\alpha_{r}$ and
$V_{g}$ are almost independent of $Re$, increasing with the increase in
$\omega_{r}$ for mode-$1$. Almost similar variation in $\alpha_{r}$ is also
noted for modes-$3$ and $6$; however, for these modes, a comparatively higher
rate of increase in $\alpha_{r}$ with $\omega_{r}$ is observed. The group
velocity $V_{g}$ is almost independent of $Re$ for moderate to high-frequency
cases ($\omega_{r}\gtrsim 0.1$) corresponding to these modes similar to
mode-$1$. However, in contrast to mode-$1$, $V_{g}$ decreases with an increase
in $\omega_{r}$, indicating signals corresponding to these modes tend to move
slower with increased frequency. The behavior displayed by modes-$1$, $3$, and
$6$ is also observed for other spatially unstable modes for $M=3$, $4$, and
$6$ (results not shown here). Following these variations, tentatively, we say
that $\alpha_{r}\simeq C\omega_{r}^{q}$, where the index $q<1$ for mode-$1$
and $q>1$ for higher order modes (mode-$2$ onwards). Considering all the modes
collectively, we also can say that for $M=4$, the spatially unstable
disturbances travel downstream at speeds varying from $0.78$ to $0.89$ times
the free-stream speed. These values are considerably higher than those
compared to subsonic cases, where low to moderate frequency disturbances are
seen to travel at speeds not exceeding $0.5$ times $U_{\infty}$.

##  6 Summary and conclusions

We revisit the spatial linear stability analysis of the $2D$ compressible
boundary layer over an adiabatic flat plate to $2D$ and $3D$ disturbances. We
apply the compound matrix method (CMM) in this regard to remove the stiffness
of the problem contrary to conventional methods using Gram-Schmidt ortho-
normalization or discretizing the governing equations by appropriate finite
difference schemes. The CMM has not been employed before to perform the
stability analysis of compressible subsonic or supersonic boundary layers for
$2D$ or $3D$ perturbations. While the former yields a $6^{th}$-order system,
the latter results in a $8^{th}$-order system governing the wall-normal
variation of spectral amplitudes as illustrated in Sec. 2.1. The corresponding
methodology to determine the eigenvalues is described in Sec. 3. The
eigenvalues representing complex wavenumber $\alpha=\alpha_{r}+i\alpha_{i}$
for a particular combination of Reynolds number $Re$, Mach number $M$,
frequency $\omega_{r}$ and spanwise wavenumber $\beta$ are found by
numerically integrating auxiliary compound matrix equations from free stream
to wall ($y=0$) and satisfying corresponding dispersion relation. Spatial
instability is determined when $\alpha_{i}<0$.

We consider the Mach number of the flow ranging from $0.1$ (subsonic regime)
to $6$ (supersonic regime) in the analysis presented here. We first validate
the methodology by comparing the spatial stability of the compressible $M=0.1$
boundary layer with the corresponding incompressible Blasius shear layer. The
variation of the spectral amplitudes for the latter is given by the Orr-
Sommerfeld equation (OSE). We show an excellent match between the spatial
stability characteristics for the primary and secondary modes corresponding to
both flows. Subsequently, we analyze the effects of flow compressibility (Mach
number $M$) and effects of spanwise wavenumber $\beta$. We note that flow
becomes increasingly stable with an increase in $M$ up to $M=2$. When
considering $3D$ disturbances on the $2D$ compressible mean flow, we also note
that for subsonic cases, flow is increasingly stable with an increase in
spanwise wavenumber $\beta$, similar to the incompressible Blasius boundary
layer. In contrast, there exists an optimum value for $\beta$ up to which
spatial stability reduces with an increase in $\beta$ for supersonic cases up
to $M=2$. Subsequent increase in $\beta$ enhances flow stability for these
cases. We also find that spatial stability of $2D$ flows to pure $2D$ and $3D$
disturbance with no spanwise variation (i.e., $\beta=0$) are not identical.
While the latter indicates the perturbation wall-normal vorticity
$\hat{\xi}_{y}$ and spanwise velocity $\hat{w}$ are trivially zero, the former
assumes that no perturbation quantities are varying along spanwise directions
as illustrated in Sec. 5.3. While the former is more spatially unstable up to
$M\simeq 1.4$, we note the opposite scenario for $1.4<M<2$. Only one primary
unstable mode is found for these cases up to $M=2$; other modes are always
spatially stable.

We next investigate the spatial stability for $M>3$ corresponding to $2D$
disturbances. Mack[3] reported the existence of two unstable zones for such
cases, which fuse to form one unstable zone with larger frequency extent at
higher Mach numbers and attributed such characteristics of supersonic flows to
enhanced inviscid acoustic instability. We note a series of unstable modes for
$M>3$, and the number of such modes is much more than two, as reported in
Mack[3]. The number and the frequency extent of the unstable zone for these
modes increase significantly with an increase in Mach number. For example, we
note $6$, $8$ and $10$ spatially unstable modes up to streamwise wavenumber of
$\alpha_{r}=0.4$ for $M=3$, $4$ and $6$, respectively. These modes display
spatial instability beyond certain critical Reynolds number $Re_{cr}$. We
calculate these spatially unstable modes to travel at a much higher speed of
$0.78$ to $0.89$ times the free-stream speed than those corresponding to
incompressible, subsonic, and low supersonic ($M<2$) cases. While the shape of
the neutral curves for the second unstable mode for $M=4$ and $6$ bears
similarity to the fused neutral curve shown in Mack[3] for $M=4.8$, the
characteristics of spatially unstable higher-order modes, to the best of our
knowledge, have not been shown or reported in any related literature so far
considering viscous stability of supersonic boundary layer.

##  7 Acknowledgments

SB acknowledges the support provided by DST(SERB) under MATRICS
(MTR/2020/000568) and CRG (CRG/2021/006377) schemes.

## References

  * Xu  _et al._ [2022] J. Xu, Y. Wang, N. Liu, J. Yang, Y. Mu, W. Liang, J. Zhang, and J. Bai, “An improved physics-informed transition-turbulence model for asymmetric transition over supersonic rotating projectiles,” Computers & Fluids 238, 105366 (2022).
  * Lee and Jiang [2019] C. Lee and X. Jiang, “Flow structures in transitional and turbulent boundary layers,” Physics of Fluids 31 (2019).
  * Mack [1984a] L. M. Mack, “Boundary-layer linear stability theory,” Tech. Rep. (California Inst of Tech Pasadena Jet Propulsion Lab, 1984).
  * Reed, Saric, and Arnal [1996] H. L. Reed, W. S. Saric, and D. Arnal, “Linear stability theory applied to boundary layers,” Annual Review of Fluid Mechanics 28, 389–428 (1996).
  * Dettenrieder and Bodony [2022] F. Dettenrieder and D. J. Bodony, “Stability analyses of compressible flat plate boundary layer flow over a mechanically compliant wall,” Theoretical and Computational Fluid Dynamics 36, 141–153 (2022).
  * Bitter and Shepherd [2015] N. Bitter and J. Shepherd, “Stability of highly cooled hypervelocity boundary layers,” Journal of Fluid Mechanics 778, 586–620 (2015).
  * Schneider [2001] S. P. Schneider, “Effects of high-speed tunnel noise on laminar-turbulent transition,” Journal of Spacecraft and Rockets 38, 323–333 (2001).
  * Malik [1990] M. R. Malik, “Numerical methods for hypersonic boundary layer stability,” Journal of Computational Physics 86, 376–413 (1990).
  * Kosinov  _et al._ [2015] A. Kosinov, A. Panina, N. Semionov, and Y. G. Yermolaev, “On the relative “receptivity” of two and three-dimensional supersonic boundary layers to stationary disturbances at Mach 2,” Progress in Flight Physics 7, 509–524 (2015).
  * Fedorov [2011] A. Fedorov, “Transition and stability of high-speed boundary layers,” Annual Review of Fluid Mechanics 43, 79–95 (2011).
  * Zhong and Wang [2012] X. Zhong and X. Wang, “Direct numerical simulation on the receptivity, instability, and transition of hypersonic boundary layers,” Annual Review of Fluid Mechanics 44, 527–561 (2012).
  * Unnikrishnan and Gaitonde [2019] S. Unnikrishnan and D. V. Gaitonde, “Interactions between vortical, acoustic and thermal components during hypersonic transition,” Journal of Fluid Mechanics 868, 611–647 (2019).
  * Reshotko [2008] E. Reshotko, “Transition issues for atmospheric entry,” Journal of Spacecraft and Rockets 45, 161–164 (2008).
  * Tollmien [1931] W. Tollmien, “The Production of Turbulence,” NACA Report-TM-609 (1931).
  * Drazin and Reid [2004] P. G. Drazin and W. H. Reid,  _Hydrodynamic stability_ (Cambridge university press, 2004).
  * Tollmien [1936] W. Tollmien, “General instability criterion of laminar velocity distributions,” Tech. Rep. (1936).
  * Schlichting [1933] H. Schlichting, “Zur Entstehung der Turbulenz bei der Plattenstr$\ddot{o}$mung,” Nachr. Ges. Wiss. G$\ddot{o}$ttingen, Math. Phys. Kl. 42, 181–208 (1933).
  * Heisenberg [1924] W. Heisenberg, “$\ddot{U}$ber Stabilit$\ddot{a}$t und turbulenz von Flussigkeitsstr$\ddot{o}$men,” Ann. Phys. Lpz. 74, 577–627 (1924).
  * Kelvin [1871] L. Kelvin, “Hydrokinetic solutions and observations,” Phil. Mag. (4) 42, 362–377 (1871).
  * Healey [2006] J. Healey, “Inviscid long-wave theory for the absolute instability of the rotating-disc boundary layer,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 462, 1467–1492 (2006).
  * Schubauer and Skramstad [1947] G. B. Schubauer and H. K. Skramstad, “Laminar boundary-layer oscillations and stability of laminar flow,” Journal of the Aeronautical Sciences 14, 69–78 (1947).
  * Mack [1975] L. M. Mack, “Linear stability theory and the problem of supersonic boundary-layer transition,” AIAA Journal 13, 278–289 (1975).
  * Mack [1984b] L. M. Mack, “Remarks on disputed numerical results in compressible boundary-layer stability theory,” Physics of Fluids 27, 342–347 (1984b).
  * Gasperas [1989] G. Gasperas, “Effect of wall temperature distribution on the stability of the compressible boundary layer,” in  _20th Fluid Dynamics, Plasma Dynamics and Lasers Conference_ (1989) p. 1894.
  * Kendall [1975] J. M. Kendall, “Wind tunnel experiments relating to supersonic and hypersonic boundary-layer transition,” AIAA Journal 13, 290–299 (1975).
  * Demetriades [1974] A. Demetriades, “Hypersonic viscous flow over a slender cone. III-Laminar instability and transition,” in  _$7^{th}$ Fluid and Plasmadynamics Conference_ (1974) p. 535.
  * Stetson  _et al._ [1983] K. Stetson, E. Thompson, J. Donaldson, and L. Siler, “Laminar boundary layer stability experiments on a cone at Mach $8$. I-Sharp cone,” in  _$16^{th}$ Fluid and Plasmadynamics Conference_ (1983) p. 1761.
  * Stetson [1992] K. F. Stetson, “Hypersonic boundary-layer transition,” in  _Advances in Hypersonics: Defining the Hypersonic Environment_ , Vol. 1 (Springer, 1992) pp. 324–417.
  * Fedorov and Tumin [2003] A. Fedorov and A. Tumin, “Initial-value problem for hypersonic boundary-layer flows,” AIAA Journal 41, 379–389 (2003).
  * Fedorov  _et al._ [2003] A. Fedorov, A. Shiplyuk, A. Maslov, E. Burov, and N. Malmuth, “Stabilization of a hypersonic boundary layer using an ultrasonically absorptive coating,” Journal of Fluid Mechanics 479, 99–124 (2003).
  * Bertolotti and Herbert [1991] F. Bertolotti and T. Herbert, “Analysis of the linear stability of compressible boundary layers using the PSE,” Theoretical and Computational Fluid Dynamics 3, 117–124 (1991).
  * Malik [1982] M. R. Malik, “Finite-difference solution of the compressible stability eigenvalue problem,” Tech. Rep. (NASA, 1982).
  * Mack [1976] L. M. Mack, “A numerical study of the temporal eigenvalue spectrum of the Blasius boundary layer,” Journal of Fluid Mechanics 73, 497–520 (1976).
  * Özgen and Kırcalı [2008] S. Özgen and S. A. Kırcalı, “Linear stability analysis in compressible, flat-plate boundary-layers,” Theoretical and Computational Fluid Dynamics 22, 1–20 (2008).
  * Ng and Reid [1985] B. Ng and W. Reid, “The compound matrix method for ordinary differential systems,” Journal of Computational Physics 58, 209–228 (1985).
  * Yiantsios and Higgins [1988] S. G. Yiantsios and B. G. Higgins, “Numerical solution of eigenvalue problems using the compound matrix method,” Journal of Computational Physics 74, 25–40 (1988).
  * Allen and Bridges [2002] L. Allen and T. J. Bridges, “Numerical exterior algebra and the compound matrix method,” Numerische Mathematik 92, 197–232 (2002).
  * Sengupta [2012] T. K. Sengupta,  _Instabilities of flows and transition to turbulence_ (CRC Press, 2012).
  * Sengupta, Bhaumik, and Bose [2013] T. Sengupta, S. Bhaumik, and R. Bose, “Direct numerical simulation of transitional mixed convection flows: Viscous and inviscid instability mechanisms,” Physics of Fluids 25 (2013).
  * Anturkar, Papanastasiou, and Wilkes [1990] N. R. Anturkar, T. C. Papanastasiou, and J. O. Wilkes, “Linear stability analysis of multilayer plane Poiseuille flow,” Physics of Fluids A: Fluid Dynamics 2, 530–541 (1990).
  * Gubernov, Sidhu, and Mercer [2006] V. Gubernov, H. S. Sidhu, and G. N. Mercer, “Generalized compound matrix method,” Applied Mathematics Letters 19, 458–463 (2006).
  * Haughton [2008] D. M. Haughton, “Evaluation of eigenfunctions from compound matrix variables in non-linear elasticity–II. Sixth order systems,” Journal of Computational Physics 227, 8960–8967 (2008).
  * Stewartson [1964] K. Stewartson,  _The theory of laminar boundary layers in compressible fluids_ , Vol. 1 (Clarendon Press Oxford, 1964).
  * Kundu, Cohen, and Dowling [2015] P. K. Kundu, I. M. Cohen, and D. R. Dowling,  _Fluid Mechanics_ (Academic Press, 2015).
  * Schmid and Henningson [2001] P. J. Schmid and D. S. Henningson,  _Stability and transition in shear flows. Applied Mathematical Sciences 142_ (Springer Nature, 2001).
  * Barker  _et al._ [2018] B. Barker, J. Humpherys, G. Lyng, and J. Lytle, “Evans function computation for the stability of travelling waves,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20170184 (2018).
  * Shapiro [1954] A. H. Shapiro,  _The Dynamics and Thermodynamics of Compressible Fluid Flow_ , Vol. II (Ronald Press, New York, 1954).
  * White [2006] F. M. White,  _Viscous fluid flow_ , Vol. 3 (McGraw-Hill New York, 2006).
  * Hilsenrath [1955] J. Hilsenrath,  _Tables of thermal properties of gases: comprising tables of thermodynamic and transport properties of air, argon, carbon dioxide, carbon monoxide, hydrogen, nitrogen, oxygen, and steam_ , Vol. 564 (US Department of Commerce, National Bureau of Standards, 1955).
  * Sengupta, Nair, and Rana [1997] T. Sengupta, M. Nair, and V. Rana, “Boundary layers excited by low frequency disturbances—Klebanoff mode,” Journal of Fluids and Structures 11, 845–853 (1997).
  * Dunn and Lin [1955] D. W. Dunn and C.-C. Lin, “On the stability of the laminar boundary layer in a compressible fluid,” Journal of the Aeronautical Sciences 22, 455–477 (1955).
  * Lees and Reshotko [1962] L. Lees and E. Reshotko, “Stability of the compressible laminar boundary layer,” Journal of Fluid Mechanics 12, 555–590 (1962).
  * Maddipati, Sengupta, and Sundaram [2021] R. Maddipati, T. K. Sengupta, and P. Sundaram, “Relevance of two-and three-dimensional disturbance field explained with linear stability analysis of Orr-Sommerfeld equation by compound matrix method,” Computers & Fluids 225, 104965 (2021).

Appendix-$I$: The formulation of the CMM

The equation for the disturbance spectral amplitudes, given by Eqs. (16-19),
are cast in the form of a system of first-order ODEs following Eq. (32), where
the state vector
$\left\\{\textbf{X}\right\\}=\left[X_{1},X_{2},X_{3},X_{4},X_{5},X_{6},X_{7},X_{8}\right]^{T}$
are defined as $X_{1}=\psi$, $X_{2}=\psi^{\prime}$, $X_{3}=\varphi$,
$X_{4}=\varphi^{\prime}$, $X_{5}=\Theta$, $X_{6}=\Theta^{\prime}$,
$X_{7}=\Omega$ and $X_{8}=\Omega^{\prime}$ for $3D$ disturbances. The elements
of the matrix $[E]$ are obtained from Eqs. (16-19). We can express Eq. (32)
using tensorial index notation as

$\displaystyle X_{j}=E_{jk}X_{k}$ (49)

There are four modes decaying in the free-stream for $3D$ disturbance, and one
can express the state vector as a linear combination of the four decaying
physical modes, following Eq. (29) as

$\displaystyle
X_{j}^{\prime}=c_{1}X_{j,1}+c_{3}X_{j,3}+c_{5}X_{j,5}+c_{7}X_{j,7}$ (50)

where $j=1,\cdots,8$. The auxiliary compound variables for CMM are defined
following Eq. (34). While defining the compound variables, we need to maintain
the following order among the indices as $1\leq j<k<l<m\leq 8$, so that there
is no repetition of the basis row-vectors
$h_{j}=[X_{j,1},X_{j,3},X_{j,5},X_{j,7}]$. We have denoted the compound
variables as
$Y_{n}=\mathcal{Y}_{j,k,l,m}=\det\left([h_{j},h_{k},h_{l},h_{m}]^{T}\right)$,
where $\det(\cdot)$ represents the determinant of a square matrix. For given
$j$, $k$, $l$, and $m$ indices corresponding to the state vector $X$, the
index $n$ of the compound variable $Y_{n}$ can be obtained by sequentially
counting it such that $1\leq j<k<l<m\leq 8$. Following this methodology, we
can mathematically express $n$ as a function of $j$, $k$, $l$, and $m$ indices
as

$\displaystyle n$ $\displaystyle=$
$\displaystyle\sum_{j_{1}=1}^{j-1}\sum_{k_{1}=j_{1}+1}^{6}\sum_{l_{1}=k_{1}+1}^{7}\sum_{m_{1}=l_{1}+1}^{8}\left(1\right)+\sum_{k_{1}=j+1}^{k-1}\sum_{l_{1}=k_{1}+1}^{7}\sum_{m_{1}=l_{1}+1}^{8}\left(1\right)$
(51)
$\displaystyle+\sum_{l_{1}=k+1}^{l-1}\sum_{m_{1}=l_{1}+1}^{8}\left(1\right)+\sum_{m_{1}=l+1}^{m-1}\left(1\right)+1$

The compound variable index $n$ by evaluating the above series is given as
$n=(m-l)+(l-k-1)(16-l-k)/2+(k-j-1)(168-23k-22j+k^{2}+j^{2}+kj)/6+(j-1)(-j^{3}+29j^{2}-306j+1344)/24$.
The inverse calculation, i.e., for a given index $n$ of the compound variable
$Y_{n}$, the indices $j$, $k$, $l$, and $m$ of the constituent basis row-
vectors can also be obtained from the above algorithm or the formulae given in
Eq. (51).

Appendix-$II$: Derivatives of the compound variable

Here, we explain the methodology involved in finding the first derivative of
the $n^{th}$ compound variable $Y_{n}$ so that elements of $[F]$-matrix (see
Eq. (35)) can be evaluated. It is quite easy to explicitly evaluate the
coefficients of the $[F]$ matrix for the Orr-Sommerfeld equation[15],
primarily because the $[E]$ matrix which relates $\textbf{X}^{\prime}$ with X
(see Eq. (32)) is significantly sparse in nature. For the present case, the
$[E]$ matrix is relatively dense for $6^{th}$\- or $8^{th}$-order cases
corresponding to $2D$ and $3D$ disturbances. Therefore, we develop an
algorithm based on several well-known theorems on the properties of the
determinant to circumvent the manual calculation of the coefficients of the
$[F]$ matrix. The proposed algorithm is generic and can be used for even
higher-order systems if one considers chemical reactions, non-equilibrium
thermodynamics, real-gas effects, or electro-magnetodynamic effects. This
algorithm is illustrated next.

First, we enumerate several well-known theorems on determinants that are used
to derive the auxiliary equations following CMM. Let $\det(\textbf{A})$
represents the determinant of $n\times n$ square matrix
$\textbf{A}=[\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{i},\cdots,\textbf{a}_{n}]^{T}$,
where, $\textbf{a}_{i}$ represents the $i^{th}$ row of A.

  1. 1.

The derivative of $\det(\textbf{A})$ is equal to the sum of $n$ auxiliary
determinants of $n\times n$ square matrices $\textbf{A}_{i}$ such that the
matrix $\textbf{A}_{i}$ is equal to A except for the $i^{th}$-row, where the
elements of A are replaced by their respective derivatives. Mathematically,
$\det(\textbf{A})^{\prime}=\sum_{i=1}^{n}\det(\textbf{A}_{i})$, where
$\textbf{A}_{i}=\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{i}^{\prime},\cdots,\textbf{a}_{n}]^{T}\right)$.

  2. 2.

If two rows of the matrix A are interchanged, the value of the determinant is
multiplied by the factor of $(-1)$, i.e.,
$\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{i},\cdots,\textbf{a}_{j},\cdots,\textbf{a}_{n}]^{T}\right)=-\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{j},\cdots,\textbf{a}_{i},\cdots,\textbf{a}_{n}]^{T}\right)$.

  3. 3.

If the $i^{th}$-row of A is multiplied by a constant scalar $c$, then
$\det(\textbf{A})$ is multiplied by the same scalar $c$, i.e.,
$\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,c\textbf{a}_{i},\cdots,\textbf{a}_{n}]^{T}\right)=c\cdot\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{i},\cdots,\textbf{a}_{n}]^{T}\right)=c\cdot\det(\textbf{A})$.

  4. 4.

The determinant is multilinear, i.e., if the $i^{th}$-row of A is written as a
linear combination as $\textbf{a}_{i}=c_{1}\textbf{u}+c_{2}\textbf{v}$, where
u and v are row vectors and $c_{1}$ and $c_{2}$ are scalars, then
$\det(\textbf{A})$ can be expressied as a similar linear combination of two
determinants as
$\det(\textbf{A})=\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{a}_{i},\cdots,\textbf{a}_{n}]^{T}\right)=\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,c_{1}\textbf{u}+c_{2}\textbf{v},\cdots,\textbf{a}_{n}]^{T}\right)=c_{1}\cdot\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{u},\cdots,\textbf{a}_{n}]^{T}\right)+c_{2}\cdot\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{v},\cdots,\textbf{a}_{n}]^{T}\right)$.

  5. 5.

If two rows of A are identical, its determinant is zero, i.e.,
$\det\left([\textbf{a}_{1},\textbf{a}_{2},\cdots,\textbf{u},\cdots,\textbf{u},\cdots,\textbf{a}_{n}]^{T}\right)=0$.

We can determine wall-normal derivative of $Y_{n}$ using these theorems and
fact that $X_{j,q}^{\prime}=\sum_{l=1}^{8}E_{jl}X_{l,q}$, where $q=1,3,5,7$
represents the number corresponding to decaying physical modes in the free-
stream for the $8^{th}$-order system. Let,
$h_{j}=[X_{j,1},X_{j,3},X_{j,5},X_{j,7}]$ where $j=1,\cdots,8$ and hence,
$Y_{n}=\mathcal{Y}_{j,k,l,m}=\det\left([h_{j},h_{k},h_{l},h_{m}]^{T}\right)$,
where $1\leq j<k<l<m\leq 8$ following Eq. (34). Therefore,

$\displaystyle
Y_{n}^{\prime}=\mathcal{Y}_{j,k,l,m}^{\prime}=\underbrace{\begin{vmatrix}h_{j}^{\prime}\\\
h_{k}\\\ h_{l}\\\
h_{m}\end{vmatrix}}_{\text{Term-1}}+\underbrace{\begin{vmatrix}h_{j}\\\
h_{k}^{\prime}\\\ h_{l}\\\
h_{m}\end{vmatrix}}_{\text{Term-2}}+\underbrace{\begin{vmatrix}h_{j}\\\
h_{k}\\\ h_{l}^{\prime}\\\
h_{m}\end{vmatrix}}_{\text{Term-3}}+\underbrace{\begin{vmatrix}h_{j}\\\
h_{k}\\\ h_{l}\\\ h_{m}^{\prime}\end{vmatrix}}_{\text{Term-4}}$ (52)

We also know $h_{j}^{\prime}=E_{jp}h_{p}$, following Eq. (32). Therefore, we
can express terms $1-4$ of Eq. (52) by putting the above relation into it,
maintaining the order of the indices to define $Y_{n}$ and using theorems
$1-5$ on the determinants as

Term-1 $\displaystyle=$
$\displaystyle\left(\sum_{p=1}^{k-1}E_{jp}\mathcal{Y}_{p,k,l,m}-\sum_{p=k+1}^{l-1}E_{jp}\mathcal{Y}_{k,p,l,m}+\sum_{p=l+1}^{m-1}E_{jp}\mathcal{Y}_{k,l,p,m}-\sum_{p=m+1}^{8}E_{jp}\mathcal{Y}_{k,l,m,p}\right)\hskip
28.45274pt$ (53) Term-2 $\displaystyle=$
$\displaystyle\left(-\sum_{p=1}^{j-1}E_{kp}\mathcal{Y}_{p,j,l,m}+\sum_{p=j+1}^{l-1}E_{kp}\mathcal{Y}_{j,p,l,m}-\sum_{p=l+1}^{m-1}E_{kp}\mathcal{Y}_{j,l,p,m}+\sum_{p=m+1}^{8}E_{kp}\mathcal{Y}_{j,l,m,p}\right)\hskip
28.45274pt$ (54) Term-3 $\displaystyle=$
$\displaystyle\left(\sum_{p=1}^{j-1}E_{lp}\mathcal{Y}_{p,j,k,m}-\sum_{p=j+1}^{k-1}E_{lp}\mathcal{Y}_{j,p,k,m}+\sum_{p=k+1}^{m-1}E_{lp}\mathcal{Y}_{j,k,p,m}-\sum_{p=m+1}^{8}E_{lp}\mathcal{Y}_{j,k,m,p}\right)\hskip
28.45274pt$ (55) Term-4 $\displaystyle=$
$\displaystyle\left(-\sum_{p=1}^{j-1}E_{lp}\mathcal{Y}_{p,j,k,l}+\sum_{p=j+1}^{k-1}E_{lp}\mathcal{Y}_{j,p,k,l}-\sum_{p=k+1}^{l-1}E_{lp}\mathcal{Y}_{j,k,p,l}+\sum_{p=l+1}^{8}E_{lp}\mathcal{Y}_{j,k,l,p}\right)\hskip
28.45274pt$ (56)

The final expression for $Y_{n}^{\prime}$ is obtained by adding Eqs. (53-56).
These linear ODEs need to be integrated from the free-stream to the wall.
Therefore, the condition at the free-stream acts as an initial value for these
equations.