1 Introduction

Interfaces and interfacial mixing govern a broad range of processes in nature and technology, and their non-equilibrium dynamics couple micro to macro scales [1, 2]. Examples include: thermonuclear flashes in supernovae, laser ablated plasmas in inertial confinement fusion, chemically reactive flows, the detonation of energetic materials, the purification of water, electro-catalysis, and nano-fabrication [3,4,5,6,7,8,9,10,11,12]. In these realistic environments, the fields change sharply and rapidly, the phases of matter are well pronounced, accelerations are strong, energy releases are high, and relaxations processes are weak [1,2,3,4,5,6,7,8,9,10,11,12]. Non-equilibrium interfacial dynamics provides nearly insurmountable obstacles to study at kinetic and at continuous scales. Its understanding is of critical importance for science, mathematics, and engineering [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. The theory of interfaces and interfacial mixing must rigorously solve a singular boundary value problem and an ill-posed initial value problem at a freely evolving discontinuity. This mathematical task is even more challenging than the Millennium problem on the Navier–Stokes equation [1, 4, 11, 16].

Our work focuses on the long-standing problem problem of stability of a phase boundary—a fluid interface with fluxes of heat and mass across it [16]. The phase boundary is broadly defined. It can be an interface between two different matters (fluids, plasmas, solid materials) or an interface between the same matters with distinct thermodynamic properties. The matter can also experience a phase transition, undergo a change in chemical composition, and be out of thermodynamic equilibrium [3,4,5,6,7,8,9,10,11,12,13,14]. While interfaces appear obvious at a first glance, they are a challenge to rigorously determine. For instance, one may employ a front with zero mass flux to describe immiscible fluids, and an interface across which mass can be transported to describe miscible fluids. One may further presume that the fronts are ‘thin’ and the interfaces are ‘thick’ [12, 13]. Realistic environments are more complex, as found in molecular dynamics simulations of energetic materials and in high resolution experiments of solvents [7,8,9, 18]. Often, a few nanometers thick interface can have strong and macroscopically significant fluxes of heat and mass across it; moreover, at microscopic scales, these transports can be essentially non-diffusive [7,8,9,10, 18,19,20].

At an interface the properties of matter experience dramatic changes at small scales [14, 17]. These changes produce microscopic interfacial transports, which, in turn, define macroscopic fields in the bulk [3, 7,8,9, 11, 17]. In order to treat rigorously the multi-phase dynamics, to provide reliable benchmarks for diagnostics, and to achieve a high predictive capability in a broad range of conditions, the theory must balance the fluxes of the conserved quantities at the freely evolving interface and must solve the boundary value and the initial value problems [4, 16, 17]. For fronts in Rayleigh–Taylor and Richtmyer–Meshkov instabilities, accurate theories were built in the last decades. Particularly, group theory approach grasped the order in Rayleigh–Taylor mixing and explained the observations [11, 19, 21,22,23,24,25,26,27]. For interfaces with interfacial mass flux, a rigorous theory of the interface stability was developed recently for ideal fluids [17, 28,29,30]. This theory discovered the inertial mechanism of interface stabilization, the instability of the conservative dynamics of accelerated interface, and the chemistry-induced instabilities [17, 28,29,30]. It also resolved the prospect of Landau 1944 and found that classical Landau’s solution for Landau–Darrieus instability is a perfect mathematical match [16, 17]. Resolutions of two other prospects of the 1962 Nobel Laureate Landau were recognized with Nobel prizes in 2003 (theory of superconductors) and 1982 (theory of phase transitions) [31].

Theory [17, 28,29,30] studied the dynamics of the interface with interfacial mass flux for ideal incompressible fluids free from thermal heat flux. It also analyzed the effect of surface tension understood as tension at the phase boundary [31, 73]. Realistic processes are often accompanied by fluxes of thermal heat across the interface and by non-ideal thermodynamics in the bulk, including thermal conductivity caused by inelastic interactions of particles at atomic scales [14, 15, 17]. The influence of the thermal heat flux and the microscopic thermodynamics on the interface stability and the flow fields’ structure requires systematic investigations. These studies are critically important to solve the corner-stone problem of modern science and mathematics and to better understand a broad range of processes in nature and technology, to which interfaces and interfacial mixing are relevant [1,2,3,4,5,6,7,8,9,10,11,12, 14, 17, 31].

In this work, we consider the inertial and accelerated dynamics of the interface with fluxes of heat and mass across it for non-ideal thermally conducting nearly incompressible fluids. We develop the rigorous theoretical framework applicable in a broad range of conditions [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. We resolve the fundamental challenges not addressed before [12,13,14]. These include: the self-consistent boundary conditions for the thermal heat flux; the perturbation waves’ structure; the dependence of the waves coupling on the system parameters. We report key discoveries regarding three regimes—advection, diffusion, and low Mach—exhibiting fluid instabilities which were never earlier discussed and which are defined by the interplay of the thermal heat flux, thermal conductivity, destabilizing acceleration and inertial stabilization. We find that the interface stability is achieved primarily through the macroscopic inertial stabilization balancing the destabilizing acceleration. The microscopic thermodynamics and the thermal conductivity lead to creation of vortical fields in the bulk. The strength of these vortical fields is set by the thermal heat flux. Our analysis defines the interface as the place where balances are achieved. Our theory directly couples the macro and micro scales in the non-equilibrium dynamics for a broad range of conditions, including matter at the extremes [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17, 19].

Potential applications of our theory span across the scales and disciplines, since interfaces and their stability impact nearly every area of science and technology. In supernovae blasts, unstable interfaces and interfacial mixing of materials of the progenitor star create conditions for nucleosynthesis of heavy and intermediate mass chemical elements (in addition to light mass elements synthesized in the star before its explosion) [3, 4, 34, 35]. The formation of Sun’s spots and solar flares is strongly influenced by downdrafts—finger type structures pushing the matter from the solar surface into the convection zone [36]. In the inertial confinement fusion, the fluid phases—the regions of the hot and cold neutral plasmas—are formed by sharply and rapidly changing fields, and the shock induced interfacial mixing of these regions can preempt the formation of hot spot [5, 37]. Interfaces and their stability are inherent in the light—matter interaction, in the water exploding plasmas, and in the plasma discharges formed in and interfacing with liquids [36, 38,39,40]. Dynamics of interfaces separating the flow heterogeneities is critical for understanding of realistic turbulent processes, including compressible turbulence, turbulent convection, turbulent combustion, and turbulent mixing [6, 12, 13, 41,42,43,44].

Interfaces with interfacial fluxes of heat and mass are essential in the processes of the atomization of liquid jets and the fluid vaporization [45,46,47], in the melting and evaporation of materials under high pressures and high strain rates, and in chemically reactive flows [7, 9, 28, 48]. A good grip on the interface dynamics is needed for efficient purification of water and for enabling the transportation security of liquefied natural gas [20, 49]. The insights on the coupling of microscopic transports at the interface and macroscopic flow fields in the bulk is critical in nano—science and nano—technology, including the detonation of energetic materials, the understanding of fluid motions in fuel cells and in micro-channels, and the in-depth comprehensions of the electro-catalysis and the properties of matter undergoing phase transitions [18,19,20, 50,51,52,53,54,55].

In these vastly different physical regimes, interfaces and interfacial mixing are observed to have similar features of their non-equilibrium dynamics [1, 2]. They can be investigated within the theoretical framework of the classical problem of stability of the phase boundary when there are heat and mass fluxes across the interface. This problem is the focus of our present study.

The novelty of our work is in the development of general rigorous theory of interface dynamics in realistic fluids applicable in a broad range of conditions. We explore the new class of fluid instabilities in the three regimes with properties not discussed before, including the growth-rates and the flow fields’ structures. A special contribution of our work is the control parameter describing transitions between the regimes of advection, diffusion and weak compressibility by varying the initial conditions. We discover that interface is a place where balances are achieved through linking micro to macro scales.

The paper has the following structure: After Introduction in Sect. 1, we present Theoretical Problem in Sect. 2, including governing equations (2.1) and linearized dynamics (2.2). We provide Results in Sect. 3: the new methodology (3.1); the joint properties of solutions (3.2); dynamics for regimes of advection (3.3), diffusion (3.4) and low Mach (3.5); the link of interface dynamics and flow fields (3.6); the range of applicability (3.7); the characteristic scales (3.8); the comparison with experiments (3.9); the properties of thermal heat flux (3.10); the design of experiments (3.11). Discussion and Conclusion are in Sects. 45. Acknowledgement, Data Availability, Author’s contribution, Conflict of interest, References, Figures and Tables are given in last sections.

2 Theoretical problem

2.1 Governing equations

In an inertial frame of reference, the equations for the conservation of mass, momentum and energy, and the heat flux equation are:

$$ \frac{\partial {\uprho } }{{\partial t}} + \frac{{\partial {\uprho } v_{i} }}{{\partial x_{i} }} = 0,\quad \frac{{\partial {\uprho } v_{i} }}{\partial t} + \frac{{\partial {\uprho } v_{i} v_{j} }}{{\partial x_{j} }} + \frac{\partial P}{{\partial x_{i} }} = 0,\quad \frac{\partial E}{{\partial t}} + \frac{{\partial \left( {E + P} \right)v_{i} }}{{\partial x_{i} }} + \frac{{\partial Q_{i} }}{{\partial x_{i} }} = 0,\quad Q_{i} + \frac{{\partial \left( {\chi e} \right)}}{{\partial x_{i} }} = 0 $$
(1a)

with spatial coordinates \({\mathbf{x}} = x_{i} = \left( {x_{1} ,x_{2} ,x_{3} } \right) = \left( {x,y,z} \right)\), time \(t\), and thermal conductivity \(\chi\), and with the fields of density \({\uprho }\), velocity \({\mathbf{v}} = v_{i}\), pressure \(P\), energy density \(E = \rho (e + {v_i}^2/2)\), specific internal energy \(e\), and heat flux \({\mathbf{Q}} = Q_{i}\). The closure equation of state relates internal energy and pressure. We presume it in the form \(P = s\;{\uprho } \;e\) with constant \(s\) [14, 17]. The inertial frame of reference is referred to the frame of reference moving with constant velocity \({\mathbf{V}}_{0}\); for definiteness \({\mathbf{V}}_{0} = \left( {0,0,V_{0} } \right)\) [17]. Equations Eq. (1a) describe non-ideal thermally conducting fluids [14]. They are reduced to the Euler equations for ideal fluids at \({\mathbf{Q}} = 0\) and \(\chi = 0\) [14], for which the interface dynamics is studied in [17, 28,29,30, 32, 33].

In engineering applications, one may apply the thermal resistance \(r\) in the equation for energy conservation and the heat flux equation, which is associated with thermal conductivity as \(r\,\sim\,\chi^{ - 1}\) [14]. In this work, we employ the thermal conductivity, as per the standards of theoretical physics [14]. For the constant thermal conductivity \(\chi = const\) and with no motion dynamics \(v_{i} = 0\), the equation for energy conservation and the heat flux equation in the system Eq. (1) are reduced to the Fourier’s law of heat conduction as \(\partial e/\partial t + \chi \partial^{2} e/\partial x_{i}^{2} = 0\) [14].

The derivation of the governing equations Eq. (1) and the explanation of non-ideal character of dynamics of thermally conducting fluids are given in classical work [14]. The critical aspect of our approach to the problem of the interface dynamics is that the governing equations Eq. (1) are derived in an inertial frame of reference moving with constant velocity \({\mathbf{V}}_{0}\) relative to laboratory frame of reference, rather than in a frame of reference moving with the interface velocity \({\tilde{\mathbf{V}}}\). This permits us to stay free from a postulate of constancy of the interface velocity (even for the stable dynamics) and to identify the inertial stabilization mechanism caused by the motion of the interface as whole. See for details below and works [17, 28,29,30].

We consider multi-phase dynamics, and introduce a continuous scalar function \({\uptheta} \left( {x,y,z,t} \right)\) describing the fluid interface, such that \({\uptheta} = 0\) at the interface and \({\uptheta} > 0\;\left( { < 0} \right)\) in the heavy (light) fluid. We present \(\left( {{\uprho } ,{\mathbf{v}},P,E,e,{\mathbf{Q}},\chi ,s} \right) = \left( {{\uprho } ,{\mathbf{v}},P,E,e,{\mathbf{Q}},\chi ,s} \right)_{h} H\left( {\uptheta} \right) + \left( {{\uprho } ,{\mathbf{v}},P,E,e,{\mathbf{Q}},\chi ,s} \right)_{l} H\left( { - {\uptheta} } \right)\) in the bulk by using the Heaviside step-function \(H\left( {\uptheta} \right)\), with subscript \(h\left( l \right)\) standing for heavy (light). At the interface, the fluxes of mass, momentum, and energy are balanced [16, 17]. We also require that thermal heat flux is normal to the interface on each side of the interface.

$$ \begin{array}{*{20}l} {\left[ {{\tilde{\mathbf{j}}} \cdot {\mathbf{n}}} \right] = 0,\quad \left[ {\left( {P + \frac{{\left( {{\tilde{\mathbf{j}}} \cdot {\mathbf{n}}} \right)^{2} }}{{\uprho } }} \right)\,{\mathbf{n}}} \right] = 0,\quad \left[ {\frac{{\left( {{\tilde{\mathbf{j}}} \cdot {\mathbf{n}}} \right)\left( {{\tilde{\mathbf{j}}} \cdot {{\varvec{\uptau}}}} \right)}}{{\uprho } }\,{{\varvec{\uptau}}}} \right] = 0,\quad \left[ {\left( {{\tilde{\mathbf{j}}} \cdot {\mathbf{n}}} \right)\;\left( {W + \frac{{{\tilde{\mathbf{j}}}^{2} }}{{2{\uprho }^{2} }}} \right) + {\mathbf{Q}} \cdot {\mathbf{n}}} \right] = 0} \hfill \\ {\left. {\left( {{\mathbf{Q}} \cdot {{\varvec{\uptau}}}} \right)} \right|_{{{\uptheta} = 0^{ + } }} = 0,\quad \left. {\left( {{\mathbf{Q}} \cdot {{\varvec{\uptau}}}} \right)} \right|_{{{\uptheta} = 0^{ - } }} = 0} \hfill \\ \end{array} $$
(1b)

Here \(\left[ {...} \right] = 0\) denotes the jump of quantities at the interface. The unit normal and tangential vectors of the interface are \({\mathbf{n}},\;{{\varvec{\uptau}}}\) with \({\mathbf{n}} = \nabla {\uptheta} /\left| {\nabla {\uptheta} } \right|,\;\left( {{\mathbf{n}} \cdot {{\varvec{\uptau}}}} \right) = 0\). Mass flux across the interface is \({\tilde{\mathbf{j}}} = {\uprho } \left( {{\mathbf{v}} + {\mathbf{n}}\dot{{\uptheta} }/\left| {\nabla {\uptheta} } \right|} \right)\). The specific enthalpy is \(W = e + P/{\uprho }\) including the enthalpy of formation [16]. Emphasize that the fluid velocity \({\mathbf{v}}\) is shear-free at the interface, as is clearly seen from the third equation in the system Eq. (1b). The boundary conditions at the outside boundaries are

$$ \left. {{\mathbf{v}}_{h} } \right|_{z \to - \infty } = {\mathbf{V}}_{h} = \left( {0,0,V_{h} } \right),\quad \left. {{\mathbf{v}}_{l} } \right|_{z \to + \infty } = {\mathbf{V}}_{l} = \left( {0,0,V_{l} } \right) $$
(1c)

The interface velocity is \({\tilde{\mathbf{V}}}\). For steady planar interface it is constant \({\tilde{\mathbf{V}}} = {\tilde{\mathbf{V}}}_{0}\) and can be chosen as the velocity of the inertial frame of reference \({\tilde{\mathbf{V}}}_{0} = {\mathbf{V}}_{0}\). For unsteady non-planar interface \({\tilde{\mathbf{V}}} \ne {\tilde{\mathbf{V}}}_{0}\) and

$$ {\tilde{\mathbf{V}}\mathbf{n}} = - \left. {{\mathbf{vn}}} \right|_{{\uptheta = 0}} = - \left. {{\tilde{\mathbf{j}}\mathbf{n}}/{\uprho } } \right|_{{\uptheta = 0}} $$
(1d)

We apply this form for the interface velocity in order to accentuate that the velocity of the inertial frame of reference \({\mathbf{V}}_{0}\), the velocity of steady planar interface \({\tilde{\mathbf{V}}}_{0}\), and the velocity of the unsteady perturbed interface \({\tilde{\mathbf{V}}}\) are distinct quantities. They may not be equal one another. The form of the interface velocity in equation Eq. (1d) permits us to reveal the unsteadiness of the interface velocity and the inertial stabilization mechanism. See for details the equation Eq. (3d) for the perturbed interface velocity in Sect. 2.2 below for thermally conducting fluids, and the works [17, 28,29,30] for ideal fluids.

Initial conditions are the initial perturbations of the interface and the flow fields. They define the dimensionality, the symmetry, the length-scale and the time-scale of the dynamics [17, 28,29,30].

Per the governing equations in the bulk and at the interface Eqs. (1a, 1b), in each fluid the thermal heat flux \(Q_{i}\) is collinear with the gradient \(\partial \left( {\chi e} \right)/\partial x_{i}\), and the interface is a level curve (level surface) of the internal energy function \(\chi e\) with constant \(\chi e\) at each side of the moving interface \({\uptheta} \left( {{\mathbf{x}},t} \right) \to 0^{ \pm }\). To our knowledge, these self-consistent boundary conditions for the thermal heat flux were never formulated and applied before.

2.2 Linearized dynamics

The general mathematical problem of interface dynamics is very challenging [12,13,14]. It can be simplified for constant \(\chi ,\;s\), by assuming that in each fluid to leading order the flow fields are uniform \(\left( {{\uprho } ,{\mathbf{v}},P,e} \right) = \left( {{\uprho } ,{\mathbf{v}},P,e} \right)_{0}\), the interface is planar \({\mathbf{n}} = {\mathbf{n}}_{0} ,\;{{\varvec{\uptau}}} = {{\varvec{\uptau}}}_{0}\), and the thermal heat flux at the interface is constant \({\mathbf{Q}}_{0}\). We slightly perturb the fields \({\uprho } = {\uprho }_{0} + \overline{{\uprho } },\;{\mathbf{v}} = {\mathbf{V}} + {\mathbf{u}},\;P = P_{0} + p,\;e = e_{0} + \overline{e}\), the interface \({\mathbf{n}} = {\mathbf{n}}_{0} + {\mathbf{n}}_{1} ,\;{{\varvec{\uptau}}} = {{\varvec{\uptau}}}_{0} + {{\varvec{\uptau}}}_{1}\), the enthalpy \(W = W_{0} + w\), and the fluxes of mass and heat \({\tilde{\mathbf{j}}} = {\mathbf{J}} + {\hat{\mathbf{j}}},\;{\mathbf{Q}} = {\mathbf{Q}}_{0} + {\mathbf{q}}\). For each quantity, the perturbation magnitude is small compared to its leading order value [17]. The leading order mass flux \({\mathbf{J}}\) is normal to the interface, with \({\mathbf{J}} \cdot {{\varvec{\uptau}}}_{0} = 0,\;{\mathbf{J}} \cdot {\mathbf{n}}_{1} = 0\).

To leading order the equations in the bulk and at the outside boundaries are obeyed; the fields are

$$ \left( {{\uprho } ,{\mathbf{v}},P,e,W,{\tilde{\mathbf{j}}},{\mathbf{Q}}} \right)_{h\left( l \right)} = \left( {{\uprho }_{0} ,{\mathbf{V}},P_{0} ,e_{0} ,W_{0} ,{\mathbf{J}},{\mathbf{Q}}_{0} } \right)_{h\left( l \right)} $$
(2a)

The boundary conditions at the interface are

$$ \begin{array}{*{20}l} {\left[ {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right] = 0,\quad \left[ {\left( {P_{0} + \frac{{\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)^{2} }}{{{\uprho }_{0} }}} \right)\,{\mathbf{n}}_{0} } \right] = 0,\quad \left[ {\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\left( {W_{0} + \frac{{{\mathbf{J}}^{2} }}{{2{\uprho }_{0}^{2} }}} \right) + \left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right)} \right] = 0} \hfill \\ {\left. {\left( {{\mathbf{Q}}_{0} \cdot {{\varvec{\uptau}}}_{0} } \right)} \right|_{{{\uptheta} = 0^{ + } }} = 0,\quad \left. {\left( {{\mathbf{Q}}_{0} \cdot {{\varvec{\uptau}}}_{0} } \right)} \right|_{{{\uptheta} = 0^{ - } }} = 0} \hfill \\ \end{array} $$
(2b)

For incompressible dynamics \(\left( {{\uprho }_{0} V^{2} /P_{0} } \right) \to 0,\;\left( {V^{2} /W_{0} } \right) \to 0\) free from thermal heat flux, the boundary conditions convert to \(\left[ {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right] = 0,\;\left[ {P_{0} \,{\mathbf{n}}_{0} } \right] = 0,\;\left[ {W_{0} } \right] = 0\) [12,13,14, 16, 17]. Thermal heat flux \(\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right) \ne 0\) seeds internal energy perturbations. We must retain \(\left[ {\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right) + \left( {{\mathbf{J}}^{2} /2{\uprho }_{0}^{2} } \right)\;\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)} \right] = - \left[ {\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\;W_{0} } \right]\) in the energy boundary condition even for sub-sonic dynamics.

