Instability in strongly stratified plane Couette flow with application to supercritical fluids

Abstract This paper addresses the stability of plane Couette flow in the presence of strong density and viscosity stratifications. It demonstrates the existence of a generalised inflection point that satisfies the generalised Fjørtoft criterion of instability when a minimum of kinematic viscosity is present in the base flow. The characteristic scales associated with this minimum are identified as the primary controlling parameters of the associated instability, regardless of the type of stratification. To support this finding, analytical stability models are derived in the long-wave approximation using piecewise linear base flows. Numerical stability calculations are carried out to validate these models and to provide further information on the production of disturbance vorticity. All instabilities are interpreted as arising from the interaction between two vorticity waves. Depending on the type of stratification, these two waves are produced by different physical mechanisms. When both strong density and viscosity stratifications are present, we show that they result from the concurrent action of shear and inertial baroclinic effects. The stability models developed for simple fluid models ultimately shed light on a recently observed unstable mode in supercritical fluids (Ren et al., J. Fluid Mech., vol. 871, 2019, pp. 831–864), providing a quantitative prediction of the stability diagram and identifying the dominant mechanisms at play. Furthermore, our study suggests that the minimum of kinematic viscosity reached at the Widom line in these fluids is the leading cause of their instability. The existence of similar instabilities in different fluids and flows (e.g. miscible fluids) is finally discussed.


Strongly stratified flows
The stability of stratified parallel shear flows, in which fluid properties such as density and viscosity vary in the direction perpendicular to that of the base flow, is a problem encountered in several geophysical phenomena (e.g., dynamics of the atmosphere and the ocean) and industrial applications.For example, fluids operating at pressure and temperature in the region of the critical point, which are employed in chemical and mechanical engineering (Brunner 2010;Liu et al. 2019), may exhibit large variations of viscosity and density in flows involving heat transfers.Stratified flows can be examined in different regimes.In this paper, the variations of fluid properties will not be restricted to small amplitudes, justifying the term of strongly stratified flow.Besides, gravity, and therefore buoyancy effects, will be neglected by assuming large Froude numbers; more details on the flow assumptions will be given in section 1.3.Density variations can nevertheless play a significant role in the flow dynamics through inertial effects.Menkes (1959) was perhaps the first to tackle the stability of such a flow, considering a hyperbolic tangent velocity profile with an exponential density stratification, which was found stabilising in this particular configuration.Soteriou & Ghoniem (1995) more comprehensively studied an incompressible mixing layer of two fluids of different densities.Depending on the density ratio, the primary instability was shown to exhibit either weaker or larger growth rates, to have its phase speed shifted, and its non-linear development altered.This last point was subsequently examined via secondary stability analysis (Reinaud et al. 2000;Fontane & Joly 2008) and direct numerical simulation (DNS) (Almagro et al. 2017).The mechanism responsible for the modified dynamics of this flow is the inertial baroclinic torque, which generates vorticity from misalignments between pressure and density gradients (Soteriou & Ghoniem 1995;Reinaud et al. 1999;Dixit & Govindarajan 2010).It is also at play in compressible flows (Lesshafft & Huerre 2007) but is classically neglected in buoyant flows modelled via the Boussinesq approximation, which ignores density variations in inertial terms (Drazin 1958;Guha & Raj 2018).
Strong viscosity stratifications will also be central in our problem, greatly affecting the base flow profile.Considering a parallel shear flow of two fluids of different viscosities separated by an interface, Yih (1967) showed that a long-wave instability exists at low Reynolds numbers.This instability does not require density gradients or surface tension effects: the jump in viscosity at the interface is sufficient to destabilise the flow.Hooper & Boyd (1983), in a similar configuration, revealed that a short-wave instability also grows at low Reynolds numbers.The mechanisms of these instabilities were discussed by Hinch (1984) and Charru & Hinch (2000).The effect of an interface of finite thickness was studied by Ern et al. (2003).The authors recovered the presence of low-Reynolds instabilities and furthermore showed that certain thicknesses could induce larger growth rates than an infinitely small one.Finally, another viscous instability exists at larger but finite Reynolds numbers (Hooper & Boyd 1987).It is fundamentally different from the previous one as its mechanism is not directly associated with the presence of the viscosity interface but, rather, of the wall.A comprehensive review of these instabilities for different flow configurations can be found in Govindarajan & Sahu (2014).
Plane Couette flow, which is linearly modally stable in the absence of stratification, was studied by Joseph (1964) in the presence of viscous heating, inducing temperature gradients and hence viscosity stratification.A linear inviscid instability was shown to develop if a liquid, rather than a gas, was considered.This observation was linked to the viscosity law, which decreases with temperature in liquids but increases in gas.While this result, as the authors themselves stressed, did not proceed from a rigorous stability analysis as the linearised energy equation was decoupled from hydrodynamic effects, this instability was recovered by numerical calculations in subsequent works (Sukanek et al. 1973;Yueh & Weng 1996).However, these studies did not consider density variations, which may arise when considering viscous heating in gases.Duck et al. (1994) carried out a stability analysis of plane Couette in a fully compressible framework.The authors mostly focused on acoustic instabilities appearing at supersonic Mach numbers, as also later studied by Malik et al. (2008) and Saikia et al. (2017).In addition to the acoustic modes, Hu & Zhong (1998) recovered the existence of a viscous mode similar to that found in the aforementioned incompressible, viscosity-stratified studies.

Recent developments in the hydrodynamics of supercritical fluids
Research on the hydrodynamics of fluids exhibiting non-ideal thermodynamic behaviour is actively progressing.A great deal of attention has recently been directed to understanding how the properties of these fluids affect turbulence, in particular turbulent heat transfer (Yoo 2013).Recent studies have investigated the statistics of turbulence in different shear flows by means of DNS, for example in channel (Nemati et al. 2015;Patel et al. 2016;Sciacovelli et al. 2017), pipe (Peeters et al. 2016;He et al. 2021), jet (Sharan & Bellan 2021) or flat-plate boundary layer flows (Kawai 2019;Sciacovelli et al. 2020).However, little is known about stability and transition to turbulence in these fluids (Robinet & Gloerfelt 2019).Gloerfelt et al. (2020) examined the linear stability of dense gas at large Mach numbers.Due to the large heat capacity of these fluids, very weak temperature gradients were observed and nearly incompressible velocity profiles were recovered.The authors showed the stabilisation of the viscous mode and the existence of radiating supersonic instabilities.From a different perspective, Ren et al. (2019a) carried out a linear stability analysis of supercritical fluids in plane Poiseuille flow.Having a lower heat capacity, significant viscous heating was present at reduced but non-negligible Mach numbers, generating temperature gradients in the base flow profile.The authors concluded that non-ideal effects may induce larger destabilisation of the flow in terms of growth rate magnitude and critical Reynolds number.In a subsequent study, Ren et al. (2019b) explored the linear stability of supercritical CO 2 in a flat-plate boundary layer flow.As viscous heating was increased, a second unstable mode, in addition to the classical Tollmien-Schlichting (TS) wave, was observed.This mode exhibits growth rates of more than one order of magnitude larger than the TS waves, which could imply new rapid modal routes of transition to turbulence in these fluids.The authors rigorously showed that this mode was not linked to the Mack's modes (Mack 1984) found in high-speed boundary layer.Bugeat et al. (2022) confirmed the inviscid nature of this instability and ruled out an acoustic origin.Recently, Ly & Ihme (2022) studied a binary compressible mixing layer at supercritical pressures and also found evidence of this instability, pointing out that its strength decreases as the reduced pressure is increased away from the critical point.But much remains to be understood about this instability as the driving parameters and the physical mechanism remain unclear.
Importantly, Ren et al. (2019b) observed that the additional mode only appears when the temperature profile of the base flow crosses the Widom line.The concept of Widom line is specific to supercritical fluids.It distinguishes the liquid-like from the gas-like region within the supercritical fluid domain.In each of these regions, fluid properties exhibit different behaviours (Simeoni et al. 2010).As such, the Widom line can be seen as the continuation of the coexistence line which separates the gas and liquid phases at sub-critical pressure, with the crucial difference that thermodynamic quantities smoothly vary across it (Banuti 2015;Banuti et al. 2017).These smooth variations can nonetheless exhibit remarkable behaviours.At constant pressure, density and dynamic viscosity, as functions of temperature, feature strong gradients near the Widom line, while the kinematic viscosity can reach a minimum; see the introduction of Ren et al. (2019a) for more details on these behaviours and their implication for hydrodynamics.Therefore, for a supercritical fluid operating at pressure and temperature near the Widom line, the presence of a temperature gradient in the flow leads to large density and viscosity variations; the flow is strongly stratified.