To first order the governing equations in the bulk and the equation of state are

$$ \begin{array}{*{20}l} {\left( {\partial {/}\partial t + {\mathbf{V}} \cdot \nabla } \right)\overline{{\uprho } } + {\uprho }_{0} \nabla \cdot {\mathbf{u}} = 0,\;\;{\uprho }_{0} \left( {\partial {/}\partial t + {\mathbf{V}} \cdot \nabla } \right){\mathbf{u}} + \nabla p = 0,\;\;{\uprho }_{0} \left( {\partial {/}\partial t + {\mathbf{V}} \cdot \nabla } \right)\overline{e} + \nabla \cdot {\mathbf{q}} + P_{0} \nabla \cdot {\mathbf{u}} = 0} \hfill \\ {{\mathbf{q}} + \nabla \left( {\chi \overline{e}} \right) = 0,\;p{/}P_{0} = \overline{{\uprho } }{/}{\uprho }_{0} + \overline{e}/e_{0} } \hfill \\ \end{array} $$
(3a)

and the boundary conditions at the interface are:

$$ \begin{array}{*{20}l} {\left[ {\left( {{\mathbf{j}} + {\overline{\mathbf{j}}}} \right) \cdot {\mathbf{n}}_{0} } \right] = 0,\;\;\left[ {\left( {p + \;\frac{{2\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\left( {{\mathbf{j}} \cdot {\mathbf{n}}_{0} } \right)}}{{{\uprho }_{0} }} + \;\frac{{\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\left( {{\overline{\mathbf{j}}} \cdot {\mathbf{n}}_{0} } \right)}}{{{\uprho }_{0} }}} \right)\,{\mathbf{n}}_{0} } \right] = 0,\;\;\left[ {\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\left( {{\mathbf{J}} \cdot {{\varvec{\uptau}}}_{1} + {\mathbf{j}} \cdot {{\varvec{\uptau}}}_{0} } \right)\frac{{{{\varvec{\uptau}}}_{0} }}{{{\uprho }_{0} }}} \right] = 0} \hfill \\ {\left[ {\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\;\left( {w + {\mathbf{J}} \cdot {\mathbf{j}}/{\uprho }_{0}^{2} } \right) + \left( {{\mathbf{q}} \cdot {\mathbf{n}}_{0} } \right)} \right] = 0,\;\;\left. {\left( {{\mathbf{Q}}_{0} \cdot {{\varvec{\uptau}}}_{1} + {\mathbf{q}} \cdot {{\varvec{\uptau}}}_{0} } \right)} \right|_{{{\uptheta} = 0^{ - } }} = 0,\;\;\left. {\left( {{\mathbf{Q}}_{0} \cdot {{\varvec{\uptau}}}_{1} + {\mathbf{q}} \cdot {{\varvec{\uptau}}}_{0} } \right)} \right|_{{{\uptheta} = 0^{ - } }} = 0} \hfill \\ \end{array} $$
(3b)

where the mass flux perturbations are \({\hat{\mathbf{j}}} = {\mathbf{j}} + {\overline{\mathbf{j}}},\;{\mathbf{j}} = {\uprho }_{0} \left( {{\mathbf{u}} + {\mathbf{n}}_{0} \dot{{\uptheta} }} \right),\;{\overline{\mathbf{j}}} = \left( {\overline{{\uprho } }{/}{\uprho }_{0} } \right){\mathbf{J}}\). The perturbations of the flow fields decay far from the interface