Objectives, method and assumptions
We aim to show that inviscid instabilities can be caused by the presence of a minimum of kinematic viscosity in strongly stratified shear flows, and that the scales associated with this minimum control the different properties of these instabilities.In particular, our objective is to provide evidence that the recently found unstable mode in supercritical fluids is related to the minimum of kinematic viscosity reached at the Widom line.We also aim to identify the driving physical mechanisms at play in these instabilities.
A differentially heated plane Couette flow will be considered.Three fluid models will first be examined, with different density and dynamic viscosity laws that strongly vary with temperature.Different types of stratification will then be observed in the flow; however, the three fluid models are designed to all feature a minimum of kinematic viscosity.In doing so, we aim to demonstrate the central role played by this minimum in the stability of the systems, regardless of the other property variations in the flow.Using piecewise linear base flow approximations, analytical results will be derived by solving the Rayleigh equation in the presence of strong density gradients, which governs the inviscid linear stability of these flows.A more realistic fluid model based on the van der Waals equation of state and diffusion laws at supercritical pressures will be used to ultimately discuss the instability in supercritical fluids.
The different hypotheses on the flow regime that we will consider in this work are summarised here.No assumption regarding the magnitude of viscosity and density variations will be made.Buoyancy will be ignored, but density variations will be retained in the inertial terms.Acoustics will not be taken into account in order to remove potential ambiguities in the physical interpretation of the results with aforementioned acoustic instabilities.The low-Mach approximation (Rehm & Baum 1978;Paolucci 1982) will be used.As a result, no viscous heating will be at play; temperature gradients will be generated in the flow by boundary conditions.Finally, only inviscid perturbations are considered.Note that this is not inconsistent with the presence of viscosity-stratification effects in the base flow which, because it is parallel, is not affected by inertia.The aforementioned instabilities induced by viscosity-stratification at low Reynolds number will therefore not be embedded in our analysis.However, it should be kept in mind that a competition may take place at finite Reynolds numbers, where the inviscid instability is damped by viscous effects.
The paper is organised as follows.The fluid and flow models, along with numerical procedures, are detailed in section 2. The condition of existence of an inviscid instability in stratified plane Couette flow is examined in section 3, leading to a criterion based on a minimum of kinematic viscosity.The base flows of the fluid models are presented in section 4. Analytical stability results, based on piecewise linear models of these base flows, are derived in section 5. Comparison with numerical calculations is provided in section 6.The generation of disturbance vorticity by different physical mechanisms is also examined, and an interpretation of the different instabilities is proposed.Section 7 eventually focuses on the stability of a supercritical fluid.A summary and a discussion on the application of these results to other fluids and flows conclude this paper (section 8).

Theoretical and numerical framework
2.1.Fluid models Four fluids will be considered throughout this paper, each of them being associated with a different equation of state and viscosity law.However, they all share the common property of assuming an extremum of kinematic viscosity  at a given temperature.Recalling that  = /, where  and  are the dynamic viscosity and the density, respectively, different ways to generate a minimum of  can be imagined.Three theoretical fluid models will be used to control and study a restricted number of parameters.A fourth more realistic model for supercritical fluids, based on the van der Waals equation of state, will also be considered.A summary of the different fluids is provided in figure 1 while a detailed description is given in the next subsections.

Fluid VB: Bump of dynamic viscosity with constant density
In this model, density is assumed constant, while the viscosity is chosen to locally exhibit a bump at a temperature  *  , using a Gaussian function (in this paper, all dimensional quantities are noted using the superscript '*').Note that, in this case, the non-dimensional kinematic viscosity is equal to the non-dimensional dynamic viscosity, both reading as reference temperature.The reference viscosities are  * ∞ and  * ∞ , which are the asymptotic values away from the bump.The parameter   controls the amplitude of the bump, and its sign determines whether the kinematic viscosity admits a minimum (  < 0) or a maximum (  > 0).The characteristic width of the bump is set through   , which is again made non-dimensional using  *  .Finally, the thermal conductivity  is assumed to be constant.

Fluid DB: Bump of density with constant dynamic viscosity
Inversely to fluid VB,  is kept constant in fluid DB while a bump is introduced in the density profile.This bump is chosen such that the resulting kinematic viscosity has the same expression as in (2.1).Hence, the density law simply reads and since, in this case,  = 1/,  is the same as in fluid VB.The conductivity  is again chosen to be constant.

Fluid HT: Hyperbolic tangent laws
In fluid HT, thermal conductivity is also kept constant, while dynamic viscosity and density are now both allowed to vary according to hyperbolic tangent laws.In order to generate an extremum of kinematic viscosity, a small shift is introduced between the two hyperbolic tangents, controlled by the non-dimensional parameter   .This choice is inspired by supercritical fluids and represents an attempt to mimic some of their features in the vicinity of the pseudo-boiling region.This will be discussed in more detail in section 2.1.4after the supercritical fluid laws are introduced.The non-dimensional governing laws for fluid HT are formally written as (2.4) The reference temperature  *  is here defined as the point of anti-symmetry of the density profile.The density at  = 1 and the viscosity at  = 1 +   are used as the reference scales.The parameter  controls the jump of density and dynamic viscosity while  , sets the temperature range over which this jump takes place.The fluid properties are shown in figure 1, where it is verified that the kinematic viscosity admits a minimum around  = 1.By analogy with fluid VB and DB, it is possible to estimate the amplitude   of this minimum, as well as the characteristic width   of the temperature range onto which it occurs.The following relations will be used in this paper:

,
(2.5) The derivation and verification of these expressions are detailed in appendix A.  Jossi et al. (1962) and Stiel & Thodos (1964) are used for the dynamic viscosity and the thermal conductivity, respectively.They provide analytical expressions for non-polar supercritical fluids based on theoretical scalings and experimental fittings.In supercritical fluids, these diffusion laws depend both on Ť and ρ.The density, dynamic and kinematic viscosity profiles are plotted in figure 1.Note that the reference temperature, here again noted  *  to maintain consistency with the previous fluids, is usually termed pseudo-boiling or pseudo-critical temperature in supercritical fluids (Banuti 2015).Density and dynamic viscosity are strongly correlated, and both exhibit strong gradients in the pseudo-boiling region.This motivated the choice of fluid HT, where density and viscosity are both defined using a hyperbolic tangent function, aiming at capturing these gradients while neglecting other variations away from them.The kinematic viscosity admits a minimum around the pseudo-critical temperature but is not localised, as opposed to the other fluids.The relatively simple model of fluid HT is found to decently reproduce this minimum as a result of the shift   introduced between the hyperbolic tangent laws, but differs away from the point where  remains strictly constant in fluid HT.
In analogy to the previous fluids, we would like to extract the characteristic scales   and   from the kinematic viscosity law.However, while  does have a minimum in fluid VdW, it is not clear that this minimum is localised over a finite, identified range   .Still, it can be observed, after calculation, that () admits two inflection points in the vicinity of the pseudo-boiling temperature -one below and one above this temperature.This can be used to define the scale   as the width between these two inflection points.From this, an amplitude   can be naturally defined.The procedure is thoroughly described in appendix B. Finally, note that the reduced pressure is the control parameter of the kinematic viscosity seen as a function of the temperature.In other words, () is different for each p, and, consequently, so are   and   .

Base flow
Linear stability analysis requires the knowledge of a base flow, defined as a steady solution of the non-linear Navier-Stokes equations.After recasting the non-linear Navier-Stokes equations given the physical assumptions associated with this flow, the equations are numerically solved.Plane Couette flow is comprised between two plates and is driven by the upper plate moving at speed  * 1 , which is used as the reference velocity scale.The streamwise and wall-normal directions are noted with  and , respectively.The flow is assumed to be parallel: the streamwise velocity  does not depend on , and the wall-normal and spanwise components of the velocity  and , respectively, are zero.The lower plate is fixed and, given the no-slip conditions, the non-dimensional streamwise velocity at the boundaries verifies (0) = 0 and (1) = 1.The distance ℎ * between the two plates is used as the reference length scale.The lower plate is kept at temperature  0 * , chosen as the reference temperature.We choose to consider the non-dimensional temperature gradient  between the two plates as an input parameter, which in turn sets the temperature of the upper plate.
The boundary conditions for the temperature are then  (0) = 1 and  (1) = 1 + .Under the assumption of a steady flow without pressure gradient -the flow is driven by the top wall, the non-dimensional Navier-Stokes equations reduce to a system of ordinary differential equations: where the superscript ′ denotes the wall-normal derivative and the overbars identify base-flow variables.The inertial terms are zero given the parallel flow assumption, and the problem does not depend on the Reynolds and Prandtl numbers.Besides, the temperature is decoupled from the velocity field.When  is constant, as it is supposed to be in fluids VB, DB and HT, the temperature profile is readily obtained as  () = 1 + .As for fluid VdW, equation (2.9) is solved using Newton's method by setting the initial guess as the aforementioned linear profile.Once  is obtained, the density profile is also known via the equation of state.The velocity profile is finally obtained by integration of equation (2.8) with the knowledge of the dynamic viscosity profile as a function of  and .Finally, note that we make the arbitrary choice to locate the extremum of kinematic viscosity at the centre line of the flow,  = 1/2.This is achieved by accordingly setting  *  / 0 * = 1 + /2 under the assumption that the temperature profile is linear -which is indeed the case for fluids VB, DB, and HT.

Rayleigh equation with density gradients
Assuming infinitely small, inviscid, two-dimensional perturbations, the linearised Navier-Stokes equations in the low-Mach approximation (Rehm & Baum 1978;Paolucci 1982) can be written: ) Perturbations of the form (, , ) = ℜ{ q()  ( −  ) } are now considered, with  = [, , ] being the state vector of the perturbations and ℜ{} the real part.These linearised equations can then be recast into the Rayleigh equation governing the linear dynamics of incompressible flows with density gradients (see also Fontane & Joly (2008)): (2.14) where  = / is the complex phase velocity.Note that the disturbance temperature does not appear in equation (2.14) since the linearised mass and momentum equations are decoupled from the energy equation (2.13).Temperature disturbances are deduced from the hydrodynamic disturbances, which can be calculated independently.Thermal effects are however at play in the velocity and density profiles of base flow, whose momentum equations (2.11) and (2.12), and ultimately the Rayleigh equation (2.14), depend on.A temporal framework is adopted: the wavenumber  is a real parameter while the frequency  is a complex number that is to be determined.The temporal growth rate is given by its imaginary part,   .A positive value corresponds to an inviscid instability.The (real) phase velocity   of the perturbation is simply   , the real part of .Equation (2.14) can be classically solved numerically as an eigenvalue problem.The boundary condition v = 0 is used at the wall.A pseudo-spectral method is employed to discretise the system and to obtain the derivative matrices (Orszag 1971).In order to avoid the singularity at the critical layer for neutral modes, a parabolic complex mapping is used, following Boyd (1985).This allows the growth rate to be computed even when it reaches small values, while a real mapping would produce spurious numerical oscillations.

Vorticity
In order to interpret some results, it can be useful to consider an alternative formulation of the problem in terms of the disturbance vorticity  = / − /.For a parallel base flow without pressure gradients,  is governed, in the physical space, by the linear equation where Ω is the vorticity of the base flow.The left-hand side represents the material derivative of  by the base flow.The right-hand side corresponds to vorticity sources, which may induce an instability.The term   is the production of vorticity responsible for shear flow instabilities.The second term,   , is the inertial baroclinic torque, which may generate vorticity when the density and pressure gradients are not aligned.In the absence of density gradients, this term is evidently zero.

Criterion of instability based on the kinematic viscosity profile
A necessary condition for an inviscid instability to exist was given by Rayleigh (1880) for constant-density flows.It requires the existence of an inflection point in the velocity profile of the base flow ū′′ = 0.In the presence of a density gradient, a generalisation of Rayleigh's theorem can be derived, often called the generalised inflection point (GIP) criterion in nonzero Mach number flow studies (Lees & Lin 1946;Mack 1984).Introducing the quantity Φ = − ρ ū′ , a necessary condition of inviscid instability is that Φ ′ = 0 somewhere in the base flow profile.The location where this condition is verified is termed GIP.Assuming that a GIP exists, an additional, more restrictive necessary condition of instability was given by Fjørtoft (1950).This criterion can be generalised to varying-density flows, stating that a region where with ū the velocity at the GIP, is required in the base flow profile in order to observe an inviscid instability.The proofs of these two results, stated in the case where density gradients are non-zero, straightforwardly follow those given in Schmid & Henningson (2001) for constant density flows by considering equation (2.14).For a monotonic velocity profile such as that of plane Couette flow, equation (3.1) must be verified everywhere (except at locations where Φ ′ ( ū − ū ) = 0).In this case, it is shown in appendix C that the generalised Fjørtoft's criterion (3.1) is equivalent to observing a maximum of |Φ| in the base flow profile.This extends the well-known interpretation of a maximum of absolute vorticity in constant-density flows.Indeed, noting that, under the parallel flow assumption, the vorticity of the base flow Ω is simply Ω = − ū′ , the quantity Φ can be interpreted as the density-weighted vorticity: For constant-density flows, the usual interpretation of the Fjørtoft's criterion is then recovered, since |Φ| = |Ω| in this case.However, in the presence of density variations, a maximum of vorticity is no longer a necessary condition of instability, and the existence of a maximum of |Φ| should instead be examined.Combining  = / and the streamwise momentum equation (2.8), it follows that ( νΦ) ′ = 0, (3.3) which, after distributing the wall-normal derivative, can be recast as The important result follows: in stratified plane Couette flow, the existence of a maximum of |Φ| is equivalent to the existence of a minimum of ν.Because of the generalised Fjørtoft's criterion, a minimum of kinematic viscosity in the base flow profile is then a necessary condition of inviscid instability.This motivated the choice of the fluid models considered in this paper (section 2.1), which all feature a minimum of  and, therefore, potentially exhibit an instability.Finally, note that equation (3.4) is specific to plane Couette flow.Different criteria may be expected for other shear flows, as discussed in appendix D.

Base flows of fluids VB, DB and HT
The base flows associated with the three fluid models VB, DB, and HT are presented in figure 2(a,b,c).The density and dynamic viscosity profiles have the same behaviour as those presented in section 2.1 -the constant temperature gradient of the base flow (section 2.2) providing a linear mapping from  to .Different velocity profiles are observed.In fluid VB, stronger gradients are present in the centre, where dynamic viscosity decreases.This is a result of the conservation of μ ū′ across the flow, yielding ū′ ∝ 1/ μ.While almost imperceptible in figure 2(a), the presence of these stronger gradients are clearly visible in figure 2(d), where the profile of |Φ|, as defined in equation (3.2), is shown.Indeed, in the case of fluid VB, density is constant and |Φ| reduces |Ω|.As for fluid DB, the velocity profile is linear (figure 2b) since viscosity is constant.Vorticity is therefore constant, but |Φ| still assumes a maximum at the centre as it now follows the density profile.Turning to fluid HT, the velocity profile features two regions of distinct gradient, which are, again, a consequence of the viscosity distribution (figure 2c).The resulting profile |Φ|, exhibiting a maximum in the central region as in the two previous fluids, is here a combination of the variations of density and vorticity.In summary, all fluids feature an excess of |Φ| in the central region.This is more generally understood because of the presence of a minimum of kinematic viscosity in each fluid as the integration of equation (3.4) leads to |Φ| ∝ 1/ ν.The normalised profiles of |Φ| all collapse (figure 2d) since identical parameters   and   are chosen for each fluid.Note that we will only consider   < 0 in order to generate a maximum of |Φ|, since no instability can occur otherwise according to the generalised Fjørtoft's criterion.