$$ \left. {\left\{ {\overline{{\uprho } },{\mathbf{u}},p,\overline{e},{\mathbf{j}},{\overline{\mathbf{j}}},{\mathbf{q}},w} \right\}_{\;h} } \right|_{z \to - \infty } = 0,\quad \left. {\left\{ {\overline{{\uprho } },{\mathbf{u}},p,\overline{e},{\mathbf{j}},{\overline{\mathbf{j}}},{\mathbf{q}},w} \right\}_{\;l} } \right|_{z \to + \infty } = 0 $$
(3c)

The perturbed velocity of the interface is \({\tilde{\mathbf{V}}} = {\tilde{\mathbf{V}}}_{0} + {\tilde{\mathbf{v}}}\), \(\left| {{\tilde{\mathbf{v}}}} \right| < < \left| {{\tilde{\mathbf{V}}}_{0} } \right|\). Up to first order, it is

$$ {\tilde{\mathbf{V}}} = {\tilde{\mathbf{V}}}_{0} + {\tilde{\mathbf{v}}},\quad {\tilde{\mathbf{v}}\mathbf{n}}_{0} = - \left. {\left( {{\mathbf{un}}_{0} + \dot{{\uptheta} }} \right)} \right|_{{\uptheta} = 0} $$
(3d)

The schematics of the dynamics is illustrated by Fig. 1 (in a far field, not to scale), with the blue color marking the planar interface (dashed line) and the perturbed interface (solid line).

Fig. 1
figure 1

Schematics of the dynamics (in a far field, not to scale). Blue color marks the planar interface (dashed line) and the perturbed interface (solid line)

3 Results

3.1 Methodology

We develop new methodology for solving the complex multi-parameter problem Eqs. (13). We consider a sample case of two-dimensional flow periodic in the \(x\) direction, motionless in the \(y\) direction and spatially extended in the \(z\) direction. The interfacial function is \({\uptheta} = - z + z^{*} \left( {x,t} \right)\). The perturbed interface is \(z^{*} = Z^{*} \exp \left( {ikx + \Omega t} \right)\), with wavevector \(k = 2\pi {/}\lambda\) and wavelength \(\lambda\) set by the initial conditions. The perturbed velocity field can have potential and vortical components, \({\mathbf{u}} = \nabla \Phi + \nabla \times {{\varvec{\Psi}}}\). For the two-dimensional flow the vortical field and vorticity are \({{\varvec{\Psi}}} = \left( {0,\Psi ,0} \right),\;\nabla \times {\mathbf{u}} = \left( {0,\Delta \Psi ,0} \right)\).

3.1.1 Solutions’ structure

To identify the structure of the perturbation waves, we represent the flow fields as \(\left( {\Phi ,p,\Psi ,\overline{{\uprho } },\overline{e}} \right) = \left( {\hat{\Phi },\hat{p},\hat{\Psi },\hat{{\uprho } },\hat{e}} \right)\exp \left( {ikx - Kz + \Omega t} \right)\), where \(K\) is the wavevector in the direction of motion, to be found [12,13,14, 17, 30]. This representation reduces the perturbed equations in the bulk to the linear system \(M_{S} {\mathbf{z}} = 0\), where \({\mathbf{z}} = \left( {\Phi ,p,\Psi ,\overline{{\uprho } },\overline{e}} \right)^{\rm T}\) is a vector, and \(M_{S}\) is a \(5 \times 5\) matrix with components dependent on quantities \(\Omega ,k,K\), and parameters \(V,{\uprho }_{0} ,P_{0} ,e_{0} ,\chi ,s\). Matrix \(M_{S}\) is provided in Table 1. The condition \(\det M_{S} = 0\) yields a fifth order polynomial equation defining the wavevector(s) \(K\) and the structure(s) of the perturbation wave(s). Equation \(\det M_{S} = 0\) is provided in Table 2. Mathematically, we have to find the five fundamental solutions for the linear system \(M_{z} {\mathbf{z}} = 0\), including their eigenvalues and eigenvectors [14, 17, 28, 30, 32, 33, 56]. These five solutions plus the interface perturbation identify the six independent waves (degrees of freedom) to further solve the boundary value problem at the interface.

Table 1 The elements of the matrix defining the solutions’ structure in the bulk
Table 2 The polynomial equation \(det {M}_{S}=0\) and its coefficients

To find the five fundamental solutions, we first derive the characteristic equation \(\det M_{S} = 0\) for the wavevector \(K\) in the motion direction, as \(\det M_{S} = \left( {K - \Omega {/}V} \right)\left( {K^{4} + c_{3} K^{3} + c_{2} K^{2} + c_{1} K + c_{0} } \right)\) with coefficient \(c_{n} = c_{n} \left( {P_{0} ,s,{\uprho }_{0} ,V,K_{m} ,k,\Omega } \right)\). See Tables 1 and 2. We next find the five roots of this equation \(K = K_{n} ,n = 1,...,5\). By substituting each of these roots to the matrix \(M_{S}\), we then identify the relations between the amplitudes \(\left( {\hat{\Phi },\hat{p},\hat{\Psi },\hat{{\uprho } },\hat{e}} \right)\) of each of the perturbation waves. The five fundamental solutions \(\left( {\Phi ,p,\Psi ,\overline{{\uprho } },\overline{e}} \right)_{n} = \left( {\hat{\Phi },\hat{p},\hat{\Psi },\hat{{\uprho } },\hat{e}} \right)_{n} \exp \left( {ikx - K_{n} z + \Omega t} \right)\), \(n = 1,...,5\) define the structures of the perturbed flow fields in the bulk [32, 56].

The characteristic equation, as a fifth order polynomial equation, can in principle be solved for any \(\left( {P_{0} ,s,{\uprho }_{0} ,V,K_{m} ,k,\Omega } \right)\). In general case the solutions for the equation \(\det M_{S} = 0\) are extremely cumbersome, except for the wave representing the vortical field with \(K = \Omega {/}V\). In this work, we find analytical solutions for the weakly compressible dynamics with \({\uprho }_{0} V^{2} /P_{0} < < 1\). For these sub-sonic waves, three regimes exist. They are dominated by the processes of advection, diffusion and weak compressibility and are defined by the interplay of the wavevector \(k\) of the initial perturbation and the wavevector \(K_{m} = {\uprho }_{0} V/\chi\) set by thermal conductivity, as well as by the relation between the parameters \(K_{m} /k\) and \({\uprho }_{0} V^{2} /P_{0}\).

To our knowledge, the accurate analytical derivations of the structures of the flow fields in the bulk were never done before, including the sub-sonic dynamics \({\uprho }_{0} V^{2} /P_{0} < < 1\) in the regimes of advection, diffusion and weak compressibility. In order to derive the solutions, we apply the rigorous methods of applied mathematics and theoretical physics, including the regular perturbation method in the regimes of diffusion and weak compressibility with \(K_{m} {/}k < < 1\) and \({\uprho }_{0} V^{2} /P_{0} < < 1\), and the singular perturbation method in the regime of advection with \(K_{m} {/}k > > 1\) and \(P_{0} {/}{\uprho }_{0} V^{2} > > 1\) [6, 12, 14, 57].

3.1.2 Three regimes

These three regimes have four perturbation waves in common. We call them mechanical waves. The four mechanical waves describe the interface, the vortical velocity field of the light fluid, and the associated fields of pressure and the potential motion of the heavy and light fluids as:

$$ \begin{array}{*{20}l} {\left( {\overline{\Phi },\overline{p}} \right)_{h} = \hat{\Phi }_{h} \;\left( {1,\; - k{\uprho }_{0} V\left( {1 - \Omega {/}kV} \right)} \right)_{h} e^{ikx + kz + \Omega t} ,\;\left( {\overline{\Phi },\overline{p}} \right)_{l} = \hat{\Phi }_{l} \;\left( {1,\;k{\uprho }_{0} V\left( {1 + \Omega {/}kV} \right)} \right)_{l} e^{ikx - kz + \Omega t} } \hfill \\ {z^{*} = Z^{*} e^{ikx + \Omega t} ,\;\Psi_{l} = \hat{\Psi }e^{{ikx - \left( {\Omega /V_{l} } \right)z + \Omega t}} } \hfill \\ \end{array} $$
(4)

The vortical field does not contribute to pressure; hence it is energetic (rather than dynamic) in nature. Emphasize that while in this work we consider weakly compressible non-ideal fluids and we solve the \(5 \times 5\) linear system \(M_{z} {\mathbf{z}} = 0\), the structures of the four mechanical waves (one wave for the interface, two waves for the potential components of the velocities of the heavy and the light fluids, and one wave for the vortical field of the light fluid) are the same as those for the conservative dynamics in the incompressible ideal fluids free from the thermal heat flux (with or without surface tension) [17, 28, 30, 32, 33]. This ensures the accuracy of our results.

The other two waves are due to internal energy perturbations. We call them energetic waves. Their structure depends on the regime.

3.1.2.1 Advection

In the advection processes, the scale set by thermal conductivity is large compared to that of the initial perturbation, \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} \to \infty\). For \(\left( {P_{0} {/}{\uprho }_{0} V^{2} } \right)_{h\left( l \right)} \to \infty\), this leads to the perturbation waves

$$ \left( {\overline{e},\;\overline{{\uprho } }_{e} } \right)_{h} = \hat{e}_{h} \left( {1,\left( {{\uprho }_{0}^{2} s/P_{0} } \right)} \right)_{h} e^{{ikx + K_{mh} \left( {1 + s_{h} } \right)z + \Omega t}} ,\;\;\left( {\overline{e},\;\overline{{\uprho } }_{e} } \right)_{l} = \hat{e}_{l} \left( {1,\;\left( {{\uprho }_{0}^{2} s/P_{0} } \right)} \right)_{l} e^{{ikx - \left( {{\Omega \mathord{\left/ {\vphantom {\Omega {V_{l} }}} \right. \kern-\nulldelimiterspace} {V_{l} }}} \right)z + \Omega t}} $$
(5a)

These waves associate the fields of internal energy and density. The fields are attached to the interface in the heavy fluid, changing sharply for \(z \to z^{*}\) and vanishing for \(kz \to - \infty\). They are correlated with the vortical field in the light fluid, \(\sim\,e^{{ikx - \left( {\Omega /V_{l} } \right)z + \Omega t}}\).

In order to validate in experiments the solution’s structure in the advection process, one needs to diagnose the fields of the internal energy (e.g., temperature) and the density and the vortical field, and to observe the following. In both fluids the fields of the internal energy and density are in-phase with one another. Yet, in the heavy and the light fluids, they are asymmetric. In the light fluid, the fields of inertial energy and density are correlated with the vortical field; in the motion direction they have the length-scale \(\left( {V_{l} {/}\Omega } \right)\) set by the vortical field. In the heavy fluid, the fields of the inertial energy and density are attached to the interface and have the length-scale \(\left( {K_{mh} \left( {1 + s_{h} } \right)} \right)^{ - 1}\) that is substantially smaller than the wavelength of the initial perturbation, \(k\left( {K_{mh} \left( {1 + s_{h} } \right)} \right)^{ - 1} \to 0\). To our knowledge, these properties of the flow fields in the advection process were not accurately diagnosed before [6, 12,13,14].

3.1.2.2 Diffusion

In the diffusion processes, the scale set by thermal conductivity is small compared to that of the initial perturbation, \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} < < 1\). In the incompressible limit with \(\left( {\left( {{\uprho }_{0} V^{2} {/}P_{0} } \right){/}\left( {K_{m} {/}k} \right)} \right)_{h\left( l \right)} \to 0\), the perturbation waves are

$$ \begin{aligned} \left( {\overline{e},\;\overline{\Phi }_{e} ,\overline{p}_{e} } \right)_{h} & = \hat{e}_{h} \left( {1,\;{\uprho }_{0} V{/}K_{m} P_{0} ,\;{\uprho }_{0} \left( {k{/}K_{m} } \right)\left( {{\uprho }_{0} V^{2} {/}P_{0} } \right)\left( {1 - \Omega {/}kV} \right)} \right)_{h} \;e^{{ikx + K_{h} z + \Omega t}} ,\;K_{h} = k - \overline{k}_{h} \\ \left( {\overline{e},\;\overline{\Phi }_{e} ,\overline{p}_{e} } \right)_{l} & = \hat{e}_{l} \left( {1,\;{\uprho }_{0} V{/}K_{m} P,\; - {\uprho }_{0} \left( {k{/}K_{m} } \right)\left( {{\uprho }_{0} V^{2} {/}P_{0} } \right)\left( {1 + \Omega {/}kV} \right)} \right)_{l} \;e^{{ikx - K_{l} z + \Omega t}} ,\;K_{l} = k + \overline{k}_{l} \\ \end{aligned} $$
(5b)