Stability models
5.1.Piecewise linear base flows Piecewise linear base flows have been extensively used to study a variety of stability problems with constant density as well as variable density -usually in the framework of the Taylor-Goldstein equation, under the Boussinesq approximation (Drazin & Howard 1966).While being simple approximations, useful analytical predictions can be derived from these models, especially predicting the linear stability of long waves (Gallaire 2015).We will here consider arbitrary large variations of density.
The base flows of fluids VB, DB, and HT are divided into three layers.The central layer, centred around  = 1/2, has a width  (figure 3).This approach follows that proposed by Rayleigh (1887) for bounded, constant-density flows.Velocity profiles are continuous at the interfaces between layers, but their gradient may not be; a vorticity jump may occur at the interface.This is the case of fluid VB, where the viscosity bump is modelled by a discontinuous jump in the central layer (figure 3a).This generates a stronger shear rate (i.e.,  vorticity) in the central layer.Note that this is precisely the configuration studied by Rayleigh (1887).In fluid DB, the shear rate is constant throughout the flow, but density exhibits a jump in the central layer (figure 3b).In fluid HT, density linearly varies in the central layer but remains constant in the two other ones (figure 3c).The same profile is used for the dynamic viscosity.All fluids feature an excess ΔΦ > 0 of |Φ| in the central layer in order to model the smooth profile of |Φ| that was observed in the previous section (figure 3d).

Expressions of 𝛿 and ΔΦ
The relations between the parameters  and ΔΦ of the sought stability model and the physical input parameters of our system   ,   , and  are now examined.The region of the smooth base flow profiles across which Φ varies is the same as that across which kinematic viscosity varies, as expressed in equation (3.4).The characteristic length of this region in the base flow profile is proportional to   , which is related to the fluid property, and inversely proportional to the temperature gradient  of the flow.The thickness of the layer  thus follows the proportionality relation: where the factor (1 + /2) results from the factor  *  / * 0 (see section 2.2) that appears when   is made dimensionless with  * 0 .A choice of a prefactor is ultimately required in order to assign a definitive value to  in equation (5.1), and will be specified for each fluid.
We define the quantity ΔΦ as the jump of |Φ| at the interface: ΔΦ = Δ(| − ρ ū′ |).The following convention is used: ΔΦ > 0 corresponds to configurations in which the magnitude of |Φ| is larger in the central layer than in the other layers.Given that ρ ū′ = μ ū′ / ν and that integrating the momentum equation (2.8) yields μ ū′ = , with  a positive constant, we can express ΔΦ as The characteristic value of ν being 1 +   in the central region and 1 elsewhere, the jump of 1/ ν at the interface reads (5.3) An excess of |Φ| in the central layer (ΔΦ > 0) is associated with a minimum of  (  < 0), consistent with section 3. The derivation of the different expressions of the constant  and the final expression of ΔΦ associated with each fluid model is detailed in appendix E.

Stability calculations in the long wave approximation 5.3.1. Derivation
The Rayleigh equation (2.14) is solved for the three piecewise linear base flows introduced in section 5.1.We will restrict our analysis to long waves ( ≪ 1).Following Rayleigh (1887) (see also Drazin & Reid (2004); Charru (2011)), equation (2.14) is first solved separately in each of the three layers of the piecewise linear base flows.As Φ is constant in each layer, the last term of the Rayleigh equation vanishes.Furthermore, assuming  ≪ 1 and writing v and  as a power series of , equation (2.14) reduces, at the order  0 , to (5.4) This equation can be solved in each layer.When density is constant across a layer, the solution is simply v =    +   , (5.5) where the index 'i' refers to the layer 1, 2 or 3 (see figure 3).In the central layer of fluid HT, where density varies, only the first derivative of v will be needed.This is because interfaces conditions, described hereafter, set the value of v2 using v1 and v3 (readily obtained from equation (5.5)).The following expression is immediately found: (5.6) At each interface between the layers, the kinematic and dynamic conditions (Charru 2011) read (5.8) Note that equation (5.7) reduces to Δ v = 0 as ū is continuous.Using equations (5.7) and (5.8) as well as the boundary conditions v = 0 at the walls lead to a linear system on the coefficients   and   .Equating the determinant to zero provides an expression of  2 .The system is unstable for  2 < 0 and stable for  2 > 0. In this section, derivations are carried out in the frame of reference moving with ū( = 1/2) = ū1/2 .

Results for fluid VB
As previously mentioned, the dispersion relation for fluid VB corresponds to that derived by Rayleigh (1887).We will use his result in the long wave regime.Using geometrical reasoning on the piecewise linear base flow of fluid VB (figure 3), Δ can be written as recalling that ΔΦ = ΔΩ for this fluid.Injecting this into Rayleigh's result yields (5.9) Note that ΔΦ was defined as ΔΦ = Δ(|Φ|) in section 5.2, and that ΔΦ > 0 corresponds to an excess of |Φ| in the central layer.The instability criterion is ΔΦ > ΔΦ  is deduced, with the instability threshold ΔΦ  : (5.10)

Results for fluid DB
The expression for  2 in fluid DB is solved by first remarking that Δ = , given that ū() =  in this case.It is also noticed that Δ = ΔΦ.Analytical calculations lead to the relation . (5.11) Because  < 1 and ΔΦ > −1 (since ρ > 0), the denominator is always strictly positive.The instability criterion is thus deduced from the sign of the numerator, leading to ΔΦ > ΔΦ  , with . (5.12)

Results for fluid HT
In fluid HT, the transformation of Δ into the parameters  and ΔΦ cannot be found.We will then use Δ = 2/(1 + ), which comes from the definition Δ = ρ(0) − ρ(1) introduced in figure 3 while using ρ * (0) as the reference scale.The stability model will then depend on  in addition to  and ΔΦ, but will become independent of  in some regimes of interest.The approximation Δ =  will be used, which is valid when the central layer or the viscosity jump are small.This will be shown to be of practical interest for the more realistic fluid VdW.Under these considerations, the following expression can be derived: (5.13) where the quantity is introduced to simplify equation (5.13).It will also be shown to be of practical interest for fluid VdW as ΔΦ 0 only depends on the jump of kinematic viscosity Δ(1/ ν) (appendix E).A criterion of instability is obtained by examining the sign of the expression under the square root in equation (5.13).Given that the factor 1 − ΔΦ 0 (1 − ) is always positive, the criterion of instability is given by the third factor, reading ΔΦ > ΔΦ  with . (5.15)

Comments and limiting case
The three different criteria all state that a certain excess of |Φ| is required in the central layer (ΔΦ > ΔΦ  ) to generate an instability.This can be seen as an improvement of the Fjørtoft's criterion, which states that an instability may occur only if ΔΦ > 0, but does not specify the magnitude of the excess of |Φ| that is sufficient to make the system unstable.It can be noticed that the derived instability thresholds always increase with .In the limit of small , fluid DB and fluid HT possess the same criterion of instability, equation (5.15) reducing to equation (5.12).Furthermore, in this regime, these equations differ from that of fluid VB (equation (5.10))only at order O ( 2 ).At order O (), all three fluids share the common criterion of instability: ΔΦ > .
(5.16)Moreover, the growth rate near the instability threshold can be calculated from the different expressions of  2 , using ΔΦ = O () ≪ 1.A general expression is obtained for the unstable modes of all fluids: =  √︁ (ΔΦ − ), (5.17)This shows the fundamental role the quantities ΔΦ and  play in modelling these instabilities, regardless of the types of stratifications.

Phase velocity
The equations on  2 obtained for each fluid also provide interesting results regarding the phase velocity of the unstable modes.For fluids VB and DB, if  2 < 0, then  is purely imaginary.Therefore, an unstable mode will have a phase velocity ū1/2 , which, by symmetry, is equal to 1/2 for these fluids.This does not hold for fluid HT, for which the phase velocity is shifted from ū1/2 by ΔΦ 0 (1 − ).By integrating the conservation of shear stress in the three layers, the following expression can be obtained for fluid HT: . (5.18) Note that the departure from 1/2 in equation (5.18) is of order , while the aforementioned additional shift is at most  2 for small ΔΦ 0 given that ΔΦ 0 ∼ |  / , | <  (appendix E).Equation (5.18) is therefore expected to be a good approximation of the phase velocity.

Numerical stability calculations for fluids VB, DB and HT
6.1.Growth rate and phase velocity Numerical stability calculations are carried out for fluids VB, DB, and HT using three different thicknesses of the central layer .Note that  is initially not an input parameter of  the problem: it is calculated from equation (5.1), which requires a prefactor.This prefactor, which is a priori different for each fluid, is set so as to yield the best agreement between the calculated and predicted stability diagram, which will be presented in the next subsection.The value 1.12 is used for fluids VB and DB.While not fully predictive -this value is not obtained by the model and requires one calculation point in order to be calibrated, it remains close to one: the model can provide order-of-magnitude predictions even without further knowledge.A prefactor equal to 1 is used for HT, requiring no external data.First, a constant value of   = −0.04,which sets the magnitude of ΔΦ for small  (E 9), is chosen for all fluids.Table 1 indicates the corresponding values of ΔΦ/ΔΦ  , providing a useful reference when comparing with figure 5, which will be presented later.As shown in figure 4(a,b,c), similar behaviours as well as close quantitative values of the growth rate are found for all fluids, despite the fundamentally different stratifications in each fluid.All fluids exhibit a long-wave instability: low wave numbers are always unstable, while a cutoff wave number   exists beyond which the system is stable.This observation justifies the restricted analysis to long waves to predict the stability of the system (section 5.3).As  increases,   decreases, as it would for a constant-density, unbounded shear layer for which   ∼ 1/ (Charru 2011).At constant   , the maximum growth rate increases as  is reduced.However, this does not hold for the growth rate of long waves, for which confinement can have a destabilising effect -increasing  being equivalent to approaching the walls closer to the central layer.This behaviour is not unexpected given that  2 depends on polynomials of  of degree larger than 1 (section 5.3).Differentiating equation (5.17) with respect to , a general estimate of the value   that yields the maximum growth rate is found to be   = ΔΦ/2.Using equation (E 9) for ΔΦ, we find that   ≃ 0.02, which is consistent with the observations.Note that a different behaviour is observed in a bounded, constant-density mixing layer, in which Healey (2009) found that confinement reduces the temporal growth rate of the instability.
The growth rate of long waves is presented in figure 5(a,b,c).The quantity   / is plotted, which corresponds to the slope of the growth rate at  = 0.For all fluids, the instability threshold is well predicted.The behaviour of the growth rate past the threshold

VB DB HT
Figure 5: Top: growth rate, normalised by , of long waves, obtained at  = 10 −2 .Bottom: associated phase velocities.Solid lines are the numerical results and dashed lines are the theoretical predictions given, for each fluid, in equations (5.9), (5.11) and (5.13).Black dashed lines correspond to the theoretical general growth rate of equation (5.17), obtained for  ≪ 1 near the instability threshold.Instability thresholds ΔΦ  are given in equations (5.10), (5.12) and (5.15).Theoretical phase velocities are given in the frame of reference associated with the lower wall using the values of ū1/2 discussed in section 5.3.6.For all figures, each colour is associated with a size  of the central layer and each column is associated with a fluid (VB, DB, and HT).
is also reasonably well captured, but piecewise linear models do not yield exact quantitative matches.The validity of the general approximation of the growth rate (5.17) is verified for  = 0.003.At  = 0.01, this approximation still captures the threshold well.Yet, significant departures from the full theoretical prediction are observed for fluids DB and HT; for these fluids, more terms are indeed neglected in the derivation leading to equation (5.17).At larger , noticeable differences appear for all fluids.Finally, the phase velocity of the unstable mode in each fluid is very well predicted by theoretical models (figure 5d,e,f).It is equal to 1/2 for fluid VB and DB, and does not depend on  and ΔΦ.The phase velocity markedly differs from 1/2 in fluid HT as the velocity of the base flow ū1/2 depends on both  and the viscosity ratio  (on which ΔΦ depends), as discussed section 5.3.6.

Stability diagrams
The stability diagram of each fluid is represented in the space (, ΔΦ) in figure 6. Neutral curves are calculated by detecting, for each ΔΦ, the value of  for which   / becomes close to zero -the threshold being set to 10 −4 .Calculations are carried out for long waves at  = 10 −2 .Theoretical predictions are generally in good agreement with the numerical results for all fluids.For fluid VB, good predictions are observed at low , but important discrepancies appear for  > 0.2.Such a mismatch is not observed in fluid DB, for which the neutral curve is still accurately predicted for  = 0.5, corresponding to a configuration in which the central layer occupies half of the flow.As for fluid HT, the prediction is also excellent.Note that, for this fluid, the range of ΔΦ that can be studied is limited by the range of   (appendix E), which is itself limited by  < 1 and the set value of   / , (appendix A).As predicted in equation (5.16), the neutral curve of all fluids collapses in the limit of  ≪ 1, following ΔΦ = .As  increases, higher-order terms in  destabilise fluid HT: the magnitude of ΔΦ required to produce an instability becomes smaller than .Conversely, higher-order terms stabilise fluids VB and DB.