with \(\bar k_h=(K_{mh}/2)(1+s_h)(1-\Omega/V_h k), \bar k_l=(K_{ml}/2)(1+s_l)(1+\Omega/V_l k)\) and wavevector \(K_{h} \approx k,K_{l} \approx k\). These waves associate the fields of inertial energy, velocity potential and pressure.

In order to validate in experiments the solution’s structure in the diffusion process, one needs to diagnose the fields of the internal energy and pressure and to observe the following. In both fluids the fields of the internal energy and pressure are correlated with one another. The fields are nearly symmetric in the heavy and the light fluids. In the direction of motion they have the length-scale \( \sim\,k^{{ - 1}} \) set by the initial perturbation. They are slightly modulated with spatial waves having length scales \(k_{h}^{ - 1}\) and \(k_{l}^{ - 1}\) in the heavy and light fluid respectively. The modulation length scales \(k_{h}^{ - 1}\) and \(k_{l}^{ - 1}\) are substantially greater than the wavelength of the initial perturbation, \(kk_{h\left( l \right)}^{ - 1} < < 1\), and they depend on the fluids thermal conductivities as \(K_{mh}^{ - 1}\) and \(K_{ml}^{ - 1}\), respectively, where \(\left( {K_{m}^{{}} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V{/}\chi } \right)_{h\left( l \right)}\) and \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} \to 0\). To our knowledge, these flow fields’ properties in the diffusion regime were not diagnosed before [6, 12,13,14].

3.1.2.3 Low Mach

In weakly compressible processes with vanishing thermal transport \(\left( {\left( {K_{m} {/}k} \right){/}{\uprho }_{0} V^{2} {/}P_{0} } \right)_{h\left( l \right)} \to 0\), in the limit \(\left( {{\uprho }_{0} V^{2} {/}P_{0} } \right)_{h\left( l \right)} \to 0\) the perturbation waves are

$$ \begin{aligned} \left( {\overline{e},\;\overline{\Phi }_{e} ,\;\overline{{\uprho } }_{e} ,\overline{p}_{e} } \right)_{h} & = \hat{e}_{h} \left( {1,\; - {\uprho }_{0} V{/}2kP_{0} + 2s{/}\left( {kV - \Omega } \right),\;{\uprho }_{0}^{2} s{/}P_{0} ,\;2{\uprho }_{0} s_{h} } \right)_{h} \;e^{{ikx + K_{h} z + \Omega t}} ,\;K_{h} = k + \tilde{k}_{h} \\ \left( {\overline{e},\;\overline{\Phi }_{e} ,\;\overline{{\uprho } }_{e} ,\overline{p}_{e} } \right)_{l} & = \hat{e}_{l} \left( {1,\;{\uprho }_{0} V{/}2kP_{0} - 2s{/}\left( {kV + \Omega } \right),\;{\uprho }_{0}^{2} s{/}P_{0} ,\;2{\uprho }_{0} s} \right)_{l} \;e^{{ikx + K_{l} z + \Omega t}} ,\;K_{l} = k + \tilde{k}_{l} \\ \end{aligned} $$
(5c)

with \(\tilde k_h = \left( {kV_{l} - \Omega } \right)^{2} \left( {{\uprho }_{0} {/}2kP_{0} } \right)_{h} ,\;\tilde{k}_{l} = \left( {kV_{l} + \Omega } \right)^{2} \left( {{\uprho }_{0} {/}2kP_{0} } \right)_{l}\) and with the wavevectors \(K_{h} \approx k,K_{l} \approx k\). These waves describe the associated fields of inertial energy, velocity potential, density and pressure.

In order to validate in experiments the solution’s structure in the weakly compressible process, one needs to diagnose the fields of internal energy, density and pressure and to observe the following. In both fluids the fields of internal energy, density and pressure are correlated with one another. The fields are nearly symmetric, with length scale \( \sim\,k^{{ - 1}} \) set by the initial perturbation, and they are slightly modulated with spatial waves having the length scales \(\tilde{k}_{h}^{ - 1}\) and \(\tilde{k}_{l}^{ - 1}\) in the heavy and the light fluid respectively. The modulation length scales \(\tilde{k}_{h}^{ - 1}\) and \(\tilde{k}_{l}^{ - 1}\) are substantially greater than the wavelength of the initial perturbation, \(k\tilde{k}_{h\left( l \right)}^{ - 1} < < 1\), and are independent of the fluids thermal conductivities. To our knowledge, these properties of the flow fields in the weakly compressible regime were not diagnosed before [6, 12,13,14].

3.1.3 Boundary value problem

In each process, to find the interface stability and the flow fields’ structure, we employ the perturbation waves representing the six independent degrees of freedom—4 common mechanical and 2 specific energetic—to solve the boundary value problem at the interface. The dynamics can be subject to a body force and an acceleration; the destabilizing acceleration \({\mathbf{g}}\) is directed from the heavy to the light fluid along the \(z\) direction of motion, \({\mathbf{g}} = \left( {0,0,g} \right)\), with constant \(g\), and modifies the pressure field.

The boundary value problem is thus reduced to the linear system \({\rm M}\,{\mathbf{r}} = 0\). The \(6 \times 6\) matrix \({\rm M}\) is defined by the interfacial boundary conditions, and the vector \({\mathbf{r}} = \left( {\Phi_{h} ,\Phi_{l} ,V_{h} z^{*} ,\Psi_{l} ,\overline{e}_{h} /kV_{h} ,\overline{e}_{l} /kV_{h} } \right)^{\rm T}\) is given by the six independent waves. The solution is \({\mathbf{r}} = C_{i} {\mathbf{r}}_{i}\), where \({\mathbf{r}}_{i}\) are fundamental solutions, and \(C_{i}\) are integration constants. By applying the condition \(\det \,{\rm M} = 0\) and by reducing the corresponding matrix to row-echelon form, we find the fundamental solution \({\mathbf{r}}_{i} \left( {\Omega_{i} ,{\mathbf{e}}_{i} } \right)\), including the eigenvalue \(\Omega_{i}\), the eigenvector \({\mathbf{e}}_{i}\), and the associated vector \({\hat{\mathbf{e}}}_{i}\) of the perturbation waves’ amplitudes.

We introduce dimensionless variables by employing \(1{/}k\) for the length scale, \(1/kV_{h}\) for the time-scale, and \(\left( {V{/}k,{\uprho }_{0} ,{\uprho }_{0} V^{2} } \right)_{h}\) for the velocity potentials, density, and pressure, respectively. We use the dimensionless values of the growth-rate \(\upomega = \Omega {/}kV_{h}\), the density ratio \(R = {\uprho }_{0h} {/}{\uprho }_{0l} ,R \ge 1\), the gravity magnitude \(G = g{/}kV_{h}^{2} ,G \ge 0\), and the thermal wavevector(s) \(\left( {k_{m} } \right)_{h\left( l \right)} = \left( {K_{m} {/}k} \right)_{h\left( l \right)}\). The thermal heat flux is scaled with the internal energy as \(\left( {{\mathbf{Q}}_{0} } \right)_{h\left( l \right)} = \left( {{\mathbf{J}}e_{0} \upvarepsilon_{0} } \right)_{h\left( l \right)}\). It defines the seeds \(\left( {\upvarepsilon_{0} } \right)_{h\left( l \right)}\) of the internal energy perturbations \(\left( {\upvarepsilon_{0} } \right)_{h\left( l \right)} = \left( {\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{J}}} \right)\;sMa^{2} /{\mathbf{J}}^{2} V^{2} } \right)_{h\left( l \right)}\). Here the modified Mach number is \(\left( {Ma^{2} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V^{2} /P_{0} } \right)_{h\left( l \right)}\). The leading order energy boundary condition \(\left[ {e_{0} \upvarepsilon_{0} + V^{2} {/}2} \right] = \left[ {W_{0} } \right]\) suggests \(\upvarepsilon_{0l} \le 0\) for \(\left[ {W_{0} } \right] = 0\) with \(\upvarepsilon_{0\,h} = 0\). Matrix \({\rm M}\) is \({\rm M} = {\rm M}\left( {\upomega ,R,G,\left( {k_{m} ,\upvarepsilon_{0} ,Ma,s} \right)_{h\left( l \right)} } \right)\). Its elements depend on \(\upomega\) and on the system parameters.

3.2 Joint dynamic properties

Advection, diffusion and weakly compressible dynamics have a number of properties in common. (1) They are degenerate. Each of them has only four fundamental solutions \({\mathbf{r}}_{1\left( 2 \right)\left( 3 \right)\left( 4 \right)}\) for six independent waves obeying six equations. This is due the thermodynamic nature of the internal energy perturbations \(\left( {\overline{e}} \right)_{h\left( l \right)}\), which are seeded by thermal heat flux, with \(\left( {\overline{e}} \right)_{h\left( l \right)} \sim\,\left( {\upvarepsilon_{0} } \right)_{h\left( l \right)}\). (2) For fundamental solutions \({\mathbf{r}}_{1\left( 2 \right)}\), the vortical and potential velocity fields are coupled with the interface perturbation and the thermal heat flux. For zero thermal heat flux, the velocity field is potential. (3) In each regime, solutions \({\mathbf{r}}_{3\left( 4 \right)}\) depend only on density ratio \(R\) and are the same as those of the conservative dynamics with zero thermal heat flux in ideal fluids [17, 28,29,30,31]. Particularly: solution \({\mathbf{r}}_{3}\) has zero perturbed velocity and pressure fields; solution \({\mathbf{r}}_{4}\) must have zero integration constant to obey the boundary conditions far from interface [17, 28,29,30,31, 73]. We omit their consideration here.

Solutions \({\mathbf{r}}_{1\left( 2 \right)}\) depend on the regime and on system parameters, with stable or unstable \({\mathbf{r}}_{1}\) and with always stable \({\mathbf{r}}_{2}\). We focus on the conservative dynamics (CD) under gravity (G) in the advection (A), diffusion (D), and low Mach (M) regimes. We set \({\mathbf{r}}_{CDG\left( A \right)\left( D \right)\left( M \right)} = {\mathbf{r}}_{1}\) with \({\hat{\mathbf{e}}}_{CDG\left( A \right)\left( D \right)\left( M \right)} = {\hat{\mathbf{e}}}_{1}\) for real \(\upomega_{CDG\left( A \right)\left( D \right)\left( M \right)}\), and \({\mathbf{r}}_{CDG\left( A \right)\left( D \right)\left( M \right)} = \left( {{\mathbf{r}}_{1} + {\mathbf{r}}_{1}^{*} } \right){/}2\), \({\mathbf{r}}_{1} = {\mathbf{r}}_{2}^{*}\), with \({\hat{\mathbf{e}}}_{CDG\left( A \right)\left( D \right)\left( M \right)} = \left( {{\hat{\mathbf{e}}}_{1} + {\hat{\mathbf{e}}}_{1}^{*} } \right){/}2\), for complex \(\upomega_{CDG\left( A \right)\left( D \right)\left( M \right)}\). While solutions \({\mathbf{r}}_{1\left( 2 \right)}\) are cumbersome, their forms are elegant for \(\left( {Ma} \right)_{h\left( l \right)} \to 0\).

3.3 Advection

In this process, the fluid velocities are \({\mathbf{u}}_{h} = \nabla \Phi_{h} ,\;{\mathbf{u}}_{l} = \nabla \Phi_{l} + \nabla \times {{\varvec{\Psi}}}_{l}\) with the field(s) \(\Phi_{h\left( l \right)} = \overline{\Phi }_{h\left( l \right)} ,\;{{\varvec{\Psi}}}_{l} = \left( {0,\Psi_{l} ,0} \right)\). The matrix is \({\rm M} = {\rm M}_{A}\). Matrix \({\rm M}_{A}\) is provided in Table 3. For \(\left( {Ma} \right)_{h\left( l \right)} \to 0\), the solution \({\mathbf{r}}_{CDGA}\) has the eigenvalue and the associated amplitude vector

$$ \upomega_{CDGA} = - F \pm i\sqrt R \sqrt {1 - G{/}G_{cr} - F^{2} {/}R + 2\;F\;k_{ml} \left( {1 + s_{l} R} \right)} ,\;{\hat{\mathbf{e}}}_{CDGA} = \left( {*,\;*,\;1,\;*,\;*,\;*} \right)^{\rm T} $$
(6a)