The different sources of vorticity production
The generation of disturbance vorticity can be examined from the structure of the unstable modes in the physical space.The wall-normal velocity perturbations are made of a plane progressive wave along  (figure 7a).A peak is observed in the central region of the flow, consistent with the linear increase from zero at the wall predicted in the outer regions in equation (5.5).The associated pressure field is a plane progressive wave with a phase shift of /2 with respect to .Note that results here presented for fluid HT, but nearly identical fields are obtained for fluids VB and DB.
In fluid VB, vorticity production only results from the shear term   in equation (2.15), given that ρ′ = 0. Given its mathematical expression,   follows the same wave structure as , multiplied by the factor −Ω ′ .The vorticity profile |Ω| of the base flow, which is equal to |Φ| for this fluid, is only non-zero around the central layer.It features strong positive and negative gradients around the lower and upper interfaces, respectively (figure 2d).As a result, the structure of   contains two out-of-phase waves that are localised around each interface of the central layer, as seen in figure 8(a).The field of total vorticity production,   +   , has the same structure as   since   = 0.
A similar reasoning can be applied to fluid DB, in which only inertial baroclinic effects are at play given that the base flow vorticity is constant.The structure of the inertial baroclinic term   (equation (2.15)) is deduced from that of  and the profile of − ρ′ / ρ2 .Given that ρ = Φ for this fluid, the   field presented in figure 8(e) is found to be similar to that of   observed for fluid VB (figure 8a).As   = 0, it follows that the total vorticity production   +   also resembles that of fluid VB (figure 8g,h).
As for fluid HT, both   and   are active in the generation of disturbance vorticity, and their structure (figure 8c,f) is markedly different from that of two previous fluids.Both feature a peak around the central region.This is again a result of the profiles of −Ω and − ρ′ / ρ2 .Moreover,   and   exhibit a phase difference of .This is readily understood from the phase difference of /2 between  and , another phase shift of /2 being added to   as it contains  / (equation (2.15)).Despite being out-of-phase and having a similar structure, the sum of   and   is not zero.Instead,   +   is composed of two out-of-phase waves around each interface (figure 8i), similar to what was observed for fluids DB and VB.This behaviour can be linked to the profile of Φ with the following arguments.Around the central region, the -momentum equation (2.11) can be approximated to Φ ≃ − /.This results from / + ū/ being much smaller than Φ in this region, given that the phase velocity of  is here close to ū, and that the growth rate   is small (this is more evident in the spectral space, where this term is simply ( ū − )).Under this approximation, the linearised vorticity equation (2.15) can be recast as where the right-hand side corresponds to   +  ; this shows how the quantity Φ encapsulates both shear and inertial baroclinic effects.The profile of Φ ′ and the structure of  ultimately explain the structure of   +   in fluid HT (figure 8i).Note that the denominator ρ, attached to Φ ′ in equation ( 6.1), modulates the amplitude of the minimum and maximum of Φ ′ ; it can be seen, in figure 8(i), that larger peaks observed around the upper interface than those around the lower interface, given that ρ is smaller in the upper region.
Overall, the structure of the total vorticity production   +   is similar for each fluid, of the fluid stratification.The associated vorticity fields are finally displayed in figure 8(j,k,l).Their structure follows that of   +   , with an alteration resulting from advection effects (left-hand side of the vorticity equation (2.15)).The final picture is two vorticity waves travelling along each interface, with a phase difference of  minus a phase shift induced by advection.It can also be noted that these waves are generated around the critical layer   , defined as ū(  ) =   .In fluids VB and DB,   = 1/2 (section 5.3.6)which, by symmetry of the base flow, leads to   = 1/2.In fluid HT, the phase speed of the mode presented in figure 8 is   ≃ 0.45; we have verified that ū(1/2) ≃ 0.45.
The relative amplitudes of the terms plotted in figure 8 are hidden by the normalisation of each field.The maximum absolute value that each field reaches within the physical space (, ) is used in the normalisation procedure.The relative values are indicated in table 2, noting   ,   and  + the maxima reached by   ,   and (  +   ), respectively.Results are straightforward for fluids VB and DB for which   = 0 and   = 0, respectively.As for fluid HT, it is shown that the shear and baroclinic effects act with similar strengths.Furthermore, the total source of vorticity (  +   ) is about six times weaker than these effects as it proceeds from the interference of the two waves, cancelling out a large part of their amplitude.6.4.Interpretation These results can be interpreted within the wave-interaction theory as reviewed by Carpenter et al. (2011).In this framework, instabilities are seen as a result of vorticity waves developing along two interfaces that are located close enough so that each wave modifies the velocity show the disturbance streamwise velocity of a fluid parcel that is induced by a streamwise gradient of disturbance pressure.The magnitude of this velocity depends on ρ, which modifies the inertia of the fluid parcel whether it is located in the lower or upper region.
field of the other.The modified velocity field further deforms each interface, which yields additional vorticity production.This forms a positive feedback loop in which the two vorticity waves are amplified, constituting an instability.In order to achieve amplification, these waves have to be phase-locked, which was indeed observed for the three fluids in figure 8(j,k,l).The appearance of vorticity waves along each interface has been shown to be driven by the structure of Φ (equation (6.1)).The physical mechanisms giving rise to these vorticity waves are now further examined.
In fluid VB, the instability arises from an excess of vorticity in a localised layer.This is the essential ingredient of the Kelvin-Helmholtz instability, which is well-known and will not be further discussed.The mechanism of this instability has indeed been given following either kinematic (Batchelor 2000) or dynamic Charru (2011) arguments.Carpenter et al. (2011) also examined this instability from a wave-interaction perspective.
In fluid DB, only inertial baroclinic effects produce vorticity disturbance.This originates from misalignments between gradients of density and pressure (Soteriou & Ghoniem 1995).More precisely, given the flow assumptions, this misalignment can only be generated via the density stratification of the base flow and the gradient of pressure perturbations in the streamwise direction.Physically, two parcels of fluids at different wall-normal locations, having two different densities (that of the base flow), do not have the same streamwise acceleration when submitted to a streamwise perturbation of pressure (Reinaud et al. 1999).This induces a wall-normal gradient of streamwise velocity, i.e., vorticity.In fluid DB, the central region has an excess of density.The above mechanism of vorticity production, therefore, occurs between the lower and central layer, as well as between the upper and central layer.This results in the two vorticity waves that have been observed in figure 8(h).
As for fluid HT, the generation of the two vorticity waves is not as straightforward.A sketch of the mechanism is shown in figure 9(a,b,c).On the one hand, disturbance vorticity is generated following the inertial baroclinic mechanism described in the previous paragraph.But contrary to fluid DB, this occurs only between two regions of the flow, the lower and the upper ones, which have different densities (figure 2c) -the central layer playing the role of an interface between them.This idealised representation is illustrated in figure 9(b).
Therefore, only one interface is felt by the baroclinic perturbations, instead of two as in fluid DB.This results in only one vorticity wave, which is generated along the central region, as previously observed in figure 8(f).On the other hand, disturbance vorticity is also produced by a shear mechanism, which consists in the wall-normal transport of base flow vorticity by the perturbations.The base flow can also be divided into two regions of vorticity -smaller and larger magnitudes of the shear rate are indeed observed in the lower and upper part of the flow, respectively (figure 2c).The central layer again plays the role of an interface between these two regions, and only one interface is felt by the shear perturbations, as sketched in figure 9(a).As a result, only one wave is generated by the shear mechanism, as previously observed in figure 8(c).This is a consequence of the plane wave structure of , which, as it takes positive and negative values along , alternatively transports parcels of fluid that contain smaller and larger magnitude of Ω towards the central region.At this point, each mechanism generates one wave that is localised in the central layer.These two waves are out-of-phase, as discussed in the previous section.Moreover, because of the existing shift between the viscosity and density profiles in fluid HT (figure 2c), a small shift Δ ∼   / also exists between the interface at play in each mechanism (figure 9).As a result, the two central waves are in fact slightly shifted from each other, such that their superposition gives rise to two waves that are localised along each side of the central layer (figure 9c).Ultimately, it is these two phase-locked, interacting vorticity waves that produce the instability in fluid HT, following the interpretation of the wave-interaction theory.

Base flow
A typical base flow of fluid VdW is shown in figure 10(a).The general behaviour is similar to that observed for fluid HT in figure 2(c), the latter being, indeed, an attempt to model some important features of the former.Two regions of markedly different shear rate can be identified.The density and dynamic viscosity profiles are not as simple as in fluid HT: strong gradients are present in the central region, but properties also exhibit weaker variations away from it.The resulting profile of |Φ| 10(b) follows that of the kinematic viscosity profile, presented in section 2.1.4.As noted, the extremum of  (and therefore |Φ|) is seemingly not as localised as in the other fluid models.Modelling it with a piecewise view of the Φ profile, as proposed in figure 3(d), does not appear as an obvious choice.However, it was also noted that scales associated with the width and amplitude of the minimum of  could be introduced for fluid VdW (appendix B).By using them, we will show that the piecewise linear model of fluid HT can indeed provide a good representation of fluid VdW, allowing some stability features to be predicted, and thereby providing useful elements to interpret the instability.Note that, contrary to the previous fluids, the minimum of  is not exactly reached at  = 1/2 but is shifted upwards.This is linked to the procedure used to set the location of the minimum through the reference temperature at the wall and based on the assumption of a constant temperature gradient (see section 2.2).This assumption is not exact for fluid VdW as thermal conductivity varies, causing the observed shift.

Numerical stability calculations: comparison with the model of fluid HT
The stability diagram of fluid VdW can be calculated in the space (, ΔΦ) by varying both the reduced pressure p and the temperature gradient .Results are compared with the following theoretical predictions obtained for fluid HT in section 5.3.4.For small , equation (5.15) leads to the theoretical criterion of instability This equation is of interest for fluid VdW as it does not involve  (since ΔΦ 0 only depends on   , as shown in appendix E), which would have had to be defined for such a fluid.Note that equation (7.1) has an error of order O ( 2 ) (and not O (), as might have been expected), widening the validity of this approximation.A generally good agreement is obtained between the numerical calculations and the theoretical prediction (figure 11).This shows the robustness of the model based on fluid HT with respect to the property variations outside the central layer and to the upward shift of this layer.Furthermore, the definition of   and   by the inflection point of () (appendix B) is proved to be relevant, leading to a quantitative prediction of the neutral curve.Note that the prefactor used in the definition of  (equation (5.1)) is here equal to 1.The different terms of the vorticity equation (2.15) for the unstable mode in fluid VdW are now examined.The shear term   (figure 12a) and the inertial baroclinic term   (figure 12b) are each composed of a unique wave, which reaches a maximum in the central region (following the discussion of section 7.1, the central region is here defined as being centred around the minimum of kinematic viscosity).This is similar to the observations made for fluid HT in figures 8(c) and 8(f).But because non-zero gradients of μ and ρ persist away from the central region in fluid VdW (figure 10a), these terms are not as localised in the centre as in fluid HT.This is particularly noticeable for   , which extends further in the upper region, as ρ′ is non-zero and ρ is much smaller than in the lower region.The sum   +   contains two out-of-phase waves (figure 12c) given that as   and   are themselves out-of-phase and are slightly shifted from each other in the wall-normal direction; see the discussion for fluid HT in 6.4.However, because of the aforementioned asymmetrical structure of   , the upper wave is not localised around the upper interface of the central layer.This constitutes a significant difference with fluid HT (figure 8i).Nevertheless, the final picture is essentially identical: after advection is accounted for, the vorticity field contains two waves localised around each interface, with an additional phase shift leading to a phase difference smaller than  between them (figure 12d).The role played by advection into the localisation of   +   in the central region can be understood from the vorticity equation (2.15), which can be recast, in the spectral space, as Therefore, | | increases as the phase velocity of the disturbance approaches that of the base flow -which occurs in the central region, as will be shown in the next subsection.
Overall, the theoretical stability model developed for fluid HT predicts important features of the stability fluid VdW.This indicates that the main ingredients of the inviscid instability developing in supercritical fluids are indeed contained in fluid HT.These ingredients are the presence of strong, localised gradients of density and viscosity and the associated existence of a localised minimum of viscosity -whose characteristic scales control the instability through the more general parameters  and ΔΦ.The inviscid instability in supercritical fluids can ultimately be interpreted, like in fluid HT, as a result of the combination of shear and inertial baroclinic mechanisms which generate two interacting waves around the central layer (section 6.4).

Growth rate and phase velocity
Having shown the correspondence between the stability of fluids HT and VdW, further numerical results of fluid VdW can now be examined through the lens of the stability model.The growth rate of the unstable mode is shown in figure 13(a) for different supercritical pressures ( p > 1) and temperature gradients .The long-wave nature of the instability is recovered.At constant p, the maximum of the growth rate and the cutoff wavenumber increase with .This is consistent with the observations in figure 4(c) for fluid HT, since increasing  alone only affects , decreasing it.Note that the potential destabilisation of long waves by confinement effects (section 6.1) is not observed in the present range of parameters.When  is fixed, the growth rate is reduced as the supercritical pressure increases.The interpretation is not straightforward since p affects both   and   (appendix B), and therefore both ΔΦ and  in the stability model.Equation (5.17), which can reasonably be invoked given that the magnitudes of ΔΦ and  are here of the order of 10 −2 , can shed light on this behaviour.Both ΔΦ, which can be approximated by |   | in this regime, and , which is proportional to   , increase with p.Therefore, the reduction of the growth rate with p proceeds from a faster decrease of (ΔΦ − ) than the increase of .
Note that negative growth rates are obtained in figure 13(a) despite carrying out an inviscid stability analysis, in which only neutral modes are usually expected in the absence of an instability.This is a consequence of the use of a complex mapping for  (see section 2.3.1),which is used to remove the singularity of the critical layer for neutral modes.Whilst these negative growth rates do not carry any physical significance, this approach allows us to compute accurate cut-off wave numbers, presented later in this section.
The phase velocity of the unstable mode is displayed in figure 13(b).It is found to be reasonably constant for all wave numbers.The velocity of the base flow at the minimum of ν, ū , provides a good prediction of the phase velocity.This is consistent with the results of section 5.3.6,substituting ū1/2 by ū because of the aforementioned upward shift of the central layer.This is also consistent with defining the location of the central layer around the minimum of ν, as used in the previous subsection.A noticeable departure from ū can however be observed as  increases, i.e.,  decreases.This behaviour is qualitatively unexpected as the additional shift predicted by the theoretical model should decrease with .
The cutoff wavenumber   can be plotted from the growth rate calculations by detecting the value of  ≠ 0 where   = 0. Results at two pressures are presented in figure 14.For each pressure, the cutoff wavenumber is calculated for a range of temperature gradients  in order to vary .The scaling   ∼  −1 is found.While this scaling was not derived from the stability models -the long wave approximation does not allow   to be examined, this result is analogous to that of constant-density mixing layers (Charru 2011).Note that the scaling is not valid as  increases towards the stability threshold.Indeed, as  is finite, it would predict finite values of   near the threshold, when   tends to zero.

Summary
The stability of strongly density-and viscosity-stratified plane Couette flow was examined.It was shown that a minimum of kinematic viscosity in the base flow profile produces a GIP, which furthermore satisfies the generalised Fjørtoft's criterion.We first consider three fluid models (termed VB, DB, and HT) which were designed to all exhibit a minimum of , while having different stratifications of density and dynamic viscosity.Modelling their base flow via piecewise linear profiles, the Rayleigh equation that accounts for strong density gradients was solved for each of them in the long-wave approximation.Expressions of the growth rate and phase velocity were derived, as well as a criterion of instability for each fluid.All these criteria express that a sufficient excess of |Φ|, the density-weighted vorticity of the base flow, shall be localised in a layer of thickness .This improves, for the specific flow studied in this paper, the generalised Fjørtoft's criterion -which does not provide such an instability threshold.A general criterion was obtained for all fluids in the limit of small , regardless of the type of stratification.
Theoretical predictions were compared to numerical stability calculations.A qualitatively good agreement was found for the growth rate, while an excellent match was observed for the phase velocity.Stability diagrams are also generally well predicted by the models, and the collapse of the neutral curves for small  was indeed observed in the calculations.The physical mechanisms responsible for the instability in each fluid were then examined.For all fluids, two phase-locked vorticity waves travelling along each interface of the central layer were found.The growth of the instability is ultimately interpreted as a result of the interaction between these two waves.These waves are generated by either shear effects in fluid VB or inertial baroclinic effects in fluid DB.In fluid HT, which features strong stratifications of both density and viscosity, the two waves were shown to result from a combination of shear and inertial baroclinic mechanisms.
The stability of a fluid governed by the van der Waals equation of state at supercritical pressure was finally examined.The stability model was found to quantitatively predict the neutral curve of this more realistic fluid.The concurrent action of shear and inertial baroclinic effects in the vorticity production was shown to produce two vorticity waves around the central layer.This suggests that the physical interpretation proposed for fluid HT also provides the essential mechanism of the instability in supercritical fluids.Ultimately, our study provides evidence that the minimum of kinematic viscosity, reached at the Widom line, is the leading cause of this instability.