with \(G_{cr} = R\left( {R - 1} \right){/}\left( {R + 1} \right),\;F = - \upvarepsilon_{0l} {/}\left( {R - 1} \right)^{2} ,F \ge 0\).Asterisks marks functions on \(R,G,\left( {\upvarepsilon_{0} ,k_{m} } \right)_{h\left( l \right)}\). For \(\left( {k_{m} {/}\upvarepsilon_{0} } \right)_{l} \to 0\), the growth-rate is \(\upomega_{CDGA} \to - F \pm i\sqrt R \sqrt {1 - G{/}G_{cr} - F^{2} {/}R}\). For \(\left( {\upvarepsilon_{0} {/}k_{m} } \right)_{l} \to 0\) it is \(\upomega_{CDGA} \to \pm i\sqrt R \sqrt {\left( {1 + 2\;F\;k_{ml} \left( {1 + s_{l} R} \right)} \right) - G{/}G_{cr} }\). The dynamics \({\mathbf{r}}_{CDGA}\) is stable \({{Re}}\;\left[ {\upomega_{CDGA} } \right] \le 0\) for the acceleration values smaller than a threshold \(G < \tilde{G}_{cr}\) and is unstable \({{Re}}\;\left[ {\upomega_{CDGA} } \right] > 0\) for \(G > \tilde{G}_{cr}\). The critical value is \(\tilde{G}_{cr} = G_{cr}\) for \(\left( {k_{m} {/}\upvarepsilon_{0} } \right)_{l} \to 0\) and \(\tilde{G}_{cr} = G_{cr} \left( {1 + 2\;F\;k_{ml} \left( {1 + s_{l} R} \right)} \right)\) for \(\left( {\upvarepsilon_{0} {/}k_{m} } \right)_{l} \to 0\).

Table 3 The elements of the matrix defining solutions of the interfacial boundary value problem for the advection process

Dynamics \({\mathbf{r}}_{CDGA}\) couples the interface perturbation with potential and vortical components of the velocity fields and with the internal energy perturbations. At the interface, the velocity is shear free. In the bulk, the strength of the vortical field is determined by the thermal heat flux. For zero thermal heat flux the velocity field is potential in both fluids. For the unstable dynamics \({\mathbf{r}}_{CDGA}\) the interface velocity \({\tilde{\mathbf{V}}}\) increases with time.

Figure 2 illustrates flow fields for dynamics \({\mathbf{r}}_{CDGA}\) in the unstable regime, including the interface perturbation \(z^{*}\) and the perturbed fields of the velocity \({\mathbf{u}}\), the streamlines \({\mathbf{s}},\left( {\left( {d{\mathbf{s}}/dt} \right) \times {\mathbf{u}}} \right) = 0\), the pressure \(p\), the vorticity \(\nabla \times {\mathbf{u}} = \left( {0,\Delta \overline{\Psi }_{l} ,0} \right)\), and the inertial energy \(\overline{e}\). The parameters are similar to those in laser ablated plasmas with \(s = 2{/}3,\;\chi = ce^{\nu } ,\;\nu = 5{/}2\) and with \(R > > 1\) [37, 38, 58,59,60,61].

Fig. 2
figure 2

Flow fields for unstable accelerated advection: plots of real parts of perturbed fields of the interface, velocity, velocity streamlines, pressure, vorticity, internal energy, with red (blue) for positive (negative) in contour plots

Our theory of the unstable interfacial dynamics in the advection regime suggests that accurate experiments can diagnose in the future the following properties of the flow fields. In the light fluid, the internal energy field strongly correlates with the vortical field and has the same length-scale in the direction of motion. In the heavy fluid, the internal energy perturbation is attached to the interface and the velocity field is potential. The strengths of the fields of the internal energy and the vorticity depend on the thermal heat flux. The fields of the vorticity and pressure are independent from one another. By measuring the pressure fields in the bulk of each fluid far from the interface, one can capture the process of formation of bubbles and spikes at the interface. Here the bubble (spike) is the portion of the light (heavy) fluid penetrating the heavy (light) fluid, and it is ‘pushed’ from the higher to the lower pressure regions, Fig. 2.

3.4 Diffusion

In this process, the fluid velocities are \({\mathbf{u}}_{h} = \nabla \Phi_{h} ,\;{\mathbf{u}}_{l} = \nabla \Phi_{l} + \nabla \times {{\varvec{\Psi}}}_{l}\) with the field(s) \(\Phi_{h\left( l \right)} = \left( {\overline{\Phi } + \overline{\Phi }_{e} } \right)_{h\left( l \right)} ,\;{{\varvec{\Psi}}}_{l} = \left( {0,\Psi_{l} ,0} \right)\). The matrix is \({\rm M} = {\rm M}_{D}\), provided in Table 4. In incompressible limit \(\left( {Ma} \right)_{h\left( l \right)} \to 0\) for solution \({\mathbf{r}}_{CDGA}\) the eigenvalue and the amplitude vector are

$$ \upomega_{CDGD} = \pm i\sqrt R \sqrt {1 - G{/}G_{cr} - 2{{\left( {\left( {\upvarepsilon_{0} \left( {1 - k_{m} } \right)} \right)_{h} + \left( {\upvarepsilon_{0} \left( {1 + k_{m} } \right)} \right)_{l} } \right)} \mathord{\left/ {\vphantom {{\left( {\left( {\upvarepsilon_{0} \left( {1 - k_{m} } \right)} \right)_{h} + \left( {\upvarepsilon_{0} \left( {1 + k_{m} } \right)} \right)_{l} } \right)} {\left( {R - 1} \right)^{2} }}} \right. \kern-\nulldelimiterspace} {\left( {R - 1} \right)^{2} }}} ,\;{\hat{\mathbf{e}}}_{CDGD} = \left( {*,\;*,\;1,\;*,\;*,\;*} \right)^{\rm T} $$
(6b)

where asterisks mark functions on \(R,G,\left( {\upvarepsilon_{0} ,k_{m} } \right)_{h\left( l \right)}\). With \(\upvarepsilon_{0l} = - F\left( {R - 1} \right)^{2} - \upvarepsilon_{0h} ,\;F \ge 0\) and with vanishing \(\left( {k_{m} } \right)_{h\left( l \right)} \to 0\), the growth-rate is \(\upomega_{CDGD} = \pm i\sqrt R \sqrt {\left( {1 + 2F} \right) - G{/}G_{cr} }\). The dynamics \({\mathbf{r}}_{CDGD}\) is stable \({{Re}}\;\left[ {\upomega_{CDGD} } \right] \le 0\) for accelerations with magnitudes smaller than critical value \(G < \hat{G}_{cr}\). It is unstable \({{Re}}\;\left[ {\upomega_{CDGD} } \right] > 0\) for \(G > \hat{G}_{cr}\), where \(\hat{G}_{cr} = G_{cr} \left( {1 - 2\left( {\left( {\upvarepsilon_{0} \left( {1 - k_{m} } \right)} \right)_{h} + \left( {\upvarepsilon_{0} \left( {1 + k_{m} } \right)} \right)_{l} } \right)/\left( {R - 1} \right)^{2} } \right)\). For \(\left( {k_{m} } \right)_{h\left( l \right)} \to 0\) the critical value of the acceleration magnitude is \(\hat{G}_{cr} = G_{cr} \left( {1 + 2F} \right)\).

Table 4 The elements of the matrix defining solutions of the interfacial boundary value problem for the diffusion process

Dynamics \({\mathbf{r}}_{CDGA}\) couples the interface perturbation with the potential and vortical components of the velocity field and with the internal energy perturbations. While the velocity is shear free at the interface, the vortical field in the bulk is set by the thermal heat flux, and it vanishes for zero thermal heat flux, similarly to the advection process. For the unstable dynamics \({\mathbf{r}}_{CDGD}\) the interface velocity \({\tilde{\mathbf{V}}}\) increases with time. Figure 3 illustrates the perturbed flow fields \(z^{*} ,{\mathbf{u}},{\mathbf{s}},p,\nabla \times {\mathbf{u}},\overline{e}\) for the solution \({\mathbf{r}}_{CDGD}\) in the unstable regime.

Fig. 3
figure 3

Flow fields for unstable accelerated diffusion: plots of real parts of perturbed fields of the interface, velocity, velocity streamlines, pressure, vorticity, internal energy, with red (blue) for positive (negative) in contour plots

Our theory of the unstable interfacial dynamics in the diffusion regime suggests that accurate experiments can diagnose in the future the following properties of the flow fields. In the light and heavy fluids, the internal energy perturbations are nearly symmetric. They have the length scale being close to that of the initial perturbation, with slight departures set by the thermal conductivities. Their length scales are very different from the length scale of the vortical field. Similarly to the advection regime, in the diffusion regime, the strengths of fields of the internal energy perturbations and vorticity depend on the thermal heat flux; they are zero for zero heat flux. The vortical field and the pressure field are independent from one another. The pressure field defines the process of formation of bubbles and spikes at the interface, Fig. 3.

3.5 Low Mach

In this process, the fluid velocities are \({\mathbf{u}}_{h} = \nabla \Phi_{h} ,\;{\mathbf{u}}_{l} = \nabla \Phi_{l} + \nabla \times {{\varvec{\Psi}}}_{l}\) with the field(s) \(\Phi_{h\left( l \right)} = \left( {\overline{\Phi } + \overline{\Phi }_{e} } \right)_{h\left( l \right)} ,\;{{\varvec{\Psi}}}_{l} = \left( {0,\Psi_{l} ,0} \right)\). The matrix is \({\rm M} = {\rm M}_{M}\), provided in Table 5. For solution \({\mathbf{r}}_{CDGA}\) in the limit \(\left( {Ma} \right)_{h\left( l \right)} \to 0\) the eigenvalue and the amplitude vector are

$$ \upomega_{CDGM} = \pm i\sqrt R \sqrt {1 - G{/}G_{cr} - 2\left( {\upvarepsilon_{0h} + \upvarepsilon_{0l} } \right)/\left( {R - 1} \right)^{2} } ,\;{\hat{\mathbf{e}}}_{CDGD} = \left( {*,\;*,\;1,\;*,\;*,\;*} \right)^{\rm T} $$
(6c)

where asterisks mark functions on \(R,G,\left( {\upvarepsilon_{0} } \right)_{h\left( l \right)}\). With \(\upvarepsilon_{0l} = - F\left( {R - 1} \right)^{2} - \upvarepsilon_{0h} ,\;F \ge 0\) the growth-rate is \(\upomega_{CDGM} = \pm i\sqrt R \sqrt {\left( {1 + 2F} \right) - G{/}G_{cr} }\). For dynamics \({\mathbf{r}}_{CDGM}\) the interface is stable \({{Re}}\;\left[ {\upomega_{CDGM} } \right] \le 0\) for the acceleration magnitudes \(G < \overline{G}_{cr}\) and is unstable \({{Re}}\;\left[ {\upomega_{CDGM} } \right] > 0\) for \(G > \overline{G}_{cr}\), with \(\overline{G}_{cr} = G_{cr} \left( {1 + 2\;F} \right)\).

Table 5 The elements of the matrix defining solutions of the interfacial boundary value problem for the low Mach process

The fields of the dynamics \({\mathbf{r}}_{CDGM}\) are similar to that of \({\mathbf{r}}_{CDGD}\) with \(\left( {k_{m} } \right)_{h\left( l \right)} \to 0\). Figure 4 illustrates the perturbed fields \(z^{*} ,{\mathbf{u}},{\mathbf{s}},p,\nabla \times {\mathbf{u}},\overline{e}\) of the solution \({\mathbf{r}}_{CDGM}\) in the stable regime. Since for \(G < \overline{G}_{cr}\) the frequency is purely imaginary, \(\upomega_{CDGM} = \pm i\sqrt R \sqrt {\left( {1 + 2F} \right) - G{/}G_{cr} }\), the wavevector \(K = \upomega_{CDGM} {/}R\) of the vortical field is also imaginary. This leads to creation of the stable vortical patterns in the light fluid bulk, which are periodic in time and in space. Such solution can be implemented only when the boundary conditions away from the interface are somewhat noisy. Otherwise the integration constant for the stable solution \({\mathbf{r}}_{CDGM}\) must be zero, leading to zero perturbation fields and to constant interface velocity \({\tilde{\mathbf{V}}} = {\tilde{\mathbf{V}}}_{0}\).

Fig. 4
figure 4