Supercritical fluids and beyond
We finally address the link between the present paper and the unstable mode recently found in supercritical fluids by Ren et al. (2019b), before discussing how this study may be relevant to other types of fluid and flow.It has been shown that, in plane Couette flow, a heated supercritical fluid features a similar instability to that developing in simpler fluids, whose only common property was to assume a minimum of kinematic viscosity.This indicates that other non-ideal thermodynamic effects, such as non-unity compressibility factor , large values of   exhibited near the Widom line or non-monotonic, large variations of the speed of sound, are unlikely to play a decisive role in the instability mechanism.In the original study of Ren et al. (2019b), the existence of the additional unstable mode was reported for a flat-plate boundary layer flow, considering supercritical CO 2 at non-negligible Mach numbers.We aimed at simplifying their configuration by considering a simpler flow (plane Couette), a simpler supercritical fluid model (VdW equation of state and simple analytical diffusion laws) and by neglecting acoustic phenomena.In this more canonical framework, the instability could be recovered and further analysed.The role played by the minimum of kinematic viscosity remains to be clearly demonstrated in flat-plate boundary layer flows as studied by Ren et al. (2019b), since the theoretical development presented in section 3 cannot directly apply to this flow.However, the authors observed the additional unstable mode only when the temperature profile crosses the Widom line, about which the presence of a minimum of  is a common property (see figure 15a).This indicates strong links with our results.Bugeat et al. (2022) furthermore pointed out how, in some cases, inertial effects can be neglected around the pseudo-boiling region, resulting in the correct prediction of a GIP in the boundary layer.In this case, the momentum equation reduces to that of Couette flow, making the aforementioned theoretical result valid for a flat-plate boundary layer flow.But a rigorous understanding of the conditions of the existence of a GIP is so far missing in this flow, and further efforts could be directed towards the search for a sufficient condition which would factor in inertial effects.The link with an excess layer of Φ should also be analysed in this case, as the criterion of instability in equation (7.1), involving ΔΦ and , may be altered for a different flow configuration.
Other types of flow may exhibit a localised minimum of viscosity, hence, potentially, a similar instability.In the present study, the scalar quantity that controls the spatial distribution of  is the temperature.Other quantities could play an analogous role.The case of a shear flow made of two miscible fluids, mentioned in the introduction, is an interesting example.The mole fraction is constant in each fluid but varies through the interface, which can be assumed to have a small, finite thickness.Ern et al. (2003) studied the instability developing in this system at low Reynolds numbers.It can be noted that certain fluid mixtures may exhibit a minimum of kinematic viscosity for intermediate mole fraction (between 0 and 1), as shown in figure 15(b).For these fluids, a minimum of  would then exist within the diffused interface, as it does for supercritical fluids in the region of the Widom line.An inviscid instability can be expected in this case, which may compete with other instability at low Reynolds numbers if viscous damping does not stabilise it too quickly.As pointed out by Ern et al. (2003), a stability analysis of such a system is valid as long as the time scale associated with the diffusion of this interface is large compared to that of the instability.Similarly, in a single (not necessarily supercritical) fluid flow, heating or cooling a small layer of liquid or gas, respectively, would produce a localised minimum of .An instability could then appear, provided that its associated time scale is small compared to that of thermal diffusion of the cooled or heated layer.Note that these time scale considerations were avoided in the present study, as the temperature field strictly was a steady solution of the Navier-Stokes equations.A comparison between the viscous time scale of the perturbations and that of the instability could however evaluate the potential of prediction of the inviscid analysis at finite Reynolds numbers.
Finally, the case of non-Newtonian fluids can be mentioned as dynamic viscosity, and hence kinematic viscosity, may vary as a function of the stress profile.In a canonical plane Couette, no variation of stress would be observed, and the flow would be linearly stable.Introducing temperature gradients through viscous heating can alter the stability of the flow (Eldabe et al. 2007).In an isothermal configuration, Nouar & Frigaard (2009) added a streamwise pressure gradient, leading to a plane Couette-Poiseuille flow.The presence of a maximum (minimum) of stress for a shear-thinning (shear-thickening) fluids may then generate a minimum of kinematic viscosity.However, the criterion of instability derived in the present paper would not apply because of the presence of a pressure gradient, and an isothermal parallel flow of non-Newtonian fluid may therefore not experience this instability.This can be interpreted as the need for the kinematic viscosity law of the fluid to contain a temperature scale, defined as Δ  = /(/), that is smaller than the temperature scale Δ of the flow.

Appendix E. Expression of 𝐾 and ΔΦ for different fluids
Integrating ū′ = / μ between 0 and 1 and using the boundary conditions on the streamwise velocity provides an expression of : (E 1)  only depends on the profile of dynamic viscosity.E.1.Fluid VB Since ρ is constant, the profile of dynamic viscosity is the same as that of kinematic viscosity.The integral in equation (E 1) can be approximated supposing that μ = 1 +   on a layer  while μ = 1 elsewhere.This leads to .
(E 2) The jump of |Φ| given by equation (5.2) can finally be expressed, for fluid VB, as .
( (E 6) We will approximate this expression in the limit of small , which will prove useful when applied to fluid VdW.Under this assumption, we simply have  = 1 − , and ΔΦ reads Note the expression of the quantity ΔΦ 0 , introduced in equation (5.14), is simply which does not depend on .Using ΔΦ 0 is then found of practical interest for fluid VdW as the parameter  does not need to be introduced and defined for this fluid.
E.4.Limiting case In the limit of  ≪ 1,   ≪ 1 and  ≪ 1, all fluids exhibit the same expression: ΔΦ ∼ −  .(E 9) This limit is of interest since it corresponds to parameters near the neutral curve when  ≪ 1, since ΔΦ = O () in this region.

Figure 1 :
Figure 1: The four fluid models considered in this paper.The definition of  ∞ , used to normalise , is given in appendix B in the particular case of fluid VdW.

Figure 2 :
Figure 2: Top: base flow profiles of fluids VB, DB and HT, for   = 10 −2 ,   = −10 −1 and  = 0.2.The inset in (a) is a close-up of the velocity profile in the central region, with a comparison to the linear function  () =  shown in black dashed line.Bottom: resulting profile of |Φ| for each fluid.

Figure 4 :
Figure 4: Growth rate, as a function of , of the fluid models VB, DB, and HT obtained for   = −0.04 and different thicknesses of the central layer .

Figure 8 :
Figure 8: Terms of the vorticity equation (2.15), plotted in the physical space (, ) at an arbitrary time , for the unstable modes presented in figure 7. The shear term, the baroclinic term, their sum, and the resulting vorticity are shown from top to bottom, in that order.Each column corresponds to one fluid.In each figure, fields are normalised by their maximum absolute value, and colour bars are the same as in figure 7. The central layer of thickness  = 10 −2 is shown in (a) in red dashed line; its location is identical in each figure, albeit not reproduced in order to ease visualisation.The location  = 1/2 is indicated in black solid line in (c) and (f), revealing an offset of the wave below and above this line, respectively.

Figure 9 :
Figure 9: Sketch illustrating the mechanism generating the two vorticity waves in fluid HT.Red and blue vortices indicate positive and negative values of the vorticity disturbance , respectively.These disturbances are generated by shear (a) and inertial baroclinic (b) mechanisms.Solid black lines represent the vorticity (a) and density (b) profiles of the base flow, Ω and ρ, respectively.Note that Ω < 0 and that a shift Δ exists between the interfaces of Ω and ρ.Vertical green arrows in (a) represent the transport of Ω by the disturbance  in and out of a control volume centred around the interface.Horizontal orange arrows in (b)show the disturbance streamwise velocity of a fluid parcel that is induced by a streamwise gradient of disturbance pressure.The magnitude of this velocity depends on ρ, which modifies the inertia of the fluid parcel whether it is located in the lower or upper region.

Figure 12 :
Figure 12: Terms of the vorticity equation (2.15) for the unstable mode in fluid VdW at  = 10 −2 , p = 1.06 and  = 0.5.Normalisation, colour bars, and annotations are identical to those detailed in figure 8.Note that, in fluid VdW, the central region is defined as being centred around the minimum of kinematic viscosity.

Figure 13 :
Figure 13: Growth rate (a) and phase velocity (b) of the unstable mode in fluid VdW.Results in solid lines are obtained for different values of  at p = 1.03.For  = 1, results at p = 1.06 and p = 1.09 are also presented in dashed line.The circles, the square and the triangle in (b) indicate the value of ū (see text) at p = 1.03, p = 1.06 and p = 1.09, respectively (these values are not a function of  and are added to this plot for comparison purposes).

Figure 14 :
Figure 14: Cutoff wavenumber obtained for fluid VdW at two reduced pressures.

Figure 15 :
Figure 15: (a) Kinematic viscosity of different real fluids at supercritical pressure p = 1.1 as a function of temperature.The subscript 'pb' refers to the pseudo-boiling point.(b) Kinematic viscosity of a mixture of methanol and toluene as a function of the mole fraction for different temperatures, from McAllister (1960).
Declaration of interests.The authors report no conflict of interest.

Figure 17 :
Figure 17: (a) Kinematic viscosity around the pseudo-boiling point of fluid VdW at supercritical pressures p = 1.03, p = 1.06 and p = 1.09.(b) Calculation of the scales   and   associated with the minimum of  in fluid VdW, for p ∈ [1.02; 1.10].
Fluid DB Dynamic viscosity is constant, which immediately gives  = 1.The jump of Φ for fluid DB is thereforeΔΦ () = −   1 +   .(E4) E.3.Fluid HT The profile of dynamic viscosity is assumed piecewise linear and equal to the profile of density of fluid HT in figure 3.In the lower wall region, μ = 1, In the upper wall region, μ = (1 − )/(1 + ), after renormalising  * with  * 0 in equation (2.4).The central region linearly matches these two regions.It follows of variable  =  − (1 + )/2 was performed for the second integral.After integration

Table 1 :
Values of ΔΦ/ΔΦ  for the different fluids and values of  presented in figure 4.