Flow fields for stable accelerated low Mach dynamics: plots of real parts of perturbed fields of the interface, velocity, velocity streamlines, pressure, vorticity, internal energy, with red (blue) for positive (negative) in contour plots

Our theory of the unstable interfacial dynamics in the low Mach regime suggests that accurate experiments can diagnose in the future the following properties of the flow fields. In the light and heavy fluids, the internal energy perturbations are nearly symmetric. They have the length scale being similar to that of the initial perturbation, with slight departures independent from the thermal conductivities. In the stable regime, stable vortical patterns may appear in the light fluid bulk, and they are periodic in space and in time. The internal energy perturbations and the vortical structures arise only in the presence of the thermal heat flux. The vortical field and the pressure field are independent from one another. The pressure field determines the process of formation of bubbles and spikes at the interface, by pushing the bubble (spike) from the light (heavy) fluid to the heavy (light) fluid from the higher to the lower pressure regions, Fig. 4.

3.6 Interface dynamics and flow fields

Figure 5 illustrates the dependence of the growth-rate on the acceleration magnitude for the conservative dynamics with zero thermal heat flux and with the thermal heat flux in the advective, diffusive and low Mach regimes, and for the accelerated Landau–Darrieus and Rayleigh–Taylor instabilities [16, 17]. For large acceleration values, the growth rates of the conservative dynamics with zero heat flux \(\upomega_{CDG}\) and of the advective, diffusive and low Mach dynamics \(\upomega_{CDG\left( A \right)\left( D \right)\left( M \right)}\) are substantially greater than the growth-rates \(\upomega_{RT}\) and \(\upomega_{LDG}\) of Rayleigh–Taylor and Landau–Darrieus dynamics. The growth-rate \(\upomega_{CDG}\) is the largest and is asymptotic for \(\upomega_{CDG\left( A \right)\left( D \right)\left( M \right)}\); it is the envelope of \(\omega_{CDG\left( A \right)\left( D \right)\left( M \right)}\), Fig. 5

Fig. 5
figure 5

Dependence of the growth-rate(s) on the acceleration at some parameters values

In the notations used in this work, for the conservative dynamics, Rayleigh–Taylor dynamics and Landau–Darrieus dynamics free from the thermal heat flux, the growth-rates and amplitude vectors are:

$$ \begin{aligned} \upomega_{CDG} & = \pm i\sqrt R \sqrt {1 - G{/}G_{cr} } ,\quad {\hat{\mathbf{e}}}_{CDG} = \left( {*,\;*,\;1,\;0,\;0,\;0} \right)^{\rm T} ; \\ \upomega_{LDG} & = {{\left( { - R + \sqrt {R^{3} + R^{2} - R + G\left( {R^{2} - 1} \right)} } \right)} \mathord{\left/ {\vphantom {{\left( { - R + \sqrt {R^{3} + R^{2} - R + G\left( {R^{2} - 1} \right)} } \right)} {\left( {R + 1} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {R + 1} \right)}},\quad {\hat{\mathbf{e}}}_{LDG} = \left( {*,\;*,\;1,\;*,\;0,\;0} \right)^{\rm T} ; \\ \upomega_{RT} & = \sqrt {G\left( {R - 1} \right){/}\left( {R + 1} \right)} ,\quad {\hat{\mathbf{e}}}_{RT} = \left( {*,\;*,\;1,\;0,\;0,\;0} \right)^{\rm T} \\ \end{aligned} $$
(7)

For the conservative and Rayleigh–Taylor dynamics the velocity fields are potential, whereas Landau–Darrieus dynamics has the vortical field in the light fluid bulk. For the conservative and Landau–Darrieus dynamics the velocity is shear free at the interface, whereas Rayleigh–Taylor dynamics has the velocity shear at the interface [17, 28, 30, 32, 33].

This work investigates the dynamics of the interface with fluxes of heat and mass across it for non-ideal thermally conducting fluids. The problem differs dramatically from that of the interface dynamics in ideal fluids, Eqs. (13), Tables 1, 2, 3, 4, 5 [17, 28,29,30, 32, 33]. In the three regimes of advection, diffusion and low Mach, we find the fluid instabilities, which were never earlier discussed and which are defined by the interplay of the thermal heat flux, thermal conductivity, destabilizing acceleration and inertial stabilization, Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5.

We reveal nevertheless that both for non-ideal thermally conducting fluids with interfacial fluxes of heat and mass and for ideal fluids with zero heat flux, the primary mechanism of the interface stabilization is the macroscopic inertial mechanism, Eqs. (18), Figs. 1, 2, 3, 4, 5 [17, 28,29,30]: When the interface is perturbed, the fluid parcels follow the perturbation causing the change of momentum and energy. To conserve the momentum and energy, the interface as whole changes its velocity. Thus the reactive force occurs; the interface stability is defined by the interplay of inertia and buoyancy. The dynamics is unstable when the gravity exceeds the reactive force, Figs. 1, 2, 3, 4, 5. The thermal heat flux and the thermal conductivity provide with additional stabilization, Eqs. (18), Figs. 1, 2, 3, 4, 5.

The thermal heat flux and thermal conductivity are associated with the particle motions at kinetic scales, and are microscopic and thermodynamic in nature [14, 15]. They impact the dynamics quantitatively, by influencing the growth-rate, and qualitatively, by creating the vortical field in the bulk. For the stable and unstable dynamics, the velocity field is shear free at the interface, and the vortical structures in the bulk are energetic in nature, because they are caused by the thermal heat flux and the energy excess at the interface.

Advective, diffusive and low Mach dynamics \({\mathbf{r}}_{CDG\left( A \right)\left( D \right)\left( M \right)}\) have a number of common properties, such as inertial stabilization, destabilizing acceleration, shear free velocity at the interface, and volumetric vortical field. They have also important distinctions, including different structures of perturbation waves and flow fields, dispersion relations, growth-rates and critical accelerations, Eqs. (16), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5. Hence, one can deduce thermodynamic microscopic properties of a fluid system by diagnosing its macroscopic dynamics. These fluid instabilities constitute the novel class of the sub-sonic dynamics \(\left( {{\uprho }_{0} V^{2} {/}P_{0} } \right)_{h\left( l \right)} < < 1\) with the control parameter \(\left( {k_{m} } \right)_{h\left( l \right)} = \left( {K_{m} {/}k} \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V{/}\chi k} \right)_{h\left( l \right)}\) in the three regimes of advection \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} > > 1\), diffusion \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} < < 1\), and low Mach \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} \to 0\). For given values of experimental parameters \(\left( {V,{\uprho }_{0} ,P_{0} ,e_{0} ,\chi } \right)_{h\left( l \right)}\), one can control the type of the dynamics by managing the initial conditions and the value of the parameter \(\left( {K_{m} {/}k} \right)_{h\left( l \right)}\). One can transition from one regime to another by varying the wavevector \(k\) of the initial perturbation, and can further transition to the conservative dynamics for the vanishing thermal heat flux \(\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right)_{h\left( l \right)} \to 0\).

3.7 Range of applicability

Our general theoretical framework has the following range of applicability. (1) It focuses on weakly compressible dynamics with \(\left( {Ma} \right)_{h\left( l \right)} < < 1\) where \(\left( {Ma^{2} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V^{2} {/}P_{0} } \right)_{h\left( l \right)}\). (2) It can be used for a broad range of scales set by the thermal conductivity \(\left( {K_{m} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V{/}\chi } \right)_{h\left( l \right)}\) and by the initial conditions \(k\), including advection \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} > > 1\), diffusion \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} < < 1\) and low Mach \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} = 0\) regimes. (3) It can be applied for general equation of state \(P = s\;{\uprho } \;e\) with some \(s = const\) and for a broad range of values of thermal conductivity \(\chi\) including fluids with \(\chi = const\) and plasmas with \( \chi \,\sim\,e^{\nu } \left( {{\uprho } \;e} \right)^{\mu } \) and with constant \(\nu ,\mu\). (4) It can be employed for a broad range of the density ratios \(R = {\uprho }_{h} {/}{\uprho }_{l}\), including fluids with very similar \(\left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) \to 1^{ + }\) and very different \(\left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) \to \infty\) densities. (5) It is applicable for a broad range of values of the acceleration strength with \(g = GkV_{h}^{2} ,\;G \ge 0\) and thermal heat flux \(\left( {{\mathbf{Q}}_{0} } \right)_{h\left( l \right)} = \left( {{\mathbf{J}}e_{0} \upvarepsilon_{0} } \right)_{h\left( l \right)}\) with \(\left( {\upvarepsilon_{0} } \right)_{h\left( l \right)} = \left( {\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{J}}} \right)\;sMa^{2} {/}{\mathbf{J}}^{2} V^{2} } \right)_{h\left( l \right)}\) and \(\upvarepsilon_{0l} + \upvarepsilon_{0h} = - F\left( {R - 1} \right)^{2} ,\;F \ge 0\), as long as the perturbed dynamics are sub-sonic.

This range of applicability is expected to meet in the design of experiments and simulations comparing the observations in fluids, plasmas, materials with our theory. For experiments in fluids at low energy density conditions, for given values of the fluid parameters \(\left( {V,{\uprho }_{0} ,P_{0} ,e_{0} } \right)_{h\left( l \right)}\), the acceleration strength \(g\) and the initial perturbation wavevector \(k\), one may employ organic liquids in order to broadly vary the thermal conductivity \(\left( \chi \right)_{h\left( l \right)}\), the wavevector \(\left( {K_{m} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V{/}\chi } \right)_{h\left( l \right)}\) and the ratio \(\left( {K_{m} {/}k} \right)_{h\left( l \right)}\). For experiments in plasmas under conditions of high energy density, for given values of the plasma parameters \(\left( {{\uprho }_{0} ,P_{0} ,e_{0} ,\chi } \right)_{h\left( l \right)}\), the acceleration strength \(g\) and the initial perturbation wavevector \(k\), one may modify the velocity \(\left( V \right)_{h\left( l \right)}\), in order to broadly vary the wavevector \(\left( {K_{m} } \right)_{h\left( l \right)} = \left( {{\uprho }_{0} V{/}\chi } \right)_{h\left( l \right)}\) and the ratio \(\left( {K_{m} {/}k} \right)_{h\left( l \right)}\). For experiments in materials, for given values of the parameters \(\left( {V,{\uprho }_{0} ,P_{0} ,e_{0} ,\chi } \right)_{h\left( l \right)}\) and the acceleration strength \(g\), one may change the wavevector \(k\) of the initial perturbation in order to broadly vary the ratio \(\left( {K_{m} {/}k} \right)_{h\left( l \right)}\). Moreover, in each of these cases, variations of the initial perturbation wavevector \(k\) can lead to transitions between the advection, diffusion and low Mach regimes. These transitions can be revealed in the (well detected) qualitative observations of the flow fields’ structures and in the (well distinguished) quantitative measures of the values of the interface growth-rate and the characteristic length-scales.

This range of applicability can be met in a broad range of realistic processes in fluids, plasmas and materials, hence ensuring the value of our work in fundamental research and in practice [1].

3.8 Characteristic scales and initial conditions

By representing the growth-rates in dimensional form, we obtain characteristic scales of dynamics \({\mathbf{r}}_{CDG\left( A \right)\left( D \right)\left( M \right)}\), including wavevector \(k_{cr}\) at which the interface is stabilized, \({{Re}} \left[ \Omega \right] = 0\), and wavevector \(k_{\max }\), at which the growth-rate achieves its maximum, \(\partial \Omega {/}\partial k = 0,\partial^{2} \Omega {/}\partial k^{2} < 0\), and the associated length scales \(k_{{cr\left( {\max } \right)}}^{ - 1}\) and time scales \(\left( {V_{h} k_{{cr\left( {\max } \right)}} } \right)^{ - 1}\). For given values of parameters \(\left( {g,{\uprho }_{0} ,V,e_{0} ,s,\chi ,Q_{0} } \right)_{h\left( l \right)}\), with \(\left( {Q_{0} } \right)_{h\left( l \right)} = - \left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right)_{h\left( l \right)}\), we first set the steady velocity of the planar interface \({\tilde{\mathbf{V}}}_{0} = {\mathbf{V}}_{0} = - {\mathbf{V}}_{h}\) with \(V_{0} = \left| {{\mathbf{V}}_{0} } \right| = V_{h}\), as per the usual convention [5, 14, 16]. We next find that for fluids with very different densities \({\uprho }_{h} {/}{\uprho }_{l} \to \infty\), the growth-rate \(\hat{\Omega }\), the wave-vectors \(k_{cr} ,k_{\max }\), the maximum growth-rate of the accelerated dynamics \(\hat{\Omega }_{\max } = \left. {\hat{\Omega }} \right|_{{k = k_{\max } }}\), and the frequency of the inertial dynamics \(\hat{\Omega }_{in} = \left. {\hat{\Omega }} \right|_{g = 0}\) are remarkably similarly in the advective, diffusive, and low Mach processes. They are well captured by those of the conservative dynamics with zero thermal heat flux, Eq. (7) [17]:

$$ \hat{\Omega } = \sqrt {gk - \left( {kV_{0} } \right)^{2} \left( {\frac{{{\uprho }_{h} }}{{{\uprho }_{l} }}} \right)} ;\;\;\hat{\Omega }_{\max } = \frac{1}{2}\frac{g}{{V_{0} }}\sqrt {\frac{{{\uprho }_{l} }}{{{\uprho }_{h} }}} ;\;\;\hat{\Omega }_{in} = \pm ikV_{0} \sqrt {\frac{{{\uprho }_{h} }}{{{\uprho }_{l} }}} ;\;\;k_{cr} = \frac{g}{{V_{0}^{2} }}\left( {\frac{{{\uprho }_{l} }}{{{\uprho }_{h} }}} \right);\;\;k_{\max } = \frac{{k_{cr} }}{2} $$
(8)

For the conservative dynamics with the thermal heat flux, the energy equation is the advection–diffusion equation [6, 41, 42]. We find that for given fluids properties and for given acceleration and heat flux values \(\left( {{\uprho }_{0} ,V,e_{0} ,s,\chi ,g,Q_{0} } \right)_{h\left( l \right)}\), one may vary the character of the dynamics from advection to diffusion and to low Mach by varying the wavevector \(k\) of the initial conditions. This allows one to rigorously quantify the observer effect on the interface dynamics, which to our knowledge, was never done before [6, 13, 41,42,43].

For instance, for laser ablated plasmas in inertial confinement fusion the equation of state is \(P = s\;{\uprho } \;e\) with \(s = {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}\), the thermal conductivity is \(\chi \;\sim\;e^{\nu } \left( {{\uprho } \;e} \right)^{\mu }\) with \(\nu = 5{/}2,\mu = 0\), and the density ratio is \(\left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) > > 1\), whereas the thermal heat flux is finite, \(F\,\sim\,O\left( 1 \right)\). Under these conditions, for weak accelerations \(\left( {g{/}kV_{h}^{2} } \right) \to 0\) the dynamics is inertial and it experiences high frequency oscillations with \(\Omega \,\sim\,\hat{\Omega }_{in} = \pm ikV_{0} \sqrt {{\uprho }_{h} {/}{\uprho }_{l} }\), whereas for strong accelerations \(\left( {g{/}kV_{h}^{2} } \right) > > 1\) the instability growth-rate is \(\;\Omega \,\sim\,\hat{\Omega }_{\max } = \left( {g{/}2V_{0} } \right)\sqrt {{\uprho }_{l} {/}{\uprho }_{h} }\) and the ratio of the critical and maximum wavevector values is \(\left( {k_{cr} {/}k_{\max } } \right)\; = 2,\;k_{cr} = \left( {{\uprho }_{l} {/}{\uprho }_{h} } \right)g{/}V_{0}^{2}\) in agreement with our results Eq. (8) and with observations [37, 58,59,60,61].

3.9 Comparison with available experiments

Here we compare in more details our theoretical results with the observations in high energy density plasmas, including the inertial confidential fusion and the laboratory astrophysics [5, 34].

3.9.1 Plasma fusion

Inertial confinement fusion and stability of laser ablated plasmas is important application of our theory [5]. In order to accurately compare our theoretical results with experiments in plasma fusion, one needs to scrupulously analyze the unprocessed raw data gathered at high power laser facilities, such as the National Ignition Facility, and the Omega and the Nike laser facilities [37, 60, 61]. Since raw data are often a challenge to directly access, we compare the functional form of the growth-rate in our theory with that in other (quasi-empirical) models of ablative Rayleigh–Taylor and Richtmyer–Meshkov instabilities [58,59,60].

To conduct such comparisons, in our theory we consider fluids with very different fluid densities, \(R \to \infty\), with vanishing thermal heat flux, \(F \to 0\), and with the interface velocity set as \({\mathbf{V}}_{0} = - {\mathbf{V}}_{h}\) with \(V_{0} = \left| {{\mathbf{V}}_{0} } \right| = V_{h}\) as per the usual convention, and we associate this velocity magnitude with the rate of mass ablation \(\dot{m}\) as \(\dot{m} = {\uprho }_{h} V_{h}\) [5, 14, 16]. To compare with the functional form of the dispersion relation in other models, we use the growth-rate \(\hat{\Omega }\), the wave-vectors \(k_{cr} ,k_{\max }\), the maximum growth-rate of the accelerated dynamics \(\hat{\Omega }_{\max } = \left. {\hat{\Omega }} \right|_{{k = k_{\max } }}\), and the frequency of the inertial dynamics \(\hat{\Omega }_{in} = \left. {\hat{\Omega }} \right|_{g = 0}\) for the conservative dynamics in Eq. (8) [17, 28, 30]. We use the growth-rate of the conservative dynamics in Eq. (8) since for strong accelerations the growth-rate \(\omega_{CDG}\) is the envelope for the growth-rates \(\omega_{CDG\left( A \right)\left( D \right)\left( M \right)}\) in the advection, diffusion, and low Mach regimes, Fig. 5.

The pioneering models of ablative Rayleigh–Taylor instability in a single fluid [58] propose for the growth-rate \(\Omega_{BAK} = - kV_{0} + \sqrt {gk}\) (sub-script marks Bodner–Anisimov–Kull). For strong acceleration it approaches \(\left. {\Omega_{BAK} } \right|_{{\left( {g/kV_{0}^{2} } \right) \to \infty }} \to \sqrt {gk}\), in agreement with our results in Eq. (8) \(\left. {\hat{\Omega }} \right|_{{\left( {g/kV_{0}^{2} } \right) \to \infty }} \to \sqrt {gk}\).

The pioneering works of ablative Richtmyer–Meshkov instability [62] propose that for the vanishing acceleration \(g{/}kV_{0}^{2} = 0\) the growth-rate \(\Omega_{N}\) (sub-script marks Nishihara) has only the imaginary part, \({\text{Im}}\left[ {\Omega_{N} } \right] \ne 0\) and \({{Re}}\left[ {\Omega_{N} } \right] = 0\). This agrees with our results in Eq. (8), \(\left. {\hat{\Omega }} \right|_{{\left( {g/kV_{0}^{2} } \right) = 0}} = \hat{\Omega }_{in} = \pm ikV_{0} \sqrt {{\uprho }_{h} {/}{\uprho }_{l} }\).

The more recent models [59] propose for ablative Rayleigh–Taylor and Richtmyer–Meshkov instabilities with \(\left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) \to \infty\) the growth-rate \(\Omega_{SPI} = - 2kV_{0} + \sqrt {gk - \left( {kV_{0} } \right)^{2} \left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) + \left( {2kV_{0} } \right)^{2} }\) (sub-script marks Sanz–Piriz–Ibanez). For fluids with very different densities \(\left( {{\uprho }_{h} {/}{\uprho }_{l} } \right) \to \infty\) in the growth-rate \(\Omega_{SPI}\), the critical wavevector is \(k_{cr} \approx \left( {g{/}V_{0}^{2} } \right)\left( {{\uprho }_{l} {/}{\uprho }_{h} } \right)\), the maximum wavevector is \(k_{max } \approx {{k_{cr} } \mathord{\left/ {\vphantom {{k_{cr} } 2}} \right. \kern-\nulldelimiterspace} 2}\), the maximum growth-rate is \(\left. {\Omega_{SPI} } \right|_{{k = k_{max } }} \approx \left( {{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}} \right)\left( {g{/}V_{0} } \right)\sqrt {{\uprho }_{l} {/}{\uprho }_{h} }\), and the inertial dynamics frequency is \(\left. {\Omega_{SPI} } \right|_{{\left( {g/kV_{0}^{2} } \right) = 0}} \approx \pm ikV_{0} \sqrt {{\uprho }_{l} {/}{\uprho }_{h} }\) in agreement with our results in Eq. (8).

The models [60] propose for ablative Rayleigh–Taylor and Richtmyer–Meshkov instabilities with \(\left( {{\uprho }_{l} {/}{\uprho }_{h} } \right) \to \infty\) the growth-rate \(\Omega_{GAV} = - 2kV_{0} + \sqrt {gk - \left( {kV_{0} } \right)^{2} \left( {{\uprho }_{l} {/}{\uprho }_{h} } \right)}\) (sub-script marks Goncharov–Aglitsky–Velikovich). In the limiting case of very large density ratio \(\left( {{\uprho }_{l} {/}{\uprho }_{h} } \right) \to \infty\) in the growth-rate \(\Omega_{GAV}\), the critical wavevector is \(k_{cr} \approx \left( {g{/}V_{0}^{2} } \right)\left( {{\uprho }_{l} {/}{\uprho }_{h} } \right)\), the maximum wavevector is \(k_{max } \approx k_{cr} {/}2\), the maximum growth-rate is \(\left. {\Omega_{GAV} } \right|_{{k = k_{max } }} \approx \left( {1{/}2} \right)\left( {g{/}V_{0} } \right)\sqrt {{\uprho }_{l} {/}{\uprho }_{h} }\) and the frequency of the inertial dynamics is \(\left. {\Omega_{GAV} } \right|_{{\left( {g/kV_{0}^{2} } \right) = 0}} \approx \pm ikV_{0} \sqrt {{\uprho }_{l} {/}{\uprho }_{h} }\). They all agree with our results in Eq. (8).

By accurately capturing dispersion relations in the existing quasi-empirical models [58,59,60, 62], our theory agrees with available experimental data to the same or even to the greater extent [37, 58,59,60,61]. Our theory is applicable in a broad range of conditions, \(\left( {{\uprho }_{l} {/}{\uprho }_{h} } \right) \in \left( {1,\infty } \right)\) and \(\left( {K_{m} {/}k} \right)_{h\left( l \right)} \in \left[ {0,\infty } \right)\), and identifies important properties of the interface dynamics not grasped by other models [57,58,59,60]. See for details Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5.

3.9.2 Laboratory astrophysics

Our theory [1, 4, 17] can be applied to explain the scaled laboratory experiment [34]. The experiments are designed to investigate the process of the interfacial mixing of matter in supernova remnants (SNR) caused by supernova blasts [34]. The experiment focuses on the effect of the thermal heat flux on the evolution of the unstable (fluid) interfaces in high energy density (neutral) plasmas, and observes significant differences between the dynamics in the high and low flux cases. By analyzing the available processed experimental data, we find that in the experiments [34] the unstable interfacial dynamics is driven by variable acceleration decaying with time as \(g\,\sim\,t^{ - 0.8}\). Furthermore, in the high flux case, when compared to low flux case, (1) the bubbles and spikes both propagate faster, (2) the growth of the interface perturbation (the difference between the bubble and spike positions) is slower, and (3) the vortical structures on the sides of evolving spikes are absent.

Our theory elegantly explain these intriguing puzzles, Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5. First, for decaying accelerations the acceleration strength is generally expected to be low, \(G\,\sim\,G_{cr}\). This can lead to higher growth-rate of the perturbation amplitude in Rayleigh–Taylor instability (corresponding to the low flux case) when compared to the interface dynamics with interfacial fluxes of mass and heat in the high flux case [34], Fig. 5. Note that the latter is stabilized by the inertial mechanism and the thermal conductivity, Fig. 5. Second, according to our theory, the inertial stabilization mechanism can also cause the growth of the interface velocity, Eqs. (18). This can explain the faster propagations of bubbles and spikes in the high flux case [34]. Third, the interface dynamics with interfacial fluxes of mass and heat is free from the interfacial shear and thus cannot have shear-driven vortical structures at the interface, in contrast to Rayleigh–Taylor instability, Eqs. (18). This can explain the absence of vortical structures in the high flux case observed in the experiment [34].

3.10 Thermal heat flux at continuous and at kinetic scales

The thermal heat flux seeds the internal energy perturbations through the interfacial boundary conditions Eqs. (13). The boundary conditions imply that in each fluid the thermal heat flux is normal to the inertial energy level curves (level surfaces), thus enabling the treatment of the interface as a discontinuity. This rigorous approach is applicable in a broad range of conditions, including very thin (a few nanometers) interfaces with non-diffusive heat and mass transports, which are observed in the experiments and investigated in the molecular dynamics simulations [7,8,9,10].

While the thermal heat flux is usually associated with the enthalpy of formation, since its values are known and tabular [12,13,14], our analysis suggests that a more cautious consideration is needed. Particularly, the thermal heat flux at the interface is set by the balance of the physics enthalpy and the specific kinetic energy as \(\left[ {\left( {{\mathbf{Q}}_{0} \cdot {\mathbf{n}}_{0} } \right)} \right] = - \left[ {\left( {{\mathbf{J}} \cdot {\mathbf{n}}_{0} } \right)\;\left( {W_{0} + \left( {{\mathbf{J}}^{2} {/}2{\uprho }_{0}^{2} } \right)} \right)} \right]\). The enthalpy \(W_{0}\) and the enthalpy of formation \({\overline{W}}_{0}\) are related as \(W_{0} = {\overline{W}}_{0} + C_{P} {\Theta}\), where the specific heat at constant pressure is \(C_{P}\) and the temperature is \({\Theta}\). To obtain the value of the thermal heat flux at the interface, one needs to accurately evaluate the amount of the specific (per unit mass) energy, which is left in the system after the deduction of the specific enthalpy of formation \(\overline{W}_{0}\), the specific thermal energy \(C_{P} {\Theta}\) and the specific kinetic energy \(\left( {{\mathbf{J}}^{2} {/}2{\uprho }_{0}^{2} } \right)\). This consideration is consistent with molecular dynamics simulations unambiguously showing complexity of thermal heat flux in energetic materials and chemically reactive systems [7,8,9, 18,19,20, 63].

Our theory finds that by linking micro to macro scales, the interface is a place where balances are achieved. For a given system, by modifying thermal heat fluxes, initial conditions and accelerations, one can impact the dynamics of flow fields and the interface velocity and stability.

3.11 Design of experiments

Stability of an interface with heat and mass fluxes is a long-standing problem in mathematics, science and engineering. While the pioneering work Landau 1994 found that the interface with interfacial mass flux is unconditionally unstable, experiments were challenged to directly observe the instability in nearly ideal incompressible fluids. A consensus was achieved that at small length scales the interface can be stabilized by microscopic mechanisms, including, e.g., viscosity, compressibility, surface tension, and thermal conductivity [6, 12,13,14, 43, 44].

Note that while in laboratory the phase boundaries are challenging to accurately study at large scales, in nature the inertial dynamics of the interface with interfacial mass flux is stable at global scales, in agreement with our theory Eqs. (18) [17, 28,29,30, 32, 33]. This stability can be observed in geophysical multi-phase flows. For instance, the waters of the Pacific and Atlantic oceans meet and not mix at global scales, and so do the waters of the Green and Colorado rivers [64].

Our theory finds that both for ideal and realistic fluids the interface stability is set primarily by the macroscopic inertial mechanism balancing the destabilizing acceleration, whereas the microscopic thermal conductivity and the thermal heat flux provide with an additional stabilization and create vortical fields. The interface velocity is not a constant value. For the unstable dynamics, the growth of the interface perturbations is accompanied by the growth of the interface velocity. For non-ideal thermally conducting fluids, the fluid instabilities constitute the novel class in the advection, diffusion and low Mach regimes, with distinct values of the growth-rates and the flow fields’ structures, and with transitions between the regimes occurring upon variations of the parameters, including the initial conditions. These theoretical results can aid to better understanding the interface dynamics in nature and technology, and to advancing the design and diagnostics of laboratory experiments.

In realistic fluids, accurate quantification of interface dynamics requires non-intrusive diagnostics of fast events and ultra-high performance in space–time resolution, bandwidth, data-acquisition rate, and control of initial and experimental conditions [1, 65]. The requirements can be implemented in integrated experimental systems incorporating the state-of-the-art in motion control, precision mechanics, optical imaging, image processing, digital signal processing, as well as data the analysis methods [65]. In order to simultaneously measure in space and in time the flow fields, experimental diagnostics can include particle image velocimetry, planar laser-induced fluorescence, holographic particle image velocimetry, and molecular tagging velocimetry [64,65,66]. By employing these metrological capabilities, one can implement and accurately measure multi-phase flows in a relatively small laboratory form factor and provide data suitable for a direct comparison with our rigorous analysis. In addition to the interface growth and growth-rate, the diagnostics can include the structure of the macroscopic fields in the bulk, the properties of microscopic transports at the interface, and the unsteadiness of the interface velocity, as suggested by our theory, Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5 [1, 17, 28,29,30, 32, 33, 64].

For plasmas with high energy densities produced at high power laser facilities, one can apply the advanced capabilities in the pulse shaping and the target fabrication, in order to control the thermal heat flux and the initial conditions and to observe the transitions between the advection, diffusion and low Mach regimes [37, 60, 61, 67]. For low energy density plasmas, one can use the advanced optical methods to quantify the microscopic structures and transports at the interface [38], and the state-of-the-art metrology to measure properties of macroscopic fields in the bulk [68]. Particularly, the recent experiments [68] on the exploding water plasma can serve as test-bed for studies of unstable hydrodynamics. These affordable experiments (with tight control of parameters, high repeatability and high reproducibility) can provide data on laboratory plasmas with unstable interfaces, and, in synergy with our theory, can complement the existing approaches on control of high and low energy density plasmas [37, 38, 60, 61, 68].

In materials, in order to better understand complex physical and chemical non-equilibrium processes and the interfacial transport at microscopic scales, one can apply highly accurate experiments [10, 63] and Lagrangian numerical methods, including reactive molecular dynamics, molecular dynamics, particle-in-cell and smoothed particle hydrodynamics [7,8,9, 18,19,20, 48, 69]. These approaches can quantify with high fidelity the microscopic processes and non-diffusive interfacial transports occurring at phase boundaries far from thermal equilibrium. In synergy with our theory they can provide the new insight into the dynamics of unstable interfaces in technological applications, such as purification of water, electro-catalysis, explosions of energetic materials, and nano-fabrication [1].

4 Discussion

Interfaces and interfacial mixing, and their non-equilibrium dynamics are a corner-stone problem of science, mathematics and engineering with a broad range of applications in nature, technology, industry [1, 2]. We focused on the theoretical problem of the interface dynamics with interfacial fluxes of heat and mass in non-ideal thermally conducting weakly compressible fluids and developed the general and original theoretical framework to rigorously solve the problem in a broad range of conditions. The interface’s stability and flow fields’ structure were identified in the three regimes not discussed before.

Our methodology addresses the long-standing challenges [12,13,14]. (1) We formulated the self-consistent boundary conditions, Eqs. (13). At the interface, in addition to the conserved fluxes of mass, momentum and energy, the thermal heat flux is normal to level curves (surfaces) of the internal energy function in each fluid. This enables the interface treating as a discontinuity in a broad range of conditions. (2) We accurately evaluated the thermal heat flux, beyond the traditional use of the enthalpy of formation, Eqs. (13) [12,13,14]. This macroscopic consideration is consistent with microscopic properties, as illustrated by the experiments and molecular dynamics simulations [8,9,10, 18, 20, 63]. (3) We rigorously derived the six independent perturbation waves in the advection, diffusion and low Mach regimes, Eqs. (46). Four of these waves are mechanical. Two others are thermodynamic and are due to internal energy perturbations seeded by the thermal heat flux. Five of these waves capture the perturbations of the flow fields in the bulk. One wave describes the interface perturbations. See Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5.

Our theory finds fundamental properties of the interface stability and the flow fields’ structure in the regimes not identified before, Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5 [12,13,14, 17]. (I) The primary stabilization is due to the macroscopic inertial mechanism caused by the motion of the interface as a whole. The interface can be destabilized by strong enough accelerations. Otherwise, the interface is stable at global scales. (II) The thermal heat flux and the microscopic thermodynamics lead to formation of the vortical field in the bulk. The strength of the vortical field is defined by the thermal heat flux. For the stable and the unstable dynamics, the volumetric vortical structures are energetic in nature, and the velocity field is shear free at the interface. (III) The interface is the place where balances are achieved. One can deduce the microscopic thermodynamics properties of the system by diagnosing at macroscopic scales the interface and the flow fields’ dynamics. By varying the thermal heat flux, thermal conductivity, initial conditions, and acceleration, one can impact the interface stability and the flow fields’ structure, Eqs. (18), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5.

Our theory explains the interface stability found in a broad range of conditions, including multiphase geophysical flows and combustible systems [6, 12,13,14, 64]. It accurately captures dispersion curve of ablative Rayleigh–Taylor instabilities in plasmas with highly contrasting fluid densities [37, 58,59,60,61]. The increase of the velocity of an unstable interface, in addition to the growth of interface perturbations, explains the quick extinction of hot spot in inertial confinement fusion observed in the experiments [5, 37]. Our results are consistent with experiments and simulations investigating the complexity of the thermal heat flux for energetic materials and for interfaces with non-diffusive transports [8,9,10, 18,19,20, 48, 63, 70].

Our analysis directly links the non-equilibrium dynamics and kinetics—which are the system macroscopic dynamics and its microscopic thermodynamic properties—and yields the qualitative and quantitative benchmarks not diagnosed earlier [12,13,14]. They include, e.g., the existence of the three regimes—advection, diffusion, and low Mach; the dependence of the interface stability, the flow fields’ and the characteristic scales on the system parameters; the interplay of the thermal heat flux and the acceleration with the initial conditions. These results indicate a need in further advancement of the methods for numerical modeling and experimental diagnostics in fluids, plasmas and materials, in high and in low energy density regimes [1,2,3,4,5,6,7,8,9,10, 18, 20, 37, 38, 48, 60, 61, 63, 65,66,67,68,69,70,71,72].

Important outcomes of our work in applied sciences are the following. (1) Fluid interfaces (phase boundaries) with interfacial mass flux are stable at global scales unless they are destabilized by strong accelerations. (2) The stabilization is primarily due to macroscopic inertial mechanism whereas the microscopic effects lead to formation of vortical structures in the bulk. (3) For the accelerated interface, the new fluid instabilities develop, which have dispersion properties not discussed before, and for which the growth of the interface perturbations is accompanied with the growth of interface velocity, Eqs. (17), Figs. 1, 2, 3, 4, 5 and Tables 1, 2, 3, 4, 5.

Our theoretical approach self-consistently defines the flow fields’ structure, identifies mechanical and thermal degrees of freedom, and partitions the velocity field into potential and vortical components. This representation allows us to systematically study the nonlinear and self-similar dynamics of unstable interfaces in a broad range of conditions. Specifically, through synergy of this approach with group theory approach [4, 21,22,23,24,25,26,27] (the latter works remarkably well for Rayleigh–Taylor (RT) and Richtmyer–Meshkov (RM) unstable fronts), we can study nonlinear and self-similar dynamics of unstable interfaces with interfacial fluxes of mass and heat in the advection, diffusion and low-Mach regimes, in two and in three-dimensional flows, for various symmetries of the initial conditions, and for constant and variable accelerations. We address the investigation and the solution of this titanic task to the future.

Our approach can be applied for studies of non-equilibrium dynamics of interfaces in a broad range of processes in nature and technology, including and not limited to the supernovae remnants, the multiphase dynamics of geophysical flows, the ablation front instabilities in fusion plasmas, the detonation of energetic materials, the D’yakov–Kontorovich instabilities of shock waves, the Stefan problem of the evolution of matter undergoing phase transition, the realistic turbulent processes, the transportation security of liquefied natural gas, the purification of water, the fluid transports in fuel cells and micro-channels, the electro-catalysis and nano-fabrication [1,2,3,4,5,6,7,8,9,10, 18, 20, 34, 36,37,38,39,40, 42, 44,45,46,47,48,49,50,51,52,53,54,55, 60, 61, 63, 68, 70,71,72,73,74,75].

5 Conclusion

This work examined the classical long-standing problem of stability of a phase boundary—the fluid interface with heat and mass fluxes across it for non-ideal thermally conducting nearly incompressible fluids. We developed a rigorous theory resolving challenges not addressed before and reported key discoveries regarding three regimes—with different flow field structures, the interplay of inertial mechanism with thermal heat flux and destabilizing acceleration, the coupling of macro to micro scales—to define interfaces where balances are achieved. We explored the novel class of the fluid instabilities and revealed the parameter controlling transitions between the regimes through varying the initial conditions. Our theory paves the path to grasp a broad range of processes in nature and technology, ranging from celestial events of supernova remnants and formation of stars to the atomic level of fusion and electro-catalysis.