This work focuses on the long-standing problem of inertial dynamics of an interface with interfacial mass flux and reports new mechanisms for the interface stabilization and destabilization. The interface is a phase boundary separating fluids of different densities and having interfacial mass flux. To analyze the interface dynamics from a far field, we develop and apply the general matrix method to rigorously solve the boundary value problem involving the governing equations in the fluid bulk and the boundary conditions at the interface and at the outside boundaries of the domain. We find the fundamental solutions for the linearized system of equations and analyze the interplay of interface stability with flow fields’ structure by directly linking rigorous mathematical attributes to physical observables. We find that the interface is stable when the dynamics conserves the fluxes of mass, momentum, and energy; the stabilization is due to an inertial mechanism causing small oscillations of the interface velocity. In the classic Landau’s dynamics, the postulate of perfect constancy of the interface velocity leads to the development of Landau–Darrieus instability. This destabilization is also linked to the imbalance of the perturbed energy at the interface. The classic Landau’s solution is found to have degeneracy; lifting of the degeneracy may lead to singularity and self-similar dynamics. Our results compare well with traditional theories of combustion and propose new experiments to study the dynamics of the interface and the flow fields in combustible systems. We further conduct reactive molecular dynamics simulations to elucidate the complexity of chemical processes, to study the destabilizing effect of energy fluctuations on the interface stability, and to illustrate the chemistry-induced instabilities. In summary, we identify the extreme sensitivity of the interface dynamics to the interfacial boundary conditions, including the formal properties of fundamental solutions and the qualitative and quantitative properties of the flow fields. This provides new opportunities for studies, diagnostics, and control of multiphase flows in a broad range of processes in nature and technology.

## I. INTRODUCTION

Interfacial mixing and transport control a broad range of natural and technological processes in fluids, plasmas, and materials from celestial to atomistic scales (Abarzhi and Goddard, 2019; Abarzhi *et al.,* 2013; Zeldovich and Raiser, 2002; and Landau and Lifshitz, 1987). Examples include supernovae and fusion, planetary convection and reactive fluids, formation of phases in liquefied natural gas, and material transformations under impact (Remington *et al.,* 2019; Zhakhovky *et al.,* 2019; Schlimpert *et al.,* 2016; Arjomandnia *et al.,* 2014; Miglani *et al.,* 2014; Abarzhi, 2010; Azechi *et al.,* 2007; Bell *et al.,* 2004a; 2004b; Stein and Norlund, 2000; Rast, 1998; Arnett, 1996; and Prosperetti and Plesset, 1984). In realistic circumstances, material transport is often characterized by sharply and rapidly changing flow fields and by relatively small effects of dissipation and diffusion. This may lead to the formation of interfaces (discontinuities) separating the flow phases (non-uniformities) at macroscopic (continuous) scales (Abarzhi *et al.,* 2019; 2015; Anisimov *et al.,* 2013; Abarzhi, 2010; Kadau *et al.,* 2010; Meshkov, 2006; and Zeldovich and Raizer, 2002). Here, we consider from a far field the evolution of a hydrodynamic discontinuity that separates ideal incompressible fluids of different densities and has the interfacial mass flux (Landau and Lifshitz, 1987). We develop the general matrix method to rigorously solve the linearized boundary value problem and directly link the interface stability to the flow fields’ structure. We find extreme sensitivity of the flow dynamics to the interfacial boundary conditions, identify new mechanisms of the interface stabilization and destabilization, and thus reveal new opportunities for studies of interfacial dynamics in multiphase flows (Abarzhi *et al.,* 2019).

This paper is organized as follows: We start with introduction in Sec. I and present the method in Sec. II, including the governing equations (Sec. II A), linearized dynamics (Sec. II B), and fundamental solutions (Sec. II C). Results in Sec. III present the investigations of the conservative system (Sec. III A), the classic Landau’s system (Sec. III B), and the dynamic Landau’s system (Sec. III C). The properties of interfacial dynamics are discussed in Sec. IV, including stable and unstable solutions (Sec. IV A), the mechanisms of the interface stabilization and destabilization (Sec. IV B), and the effect of energy fluctuations (Sec. IV C). Comparison of our results with traditional studies of combustible systems is given in Sec. V, including the outlines of traditional studies (Sec. V A) and the new experiments (Sec. V B). Chemically reactive systems are analyzed in Sec. VI, including reactive molecular dynamics (Sec. VI A), the RMD2Kin method (Sec. VI B), and the chemistry-induced instability (Sec. VI C). Interface dynamics influenced by acceleration is briefly compared with our results for inertial dynamics in Sec. VII. Multiphase flows outcomes of our theory are outlined in Sec. VIII. Discussion and conclusion are provided in Secs. IX and X, respectively.

In hydrodynamics, the evolution of an interface with mass flow across it is usually studied in premixed combustion and in fluid vaporization (Abarzhi *et al.,* 2019; Zhakhovsky *et al.,* 2019; Prosperretti, 2017; Klimenko and Class, 2000; Class *et al.,* 2003a; 2003b; Peters, 2000; Mayo *et al.,* 1990; Sivashinsky, 1983; 1980; Williams, 1965; and Zeldovich, 1944). Similar problems are considered in plasmas (stability of ablation fronts in inertial confinement fusion), astrophysics (thermonuclear flashes on the surface of compact stars), materials science (melting and evaporation of materials under high pressures and high strain rates), and industrial applications (atomization of liquid jets in scramjets) (Remington *et al.,* 2019; Zhakhovsky *et al.,* 2019; 2014; Hurricane *et al.,* 2014; Anisimov *et al.,* 2013; Abarzhi, 2010; Kadau *et al.,* 2010; Arnett *et al.,* 2009; Gorokhovski and Herrmann, 2008; Azechi *et al.,* 2007; Marmottant and Villermaux, 2004; Bell *et al.,* 2004a; 2004b; Hillebrandt and Niemeyer, 2000; Cherkov and Yakhot, 1998; Piriz *et al.,* 1997; Wu *et al.,* 1997; and Sanz, 1994). The classic theoretical framework for the problem was developed by Landau (1944). It analyzed the dynamics of a discontinuous interface separating ideal incompressible fluids of different densities. By balancing at the interface the fluxes of mass and momentum and by implementing a special condition for the perturbed mass flux, Landau (1944) found that the interface was unconditionally unstable: small perturbations of the interface grew exponentially with time, leading to the development of Landau–Darrieus Instability (LDI) (Darrieus, 1938).

To connect the classical framework (Landau, 1944) to realistic environments, several theoretical approaches have been developed. In high energy density plasmas, significant deviations have been detected between the growth-rates of ablative Rayleigh–Taylor (RT) and accelerated Landau–Darrieus (LD) instabilities (Hurricane *et al.,* 2014; Sanz 1944; and Kull and Anisimov, 1985). In reactive and super-critical fluids, stabilizing influences have been found for dissipation, diffusion, and finite interface thickness on the interfacial dynamics at small scales (Class *et al.,* 2003a; 2003b; Peters, 2010; 2000; Clavin, 2000; Prosperetti and Plesset, 1984; Pelce and Clavin, 1982; Matkowsky and Sivashinsky, 1979; Williams, 1971; Markstein, 1951; and Zeldovich, 1944). More recently, it has been found that the interface stability is sensitive to the flux of energy fluctuations produced by the perturbed interface (Abarzhi *et al.,* 2015). These approaches successfully expanded the classical framework (Landau, 1944; Zeldovich, 1944) to explain the observations (Remington *et al.,* 2019; Hurricane *et al.,* 2014; Aglitsky *et al.,* 2010; Arnett *et al.,* 2009; Azechi *et al.,* 2007; Bell *et al.,* 2004a; 2004b; Class *et al.,* 2003a; 2003b; Peters, 2000; Sivashinsky, 1980; Williams, 1971; Buckmaster, 1979; Markstein, 1951; and Zeldovich, 1944). Some fundamental challenges remain.

First, the general theoretical framework (Landau, 1944) describes the evolution of a phase boundary and is relevant to phenomena far beyond the processes with gradually changing flow fields (Zhakhovsky *et al.,* 2014; Hurricane *et al.,* 2014; Anisimov *et al.,* 2013; Abarzhi, 2010; Kadau *et al.,* 2010; Arnett *et al.,* 2009; Meshkov, 2006; Class *et al.,* 2003a; and 2003b). We have yet to understand whether the interface is stable when the flow quantities experience sharp changes and the effects of dissipation and diffusion are negligible.

Second, direct diagnostics of various physical effects on the dynamics of a multi-phase flow requires detailed information of the interface structure. Such information is often a challenge to directly obtain in experiments and simulations (Remington *et al.,* 2019; Zhakhovsky *et al.,* 2019; Hurricane *et al.,* 2014; Azechi *et al.,* 2007; Bell *et al.,* 2004a; 2004b; Class *et al.,* 2003a; 2003b; Peters, 2000; Clavin, 2000; Rast, 1998; Wu *et al.,* 1997; and Williams, 1971). Besides, one usually observes the flow evolution from a far field and at time scales and length scales that are substantially greater than characteristic scales induced by, e.g., diffusion, dissipation, and finite interface thickness (Remington *et al.,* 2019; Zhakhovsky *et al.,* 2019; Abarzhi *et al.,* 2015; Kadau *et al.,* 2010; Landau and Lifshitz, 1987; and Sivashinsky, 1983). Thus, there is a need to directly connect the interface stability to the flow fields’ structure and to associate the mechanism of (de)stabilization of the interfacial dynamics with qualitative and quantitative properties of the vector and scalar fields of the flow (Abarzhi *et al.,* 2019).

Third, there are still formal questions to clarify in the classic theory of Landau (1944) for ideal incompressible fluids. For instance, Landau’s framework involves consideration of four independent variables: velocity potential of the heavy fluid, velocity potential of the light fluid, vortical potential of the light fluid, and interface perturbation (Landau and Lifshitz, 1987; Landau, 1944). These variables obey the principles of conservation of mass and normal and tangential components of momentum and a special condition for the perturbed mass flow at the interface. However, only two eigenvalues of the linearized dynamics are usually derived (Peters, 2010; Class *et al.,* 2003a; 2003b; Clavin, 2000; Landau and Lifshitz, 1987; Zeldovich, 1994; Williams, 1971; Landau, 1944). These eigenvalues are the unstable positive real eigenvalue (leading to LDI) and its stable negative counterpart. It is unclear whether the other eigenvalues exist and, for each of the fundamental solutions, what are the flow fields, how these flow fields link to interface stability, and what is the effect of fluctuations on the flow stability and fields.

In this work, we address these questions by systematically studying, in a far field approximation, the evolution of the discontinuous interface separating fluids of different densities and having inertial mass flux. We consider a two-dimensional flow of ideal incompressible fluids. By expanding the classic framework (Landau and Lifshitz, 1987), we implement and apply a general matrix method to solve the boundary value problem involving the governing equations in the fluid bulk and the boundary conditions at the interface and at the outside boundaries of the domain. We find formal fundamental solutions of the dynamics, analyze the interplay of interface stability with flow fields’ structure, and directly link rigorous mathematical attributes to physical observables (Abarzhi *et al.,* 2019; 2015). New mechanisms are of the interface stabilization and destabilization are identified. We find that the interfacial dynamics is stable when it conserves the fluxes of mass, momentum, and energy; the stabilization is due to inertial effects, leading to small oscillations of the interface velocity. This mechanism is absent in the classic Landau’s dynamics due to the postulated constancy of the interface velocity, causing LDI to develop. The interface destabilization is also associated with the imbalance of the perturbed energy at the interface, which is in agreement with the classic results (Landau, 1944). We identify the extreme sensitivity of the interface dynamics to the interfacial boundary conditions, including the formal properties of fundamental solutions and the qualitative and quantitative properties of the flow fields.

## II. METHOD

In this section, we outline the equations governing the interface dynamics, formulate the boundary value problem, and introduce a general matrix method to solve it.

### A. Governing equations

#### 1. Conservation laws and interfacial boundary conditions

The flow dynamics is governed by the conservation of mass, momentum, and energy. For ideal fluids, these are the Euler equations. In an inertial frame of reference,

where *x*_{i} are the spatial coordinates with (*x*_{1}, *x*_{2}, *x*_{3}) = (*x*, *y*, *z*), *t* is time, (ρ, **v**, *P*, *E*) are the fields of density ρ, velocity **v**, pressure *P*, and energy density *E* = ρ(*e* + **v**^{2}/2), with the specific internal energy *e* in its physics definition. The specific internal energy refers to the energy per unit mass that is contained in the system while excluding the kinetic energy of motion of the system as a whole and the potential energy of the system as a whole due to external force fields. The inertial frame of reference refers to a frame of reference moving with a constant velocity **V**_{0} = (0, 0, *V*_{0}). The Euler equations [Eqs. (1)] are augmented with the closure equation—the equation of state associating the internal energy and the pressure (Landau and Lifshitz, 1987).

We introduce a continuous local scalar function θ(*x*, *y*, *z*, *t*), whose derivatives $\theta \u0307$ and ∇θ exist (with dot marking a partial time-derivative), such that *θ* = 0 at the interface, and the heavy (light) fluid is located in the region *θ* > 0 (*θ* < 0) (Abarzhi *et al.,* 2015). In the bulk of the heavy (light) fluid, the flow fields are (ρ, **v**, *P*, *E*) → (ρ, **v**, *P*, *E*)_{h(l)}, and they are marked with subscript *h*(*l*). We represent the flow fields in the domain as (ρ, **v**, *P*, *E*) = (ρ, **v**, *P*, *E*)_{h}*H*(θ) + (ρ, **v**, *P*, *E*)_{l}*H*(−θ) by using the Heaviside step-function *H*(θ) (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018; and Abarzhi, 2010).

The balances of fluxes of mass and normal and tangential components of momentum and energy across the interface obey the boundary conditions (Landau and Lifshitz, 1987),

Here, [⋯] = 0 denotes the jump of functions across the interface, **n** and **τ** are the unit normal and tangential vectors at the interface with **n** = ∇θ/|∇θ| and (**n** · **τ**)=0, $j\u0303=\rho n\theta \u0307/\u2207\theta +v$ is the mass flux across the moving interface, and *W* is the specific enthalpy in its physics definition (Landau and Lifshitz, 1987). The value of specific enthalpy *W* includes the enthalpy of formation (Landau and Lifshitz, 1987).

Note that the interfacial boundary conditions [Eqs. (2)] are derived from the Euler equations [Eqs. (1)] independently of the equation of state and hold true for ideal fluids with any equation of state, both compressible and incompressible (Landau and Lifshitz, 1987).

Emphasize that in Eqs. (2), since $j\u0303\u22c5n=0$ and $j\u0303=\rho n\theta \u0307/\u2207\theta +v$, the condition of continuity of the tangential component of momentum at the interface $j\u0303\u22c5n\u2009j\u0303\u22c5\tau /\rho \u2009\tau =0$ is equivalent to the condition of continuity of the tangential component of velocity at the interface $v\u22c5\tau =0$ and the absence of shear at the interface, in agreement with the classic results (Landau and Lifshitz, 1987).

#### 2. Derivation of the interfacial boundary conditions

For readers’ convenience, we provide some details on how the interfacial boundary conditions are derived from the governing equations. As a sample case, we consider first of these conditions—the mass conservation condition (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018).

In the domain, the density of the fluid is represented as ρ = ρ_{h}*H*(θ) + ρ_{l}*H*(−θ), leading to

Here, we account for *∂H*(θ)/*∂*θ = *δ*(θ) and *δ*(θ) = *δ*(−θ), with *δ*(θ) being the Dirac delta-function. The other term is calculated as

We substitute these expressions in the equation for mass conservation, introduce the unit normal and tangential vectors at the interface **n** and **τ** as

and find

The expression is further reduced to the conditions in the fluids’ bulk and at the interface as in Eqs. (2),

The other interfacial boundary conditions in Eqs. (2) are derived likewise.

Note that the interfacial function θ is applied in numerical modeling of the multiphase flows, including the level set method and the in-phase field method (Osher and Fedkiw, 2002; Anderson *et al.,* 1998; and Sethian, 1996).

There is an important particular case in Eqs. (2), when the mass flux is conserved at the interface, $j\u0303\u22c5n=0$, and when the conserved mass flux is zero at the interface, $j\u0303\u22c5n\theta =0=0$. This special condition transforms the boundary conditions in Eqs. (2) to the continuity of the normal component of velocity [**v** · **n**] = 0, the continuity of the pressure [*P*] = 0, and the arbitrariness of the jumps of the tangential component of velocity [**v** · **τ**] = *arbitrary* and enthalpy[*W*] = *arbitrary* (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018). This case corresponds to the dynamics of contact discontinuity (in, e.g., immiscible fluids), such as Rayliegh–Taylor and Richtmyer–Meshkov instabilities (Meshkov and Abarzhi, 2019; Dell *et al.,* 2017; Abarzhi, 2010; Meshkov, 2006; Davies and Taylor, 1950; and Lord Rayleigh, 1883).

#### 3. Inertial reference frame, interface velocity, and outside boundaries

The governing equations and the boundary conditions [Eqs. (1) and (2)] are derived in an inertial frame of reference that moves with a constant velocity **V**_{0} (Fig. 1). Some discussion is required on the relation between the velocity of the inertial frame of reference **V**_{0} and the velocity $V\u0303$ of the interface as a whole. For a planar and steady interface, the interface velocity $V\u0303$ is constant ($V\u0303=V\u03030$), and it may be equal to the velocity **V**_{0} of the inertial frame of reference ($V\u03030=V0$). In the general case of a non-steady non-planar interface $V\u0303\u2260V\u03030$, the interface velocity $V\u0303$ and the velocity of the inertial frame of reference are distinct quantities ($V\u0303\u2260V\u03030$), and the interface velocity $V\u0303$ obeys the condition $V\u0303n=\u2212vn\theta =0=\u2212j\u0303/\rho n\theta =0$ (Peters, 2010; Landau and Lifshitz, 1987).

In our general theoretical framework, we consider the velocity **V**_{0} as the constant velocity of the inertial frame of reference. To study the effect of the interfacial boundary conditions on the interface velocity, we set that the velocity of the steady planar interface equals the velocity of the inertial frame of reference, as $V\u03030=V0$.

In the inertial frame of reference, the boundary conditions at the outside boundaries of the domain are

with constant *V*_{h(l)} (Fig. 1). The initial conditions are the initial flow fields at the interface and in the bulk. They define the characteristic length-scale and time-scale of the dynamics (Zeldovich, 1944).

#### 4. The flow dimensionality and the interfacial function

We consider a spatially extended two-dimensional flow propagating in the *z* direction, being periodic in the *x* direction, *x*+*λ*→*x*, and having no motion in the *y* direction (Fig. 1). We define the local scalar function θ as

### B. Linearized dynamics

#### 1. Interfacial boundary conditions to the leading and first orders

In Eqs. (1) and (2), the unperturbed flow fields are $j\u0303,v,P,W=J,V,P0,W0$, and the normal and tangential unit vectors of the unperturbed interface are {**n**, **τ**} = {**n**_{0}, **τ**_{0}}. We slightly perturb the flow fields as $j\u0303=J+j$, $v=V+u$, *P* = *P*_{0} + *p*, and *W* = *W*_{0} + $w$, with |**j**| ≪ |**J**|, |**u**| ≪ |**V**|, |*p*| ≪ |*P*_{0}|, and |$w$| ≪ |*W*_{0}|, respectively. We slightly perturb the fluid interface as **n** = **n**_{0} + **n**_{1} and **τ** = **τ**_{0} + **τ**_{1}, with |**n**_{1}| ≪ |**n**_{0}| and |**τ**_{1}| ≪ |**τ**_{0}|, and with $\theta \u0307/\u2207\theta \u226aV$. The fluid density is perturbed as ρ → ρ + *δ*ρ with |*δ*ρ| ≪ |ρ|. In the laboratory frame of reference, the perturbed velocity of the interface is $V\u0303=V\u03030+v\u0303$, with $v\u0303\u226aV\u03030$.

To the leading order, the boundary conditions at the interface are

To the first order, the boundary conditions at the interface are

Small perturbations of the flow fields decay away from the interface,

The boundary conditions in Eqs. (1) and (2) and (5a)–(5c) are valid for both compressible and incompressible ideal fluids with any equation of state, for two- and three-dimensional flows, and for arbitrary positioning of the interface relative to the mass flux (Landau and Lifshitz, 1987).

#### 2. Directionality, dimensionality, and density variations

The conditions [Eqs. (5a) and (5b)] can be simplified by applying the conditions of the directionality of the mass flux, the dimensionality of the flow, and the density variations in the fluids.

Indeed, to the leading order, the mass flux is **J** = ρ**V**, the flow fields are uniform in the bulk, (ρ, **v**, *P*, *W*)_{h(l)} = (ρ, **V**, *P*_{0}, *W*_{0})_{h(l)}, and obey conditions [Eqs. (3)] at the outside boundaries of the domain. The components of mass flux normal and tangential to the interface are *J*_{n} = **J** · **n**_{0} and *J*_{τ} = **J** · **τ**_{0}, respectively. We presume that the leading order mass flux is normal to the planar interface; hence, its tangential component is zero. This simplifies the interfacial boundary conditions to the leading order in Eqs. (5a) as

For a two-dimensional flow [Eqs. (4)], to the leading order, the normal and tangential vectors of the interface are **n**_{0} = (0, 0, −1) and **τ**_{0} = (1, 0, 0), respectively. The first order perturbations of the normal and tangential vectors of the interface are **n**_{1} = (*∂Z*^{*}/*∂x*, 0, 0) and **τ**_{1} = (0, 0, *∂Z*^{*}/*∂x*), respectively. This leads to **J** · **n**_{1} = 0 and further simplifies the interfacial boundary conditions to the first order in Eqs. (5b) as

For small density variations, |*δ*ρ/ρ| ≪ 2|(**j** · **n**_{0})/(**J** · **n**_{0})| and |*δ*ρ/ρ| ≪ 2|(**J** · **j**)/(**J**^{2})|, with |*δ*ρ/ρ| ≪ 2|**u**|/|**V**| and |*δ*ρ/ρ| ≪ |**u**|/|**V**|, the expressions in Eqs. (5e) are estimated as

For |*δ*ρ/ρ| ≪ |**u**|/|**V**|, the first order interfacial boundary conditions [Eqs. (5e)] are further simplified as

#### 3. Incompressible dynamics

In this paper, we consider the limiting case of incompressible dynamics in which the speed(s) of sound in the fluid(s) is substantially greater than other velocity scales, and the effect of density variations is negligible, (*δ*ρ/ρ) = 0. For ideal gases, the speed of sound is $c=\gamma P/\rho $, where *γ* is the fluid’s adiabatic index. For (*c*/*V*)_{h(l)} → ∞, the values approach $P0+Jn2/\rho hl\u2192P0hl$ and $W0+J2/2\rho 2hl\u2192W0hl$. Equations (5a) and (5d) are transformed to

By introducing temperature (*T*)_{h(l)} and the specific heat at constant pressure (*c*_{P})_{h(l)}, we present the physics enthalpy as $W0hl=W\xaf0+cPThl$. The enthalpy $W\xaf0$ is often used in engineering applications, and it has a jump at the interface, with $W\xaf0=Q\xaf$, where $Q\xaf$ is the heat release per unit mass. In reactive fluids, the value $Q\xaf$ is considered as the specific heat release of a chemical reaction at absolute zero temperature (Landau and Lifshitz, 1987).

Note the difference between the definitions of enthalpy *W*_{0} in physics sense accounting for the enthalpy of formation and the enthalpy $W\xaf0$ commonly applied in reactive fluids and engineering. Their difference is discussed in detail in Landau and Lifshitz (1987). In full consistency with this discussion, in the incompressible limit, the physics enthalpy *W*_{0} is continuous at the interface, [*W*_{0}] = 0, whereas the enthalpy $W\xaf0$ is discontinuous at the interface and its jump at the interface is $W\xaf0=Q\xaf$.

To the first order, the boundary conditions [Eqs. (5f)] at the interface take the form

where the first order perturbation of the mass flux is $j=\rho u+n0\theta \u0307$. Its normal component is $jn=j\u22c5n0$.

In addition to the boundary conditions at the interface, the governing equations in the bulk and at the outside boundaries of the domain to the first order are

In the laboratory frame of reference, the interface velocity is $V\u0303=V\u03030+v\u0303$, and the velocity perturbation is $v\u0303$,

### C. Fundamental solutions

#### 1. Solution structure

We seek solutions for the boundary value problem [Eqs. (7) and (8)] with the structure in which the velocity field of the heavy fluid is potential and the velocity field of the light fluid has both potential and vortical components (Landau and Lifshitz, 1987),

The rationale for this structure is the following: In the sub-domain of the heavy fluid, the flow is potential in accordance with the Kelvin theorem. In the sub-domain of the light fluid, the flow can be a superposition of the potential and vortical components. This structure is applied in the theoretical framework (Landau, 1944) and agrees with observations (Peters, 2000; Zeldovich, 1944). It is expected that the solution structure is established for any initial conditions (Landau, 1944; Zeldovich, 1944).

The scalar and vector potentials of the fluid and the interface perturbation are

Here, Ω is the growth-rate, i.e., the characteristic frequency or the eigenvalue of the system of linear differential equations, *k* = 2π/*λ* is the wavevector, and *λ* is the wavelength (spatial period). The pressure perturbations *p*_{h(l)} and the vortex length-scale $\lambda \u0303=2\pi /k\u0303$ are identified from Eqs. (8) as follows:

The vortical field does not contribute to pressure perturbations (Landau and Lifshitz, 1987). This leads to

In the system of Eqs. (7) and (8), the initial conditions set the characteristic length-scale 1/*k* and the time-scale 1/*kV*_{h}. We use dimensionless values for the growth-rate ω = Ω/*kV*_{h} and for the density ratio *R* = ρ_{h}/ρ_{l} with *R* ≥ 1. This leads to *V*_{l}/*V*_{h} = *R* and $k\u0303/k=\omega /R$. We use dimensionless values for the flow fields and the interface as φ = Φ/(*V*_{h}/*k*), $\phi \u0303=\Phi \u0303/Vh/k$, *ψ* = Ψ/(*V*_{h}/*k*), $z\xaf=kZ$ and for the variables as *kx* → *x*, *kz* → *z*, *kV*_{h}*t* → *t*. In dimensionless units, the fluid potentials are $\phi h=\phi eix+z+\omega \u2009t,\u2009\phi l=\phi \u0303eix\u2212z+\omega \u2009t$, and **ψ**_{l} = (0, ψ_{l}, 0) with $\psi l=\psi eix\u2212k\u0303/kz+\omega \u2009t$, the fluid velocities are **u**_{h} = ∇φ_{h}, **u**_{l} = ∇φ_{l} + ∇ × **ψ**_{l}, and the interface perturbation is $z*=z\xafeix+\omega \u2009t$.

#### 2. Eigenvalues and eigenvectors of fundamental solutions

With Eqs. (10), the system of differential equations governing the interface dynamics [Eqs. (7) and (8)] is reduced to the linear system M**r** = 0, where the matrix M is defined by the boundary conditions at the interface [Eq. (7)]. In dimensionless form, the elements of the matrix M are the functions of the growth-rate (the frequency, the eigenvalue) ω and the density ratio *R*, and the vector **r** is **r** = (φ_{h}, φ_{l}, *z*^{*}, ψ_{l})^{τ}. Traditionally, the linear system M**r** = 0 is solved by applying the compatibility condition: its eigenvalues are derived, and the (in)stability of the associated solutions is studied (Peters, 2010; Class *et al.,* 2003a; 2003b; Peters, 2000; Clavin, 2000; Landau and Lifshitz, 1987; Williams, 1971; Landau, 1944; Zeldovich, 1944). In this work, we proceed in a more formal way. Particularly, we develop and apply the general matrix method to analyze the interfacial dynamics [Eqs. (1)–(10)]. We find the fundamental solutions for the system of linear ordinary differential equations, including their eigenvalues and eigenvectors, and analyze, for each of these solutions, the stability and the structure of the flow fields. This formal analysis is required to better understand the interfacial dynamics. To our knowledge, it has not been done before.

We derive the eigenvalues ω_{i} by using the condition *det*Μ(ω_{i}, *R*) = 0 and identify the associated eigenvectors $e\u0303i$ by reducing the matrix Μ(ω_{i}, *R*) to row-echelon form. The matrix Μ is a 4 × 4 matrix. For a non-degenerate 4 × 4 matrix, the expected rank is 4, and the number of eigenvalues ω_{i} and associated eigenvectors $e\u0303i$ is 4, with *i* = 1,..., 4, corresponding to four degrees of freedom and four independent variables obeying four equations [Eqs. (7) and (10)].

Solution **r** for the system given by Eqs. (7) and (10) is a linear combination of fundamental solutions **r**_{i},

with *C*_{i} being the integration constants and $ri=ri\omega i,e\u0303i$ being the fundamental solutions with $ri=e\u0303ie\omega it$ and $e\u0303i=\phi eix+z,\u2009\phi \u0303eix\u2212z,\u2009z\xafeix,\u2009\psi eix\u2212k\u0303/k\u2009z\u2009i\u2009T$, and the associated vector $ei=\phi ,\u2009\phi \u0303,\u2009z\xaf,\u2009\psi \u2009i\u2009T$.

The linear system Μ**r** = 0 with Μ = Μ(ω, *R*) results from a linear system $Pr\u0307=Sr$, where P, S are 4 × 4 matrices, P = P(R), S = S(R), under the assumption that vector **r** varies in time as **r** ∼ *e*^{ωt}. This leads to Μ = (S − ωP). In a non-degenerate case, the inverse P^{−1} exists, and the system $Pr\u0307=Sr$ can be reduced to a standard form $r\u0307=P\u22121Sr$. The eigenvalues ω_{i} of the dynamics can then be found from the condition $detP\u22121S\u2212\omega I=0$, ω = {ω_{i}}, where the unit matrix is I, and index *i* marks the independent degrees of freedom. Equations *det*(P^{−1}S − ωI) = 0 and *det* M(ω, *R*) = 0 have the same eigenvalues ω_{i} because a linear combination of differential equations preserves the dynamic properties (Beklemishev, 1971).

By using the general matrix method for solving the boundary value problem [Eqs. (1)–(11)], we directly link the microscopic transport at the interface to macroscopic flow fields in the bulk, and we conduct a systematic study of the interplay of the interface stability with the structure of flow fields.

## III. RESULTS

In this section, we find the fundamental solutions for the boundary value problem in the conservative system and the classic and dynamic Landau’s systems. We identify the stability, the flow fields’ structure, the interface velocity, and the formal properties of these fundamental solutions.

### A. Conservative system

#### 1. Fundamental solutions

We consider the properties of dynamics balances of the fluxes of mass and normal and tangential components of momentum and energy across the interface [Eqs. (6)]. We call this dynamics the conservative dynamics (CD). For the conservative dynamics, the matrix Μ has the form Μ = *M*,

Its rank is 4. Its determinant is $detM=iR\u221212/R\omega \u2212R\omega +R\omega 2+R$, and the eigenvalues ω_{i} and eigenvectors **e**_{i} (derived as described before) are

#### 2. Stability of fundamental solutions

For the conservative system [Eqs. (12)], there are four fundamental solutions. Solutions $r12\omega 12,e12$ are stable, with $\omega 1=iR$, $\omega 2=\omega 1*$, and $Re\omega 12=0$. Solution **r**_{3}(ω_{3}, **e**_{3}) is unstable, with ω_{3} = *R* and Re[ω_{3}] > 0. Solution **r**_{4}(ω_{4}, **e**_{4}) is stable, with ω_{4} = −*R* and Re[ω_{4}] < 0. Figure 2 shows the dependence of real and imaginary parts of the eigenvalues on the density ratio.

#### 3. Flow fields of fundamental solutions

By substituting fundamental solutions **r**_{i} = **r**_{i}(ω_{i}, **e**_{i}) from Eqs. (12) to Eqs. (10), we find the flow fields for the conservative dynamics. For each of the fundamental solutions, Figs. 3–6 illustrate the interface perturbation *z*^{*}, the perturbed velocity fields **u**_{h(l)}, the perturbed velocity streamlines **s**_{h(l)} defined as $dshl/dt\xd7uhl=0$, the contour plots of the perturbed pressure *p*_{h(l)}, and (when applicable) the vortical field of the light fluid velocity, including the velocity vortical component ∇ × **ψ**_{l} and the contour plot of vorticity ∇ × **u**_{l} with $\u2207\xd7ul=0,k2\u2212k\u03032\psi l,0$. The real parts of the functions are shown at some density ratio and some instance of time in the motion plane. In contour plots, red (blue) marks the positive (negative) values. Each plot has its own range of values.

Fundamental solution **r**_{1}(ω_{1}, **e**_{1}) in Eqs. (12) describes stable oscillations of the flow fields in the vicinity of the interface, which decay away from the interface. For this solution, the velocity fields are potential, and all boundary conditions are satisfied, including balances of mass, momentum, and energy at the interface and the conditions at the outside boundaries of the domain. Note that for solution **r**_{1}(ω_{1}, **e**_{1}), the velocity field of the light fluid **u**_{l} has a phase shift π/2 relative to that of the heavy fluid **u**_{h}. That is, the field **u**_{l} is shifted relative to **u**_{h} by *λ*/4 along the *x* direction of translation [Eqs. (12)]. The velocity fields are also shifted relative to the interface perturbation. These shifts are needed to balance the perturbed mass flux across the perturbed interface in the incompressible approximation (Fig. 3).

For solution **r**_{1}, the perturbations of pressure *p*_{h(l)} oscillate in the vicinity of the interface and decay away from the interface. The pressure fields *p*_{h(l)} are shifted relative to the interface perturbation. The fields of pressure *p*_{h} and *p*_{l} are in anti-phase with one another, thus causing stable oscillations in the *z* direction. These fields span the same range of values with |*p*_{l}|_{max(min)} = |*p*_{h}|_{max(min)}. For solution **r**_{1}(ω_{1}, **e**_{1}), the vortical field is zero, ψ_{l} = 0, ∇ × **ψ**_{l} = 0, ∇ × **u**_{l} = 0 (Fig. 3).

The properties of the fundamental solution **r**_{2}(ω_{2}, **e**_{2}) in Eqs. (12) are similar to those of **r**_{1}(ω_{1}, **e**_{1}). This is because these solutions are complex conjugates, with $\omega 2=\omega 1*,e2=e1*,C2=C1*$ (Fig. 3).

Note that by their very meaning, the solutions $r12$ are traveling waves. The interference of these traveling waves results in the appearance of standing waves, whose flow fields experience stable oscillations, particularly solutions **r**_{CD} = (**r**_{1} + **r**_{2})/2 and $r\u0303CD=r1\u2212r2/2$. Figure 4 illustrates the flow fields for solution **r**_{CD}, including the perturbed velocity, the perturbed velocity streamlines, the pressure perturbations, and the interface perturbation at some density ratio and some instance of time in the motion plane.

Fundamental solution **r**_{3}(ω_{3}, **e**_{3}) obeys all the boundary conditions and is thus a physical solution. Because Re[ω_{3}(*R*)] > 0 for any density ratio *R* ≥ 1, it appears at first glance to describe unstable dynamics. However, a more detailed consideration finds that for this solution, with ω_{3} = *R* and **e**_{3} = (0, *i*, 0, 1), the velocity fields of the light and heavy fluids are identically zero, with **u**_{l} = 0 and **u**_{h} = 0. Furthermore, the interface perturbation and the pressure fields are also zero, *z*^{*} = 0 and *p*_{h(l)} = 0 (Fig. 5). Note that while, for this solution, the vortical velocity component is non-zero, ∇ × **ψ**_{l} ≠ 0, its vorticity is zero, ∇ × **u**_{l} = 0, because for the vorticity field with $\u2207\xd7ul=0,k2\u2212k\u03032\psi l,0$, the values are $k\u0303/k2=\omega 3/R2=1$ (Fig. 5). These properties hold true in the entire domain, at any time, and for any integration constant *C*_{3}. Therefore, solution **r**_{3} is unstable because its eigenvalue is ω_{3} > 0, and it corresponds to the unperturbed fields of the velocity, the pressure, and the interface because its eigenvector is **e**_{3} = (0, *i*, 0, 1) (Fig. 5).

Note that since the solution **r**_{3} has zero perturbed fields of the velocity, the vorticity, and the pressure, one may set, for this solution, the integration constant equal to zero. The solution **r**_{3} may also reflect the symmetry of the dynamics, associated with the arbitrariness of choice of the velocity of the inertial reference frame. Further investigations are required to clarify this issue.

Fundamental solution **r**_{4}(ω_{4}, **e**_{4}) appears at first glance to be stable and physical (Fig. 6). A more detailed consideration finds that the required choice of the integration constant *C*_{4} is *C*_{4} = 0 for this solution. Indeed, according to the expression for ω_{4}, **e**_{4}, at the outside boundaries of the domain, the velocity field of the heavy fluid decays (**u**_{h}|_{z→−∞} = 0) at any time. For the velocity field of the light fluid **u**_{l}, the potential component is vanishing away from the interfaces, ∇φ_{l}|_{z→+∞} = 0, whereas its vortical component ∇ × **ψ**_{l}, while decaying in time, as ∼*e*^{−Rt}, increases away from the interface, as ∼*e*^{z}. Note that while, for this solution, the vortical component of velocity is non-zero, ∇ × **ψ**_{l} ≠ 0, its vorticity is zero, ∇ × **u**_{l} = 0, because for the vorticity field with $\u2207\u2009\xd7\u2009ul=0,k2\u2212k\u03032\psi l,0$ the value is $k\u0303/k2=1$. Finally, note that, for this solution, the pressure fields *p*_{h(l)} are anti-phase with one another and decay away from the interface. Thus, for solution **r**_{4} to satisfy boundary condition $ulz\u2192+\u221e=0$, we must have *C*_{4} = 0 (Fig. 6). Note that due to the exponential growth of the vortical component of the velocity ∇ × **ψ**_{l} away from the interface, ∼$ez$, the integration constant must remain zero *C*_{4} = 0 even when the boundary conditions away from the interface are slightly perturbed (by, e.g., some external noise).

#### 4. Interface velocity for fundamental solutions

For fundamental solutions **r**_{1(2)}, which describe stable traveling waves, as well as for solutions $rCD,r\u0303CD$, which describe stable standing waves resulting from interference of the traveling waves **r**_{1(2)}, the interface velocity is $V\u0303=V\u03030+v\u0303$ in the laboratory frame of reference with $v\u0303n0=\u2212uhn0+\theta \u0307\theta =0$. Because of the oscillatory terms $uhn0+\theta \u0307\u223ce\xb1iRt$, the interface velocity experiences stable oscillations near the steady value $V\u03030$. The amplitude of the oscillations is small, $v\u0303\u226aV\u03030$, and is set by the initial conditions and the integration constant(s) *C*_{1(2)}, whereas the period of the oscillations $1/R$ is substantially smaller than the characteristic time-scale of the dynamics.

For fundamental solution **r**_{3} in the laboratory frame of reference, the interface velocity is constant, $V\u0303=V\u03030$, because its perturbed velocity fields and perturbed interface are zero for any integration constant *C*_{3}. For solution **r**_{4}, in the laboratory frame of reference, the interface velocity is constant $V\u0303=V\u03030$ because its integration constant is zero, *C*_{4} = 0.

#### 5. Formal properties of the conservative system

To study the formal properties of the conservative system, we find, for matrix M = *M*, the associated matrices S = *S*_{M} and P = *P*_{M},

The solutions of equations *det*(P_{M}^{−1}S_{M} − ωI) = 0 and *detM* = 0 yield the same eigenvalues: $\omega 12=\xb1iR$ and ω_{3} = *R*, ω_{4} = −*R*. The conservative dynamics has four eigenvalues, four fundamental solutions, and four independent degrees of freedom. It is non-degenerate.

#### 6. Summary of properties of solutions for the conservative system

For the conservative system, there are four fundamental solutions **r**_{i}(ω_{i}, **e**_{i}), *i* = 1, 2, 3, 4 [Eqs. (12), Figs. 2–6, and Table I]. Solutions **r**_{1} and **r**_{2} with integration constant $C1,C2=C1*$ describe the stable oscillations of the flow fields. These solutions are stable traveling waves; their interference results in stably oscillating standing waves **r**_{CD} and $r\u0303CD$. For these solutions, the interface velocity $V\u0303$ experiences small stable oscillations near the steady value $V\u03030$.

. | r_{1}
. | r_{2}
. | r_{3}
. | r_{4}
. |
---|---|---|---|---|

Re[ω] | 0 | 0 | >0 | <0 |

Im[ω] | >0 | <0 | 0 | 0 |

C | * | * | *(0) | 0 |

u_{h} | * | * | 0 | 0 |

u_{l} | * | * | 0 | 0 |

∇×u_{l} | 0 | 0 | 0 | 0 |

z^{*} | * | * | 0 | 0 |

. | r_{1}
. | r_{2}
. | r_{3}
. | r_{4}
. |
---|---|---|---|---|

Re[ω] | 0 | 0 | >0 | <0 |

Im[ω] | >0 | <0 | 0 | 0 |

C | * | * | *(0) | 0 |

u_{h} | * | * | 0 | 0 |

u_{l} | * | * | 0 | 0 |

∇×u_{l} | 0 | 0 | 0 | 0 |

z^{*} | * | * | 0 | 0 |

Solution **r**_{3}, corresponds to the unperturbed fields of the velocity, the pressure, and the interface at any constant *C*_{3}. For solution **r**_{4}, the value *C*_{4} must be zero to obey all the boundary conditions, as discussed in the foregoing, *C*_{4} = 0. For solutions **r**_{3(4)}, the interface velocity is constant, $V\u0303=V\u03030$.

### B. Classic Landau’s system

#### 1. Fundamental solutions

Landau (1944) analyzed the stability of the interface by considering the boundary conditions balancing the transport of mass and normal and tangential components of momentum at the interface and by omitting the equation for the energy balance at the interface (Peters, 2010; Landau, 1944; Zeldovich, 1944; and Mallard and Le Chatelier, 1883). Since four degrees of freedom and four independent variables require four conditions, Landau employed an additional condition for the continuity of the normal component of the perturbed velocity at the interface $u\u22c5n0=0$ or $jn\theta =0=0$ at the interface (Landau, 1944).

With this condition, the boundary conditions have the form

Recall that the perturbed mass flux is $j=\rho u+n0\theta \u0307$ with $jn=j\u22c5n0=\rho u\u22c5n0+\theta \u0307$. The condition *j*_{n}|_{θ=0} = 0 leads to $u\u22c5n0=\u2212\theta \u0307$ for any ρ and thus to $u\u22c5n0=0$ in Eqs. (13). The condition for the normal component of momentum is transformed to [*p***n**_{0}] = 0, and the field of pressure perturbation is continuous at the interface (Landau, 1944).

For the classic Landau’s system, the matrix Μ is Μ = *L*,

Its rank is 4. Its determinant is $detL=iR\u22121/R\omega \u2212RR+1\omega 2+2R\omega \u2212RR\u22121$. Its eigenvalues ω_{i} and eigenvectors **e**_{i} are

#### 2. Stability of fundamental solutions

For the classic Landau’s dynamics, there are three fundamental solutions. Solution **r**_{1}(ω_{1}, **e**_{1}) is unstable, Re[ω_{1}] > 0. Solution **r**_{2}(ω_{2}, **e**_{2}) is stable, Re[ω_{2}] < 0. Solution **r**_{3}(ω_{3}, **e**_{3}) is unstable, Re[ω_{3}] > 0, and is identical to that for the conservative system. Figure 7 illustrates the dependence of eigenvalues on the density ratio.

#### 3. Flow fields of fundamental solutions

By substituting solutions **r**_{i} = **r**_{i}(ω_{i}, **e**_{i}) in Eqs. (14) to Eqs. (10) for each *i*, we find the flow fields for fundamental solutions describing the classic Landau’s dynamics.

For each of fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) for the classic Landau’s dynamics, Figs. 8 and 9 illustrate the interface perturbation *z*^{*}, the perturbed velocity fields **u**_{h(l)}, the perturbed velocity streamlines **s**_{h(l)} defined as $dshl/dt\xd7uhl=0$, the contour plots of the perturbed pressure fields *p*_{h(l)}, and the vortical field of the light fluid velocity, including the vortical velocity component ∇ × **ψ**_{l} and the contour plot of vorticity ∇ × **u**_{l} with $\u2207\xd7ul=0,k2\u2212k\u03032\psi l,0$. The real parts of the functions are shown at some density ratio and some instance of time in the motion plane. In contour plots, red (blue) marks the positive (negative) values. Each plot has its own range of values.

Fundamental solution **r**_{1}(ω_{1}, **e**_{1})in Eqs. (14) corresponds to the Landau–Darrieus instability (Landau, 1944) (Fig. 8). This solution, although unstable, is a physical solution. It satisfies all boundary conditions, including the assigned boundary conditions at the interface and the conditions at the outside boundaries of the domain [Eqs. (3) and (13)]. This solution describes dynamics in which the interface perturbation and the vortical and potential components are strongly coupled. The perturbed velocity fields **u**_{h(l)} and the perturbed velocity streamlines **s**_{h(l)} illustrate the formation of vortical structures near the interface and in the bulk of the light fluid. For solution **r**_{1}(ω_{1}, **e**_{1}), the vortical component, ∇ × **ψ**_{l}, and the vorticity ∇ × **u**_{l} of the velocity of the light fluid, while increasing in time, decay away from the interface. The fields of pressure perturbations *p*_{h} and *p*_{l} are equal to one another at the interface and decay away from the interface. They are in phase with one another and in anti-phase with the interface perturbation. The fields of pressure *p*_{h(l)} are symmetric and span the same range of values with |*p*_{l}|_{max(min)} = |*p*_{h}|_{max(min)} (Fig. 8).

The vortical field for solution **r**_{1}(ω_{1}, **e**_{1}) in Eqs. (14) has the following properties (Fig. 8): (1) Consistent with expressions in Eqs. (14), the potential components of **u**_{h} and **u**_{l} are in anti-phase with one another, whereas the vortical component of velocity **u**_{l} is shifted relative to the potential components of velocities **u**_{h} and **u**_{l} and the interface perturbation *z*^{*}. (2) For the length-scale of the vortical field, we have $k\u0303/k=\omega 1/R$, $k\u0303/k=\u2212R+\u2212R+R2+R3/R1+R$. The length-scale of the vortical field is large, with $k\u0303/k\u223cR\u22121/2\u21920$ for *R* → 1^{+} and with $k\u0303/k\u223cR\u22121/2\u21920$ for *R* → ∞. The maximum value of $k\u0303/k$ (minimum value of $\lambda \u0303/\lambda $) is achieved at $R=2+5\u22484.24$. (3) For *R* → 1^{+}, the values are $\psi /\phi \u0303\u223cR\u22121$ and $\psi /\phi \u0303\u223cR\u22121$, whereas for *R* → ∞, the values are $\psi /\phi \u0303\u223cR3/2$ and $\psi /\phi \u0303\u223c1$. Therefore, the contribution of the vortical field to the dynamics is quantitatively large for fluids with very different densities and is quantitatively small for fluids with similar densities. (4) To quantify the strength of the vorticity field for given ψ, we note that $\u2207\xd7ul/\psi k2=1\u2212k\u0303/k2$. For solution **r**_{1}(ω_{1}, **e**_{1}), this value is ∼1 for all *R* ≥ 1.

The other fundamental solution in Eqs. (14) is solution **r**_{2}(ω_{2}, **e**_{2}) (Fig. 9). Solution **r**_{2} describes the dynamics for which the interface perturbation and the vortical and potential components are strongly coupled. It appears stable and physical at first glance. Its fields of pressure *p*_{h(l)} are equal to one another and are in-phase at the interface and decay away from the interface. For this solution, the field of the perturbed velocity of the heavy fluid vanishes at the outside boundaries of the domain **u**_{h}|_{z→−∞} = 0 at all times. The potential component of the perturbed velocity of the light fluid **u**_{l} also vanishes away from the interface, ∇φ_{l}|_{z→+∞} = 0. Its vortical component, ∇ × **ψ**_{l}, decays exponentially in time but increases far from the interface. The strength of the vorticity field ∇ × **u**_{l} for this solution also increases far from the interface. For fundamental solution **r**_{2}(ω_{2}, **e**_{2}) in Eqs. (14) to satisfy the boundary condition **u**_{l}|_{z→+∞} = 0 at any time, the value of integration constant *C*_{2} should be zero (*C*_{2} = 0). This integration constant must remain zero *C*_{3} = 0 even when the boundary conditions away from the interface are slightly perturbed (due to the exponential growth of the vortical component of the velocity).

Fundamental solution **r**_{3}(ω_{3}, **e**_{3}) in Eqs. (14) is identical to that in Eqs. (12) (Fig. 5). Note that solution **r**_{3} exists in the classic framework (Landau and Lifshitz, 1987; Landau, 1944) and can be found similarly to Eqs. (13) and (14).

#### 4. Interface velocity for fundamental solutions

For fundamental solution **r**_{1}, in the laboratory frame of reference, the interface velocity is $V\u0303=V\u03030+v\u0303$ with $v\u0303n0=\u2212uhn0+\theta \u0307\theta =0$. While in this expression, each of the terms **u**_{h} and $\theta \u0307$ grows exponentially in time, $uh,\theta \u0307\u223ce\omega 1t$, the exact balance in the boundary conditions $u\u22c5n0=0$ in Eqs. (13) leads to $uhn0+\theta \u0307\theta =0=0$ and $v\u0303=0$ and thus to the constancy of the interface velocity $V\u0303=V\u03030$. The classic Landau’s dynamics is the perfect mathematical match.

For fundamental solutions $r23$, in the laboratory frame of reference, the interface velocity is constant $V\u0303=V\u03030$, similarly to the foregoing.

#### 5. Formal properties of the classic Landau’s system

To study the formal properties of the classic Landau’s system, we find, for matrix Μ = *L*, the associated matrices S = *S*_{L} and P = *P*_{L},

In matrix *P*_{L}, the fourth row only has null elements and *detP*_{L} = 0. Thus, the inverse matrix $PL\u22121$ does not exist. This suggests the degeneracy of the classical Landau’s system and a singular character of the Landau–Darrieus instability.

#### 6. Summary of properties of solutions for the classic Landau’s system

For the classic Landau’s system, there are three fundamental solutions **r**_{i}(ω_{i}, **e**_{i}), *i* = 1, 2, 3 [Eqs. (13) and (14), Figs. 5 and 7–9, and Table II]. Solution **r**_{1}(ω_{1}, **e**_{1}) with *C*_{1} ≠ 0 describes the Landau–Darrieus instability. For solution **r**_{2}(ω_{2}, **e**_{2}), the integration constant must be zero, *C*_{2} = 0. Solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed fields of velocity and pressure at any *C*_{3}. For each of these solutions, in the laboratory frame of reference, the interface velocity is constant $V\u0303=V\u03030$.

### C. Dynamic Landau’s system

#### 1. Fundamental solutions

The paradox of the classic Landau’s system (which to our knowledge has not been discussed in detail before) is that it has a smaller than expected number of fundamental solutions. In physics sense, this paradox can be understood from the following considerations: The boundary conditions in the classic Landau’s system [Eqs. (13)] imply dynamic conditions for the conservation of mass and normal and tangential components of momentum at the interface because they contain partial time-derivatives of components of vector **r**. Landau’s special condition in Eqs. (13) is a “static” condition of continuity of the normal component of the perturbed velocity at the interface [**u** · **n**_{0}] = 0. This condition implies that the perturbed mass flux is a constant value at the interface, $jn|\theta =0$ = *const.*, and furthermore, this value is zero, $jn|\theta =0$ = 0. In order to resolve the paradox, we modify the special Landau’s condition $jn|\theta =0$ = 0 to a more general form $jn|\theta =0$ = *const.* and implement it in a dynamic condition, $\u2202u\u22c5n0/\u2202t=0$, at the interface. While other implementations are possible, we note that the condition of constant mass flux at the interface is met in realistic environments (Azechi *et al.,* 2007; Peters, 2000).

For this dynamic implementation of Landau’s system, the modified boundary conditions are

For the dynamic Landau’s system, the matrix Μ has the form $M=L\u0303$,

Its rank is 4. Its determinant is $detL\u0303=iR\u22121/R\omega \omega \u2212RR+1\omega 2+2R\omega \u2212RR\u22121$. It has four eigenvalues ω_{i} and four eigenvectors **e**_{i}. For *i* = 1, 2, 3, the eigenvalues ω_{i} and eigenvectors **e**_{i} are the same as in Eqs. (14). For *i* = 4, the fundamental solution is $r4o\omega 4o,e4o$ with

Symbol “o” in $r4o\omega 4o,e4o$ emphasizes the presence of the null eigenvalue.

#### 2. Stability of fundamental solution

#### 3. Flow fields of fundamental solutions

Fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 3 in Eqs. (16) are the same as solutions for system [Eqs. (14)] and have the same corresponding fields of the velocity, the pressure, and the interface perturbation (Figs. 5, 8, and 9). For fundamental solution $r4o\omega 4o,e4o$, Fig. 10 presents the fields of the perturbed velocity **u**_{h(l)}, the perturbed streamlines **s**_{h(l)}, the perturbed pressure *p*_{h(l)}, the perturbed vortical field ∇ × **ψ**_{l}, the vorticity $\u2207\xd7ul=0,k2\u2212k\u03032\psi l,0$, and the interface perturbation $z*$ at some instance of time for some density ratio in the motion plane. Each plot has its own range of values.

Solution $r4o\omega 4o,e4o$ is independent from fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 3 in Eqs. (14) and has the following properties: (1) The vortical component of velocity **u**_{l} is shifted relative to the potential components of velocities **u**_{h}, **u**_{l} and the interface perturbation *z*^{*}, whereas the potential components of **u**_{h} and **u**_{l} are in anti-phase with one another. (2) The vortical field length-scale is infinitely large, with $k\u0303/k=0$ and $\lambda \u0303/\lambda =\u221e$, due to the null eigenvalue $\omega 4o=0$. (3) The pressure fields of *p*_{h} and *p*_{l} are in phase and decay away from the interface. The pressure fields are no longer symmetric and |*p*_{l}|_{max(min)} ≪ |*p*_{h}|_{max(min)}. (4) The vortical component of the velocity is ∇ × **ψ**_{l} ≠ 0, and the vorticity is ∇ × **u**_{l} ≠ 0. They change periodically in the *x* direction and are uniform in the *z* direction. (5) The values are |ψ/φ| ∼ *R* and $\psi /\phi \u0303\u223c2$ for *R* → ∞ and $\psi /\phi \u223c\psi /\phi \u0303\u223cR\u22121$ for *R* → 1^{+}, and |∇ × **u**_{l}|/|ψ|*k*^{2} = 1 for any *R*.

For solution $r4o\omega 4o,e4o$, the velocity **u**_{h} vanishes at *z* → −∞, whereas for the velocity **u**_{l} at *z* → +∞, the potential component decays and the vortical component changes periodically. The boundary conditions **u**_{l}|_{z→+∞} = 0 require us to set, for this solution, the integration constant equal to zero. Note that for solution $r4o\omega 4o,e4o$, some slight perturbations of the boundary conditions away from the interface *z* → +∞ (due to some external noise) may lead to non-zero integration constant *C*_{4}.

#### 4. Interface velocity for fundamental solutions

For fundamental solutions **r**_{i} with *i* = 1, 2, 3 in Eqs. (16), the interface velocity in the laboratory frame of reference is constant, $V\u0303=V\u03030$, as discussed in the foregoing.

For fundamental solutions $r4o$ in Eqs. (16) in the laboratory frame of reference, the interface velocity is $V\u0303=V\u03030+v\u0303$ with $v\u0303n0=\u2212uhn0+\theta \u0307\theta =0$. Due to the zero eigenvalue, $\omega 4o=0$, the value is $\theta \u0307=0$, whereas the value is $uhn0\theta =0=const$. Thus, for solution $r4o$ with a constant perturbed mass flux across the interface, $\u2202u\u22c5n0/\u2202t=0$, for non-zero integration constant *C*_{4}, the interface velocity $V\u0303=V\u03030+v\u0303$ is slightly shifted from its steady value $V\u03030$, $v\u0303\u226aV\u03030$, and the shift is stationary, $v\u0303\u0307=0$. For zero integration constant *C*_{4} = 0, the interface velocity is constant.

#### 5. Formal properties of the dynamic Landau’s system

##### a. Non-degeneracy of the dynamic Landau’s system.

To study the formal properties of the dynamic Landau’s system, we find for matrix $M=L\u0303$ the associated matrices $S=SL\u0303$ and $P=PL\u0303$,

The equations $detPL\u0303\u22121SL\u0303\u2212\omega I=0$ and $detL\u0303=0$ yield the same eigenvalues. This implies that the dynamic Landau’s system is non-degenerate. It has four fundamental solutions for four degrees of freedom, thus resolving the long-standing formal paradox.

##### b. Singular character of the classic Landau’s system.

An important consequence of the zero eigenvalue $\omega 4o=0$ is the singular nature of the LDI and the classic Landau’s solution (Kadanoff *et al.,* 1967; Barenblatt *et al.,* 1962). Indeed, differential equations [Eqs. (13)–(16)] govern the dynamics of slight perturbations **r** near the equilibrium state of the uniform flow fields in Eqs. (6) and (7)–(10) for which **r** = 0. It is presumed that these deviations are exponential in time, with **r** varying as **r** ∼ *e*^{ωt}, and with characteristic time-scale being ∼|ω^{−1}|. For ω = 0, formal substitution leads to **r** ∼ *e*^{0·t} ∼ *const*. This implies that slight deviations from the equilibrium dynamics are actually power-law functions of time, **r** ∼ *t*^{a}, rather than exponential functions, **r** ∼ *e*^{ωt} (Kadanoff *et al.,* 1967; Barenblatt *et al.,* 1962). The requirement of the asymptotic stability of fundamental solution $r4o\omega 4o,e4o$ leads to condition Re[*a*] < 0 (Landau and Lifshitz, 1987). For Re[*a*] < 0, slight initial deviations from the equilibrium, **r** ∼ *t*^{a}, vanish asymptotically with time, **r** → 0 for *t* → ∞. To derive the exponent *a*, one has to account for nonlinear terms. We address this problem to future research.

Note that the power-law functions are scale-invariant and describe processes with no characteristic scales. The power-law dependence **r** ∼ *t*^{a} may also imply that the system may have a singularity and experience a phase transition (Kadanoff, 2000; Kadanoff *et al.,* 1967). Note also that with $\omega 4o=0$, solution $r4o\omega 4o,e4o$ in Eqs. (16) can be viewed as a steady vortical field that is imposed at some instance of time (by some external noise) and is neutrally stable (Figs. 8 and 10). This vortical field can “seed” the vortical field in Landau’s solution **r**_{1}(ω_{1}, **e**_{1}) in Eqs. (14) and trigger the unstable dynamics.

#### 6. Summary of properties of fundamental solutions for the dynamic Landau’s system

For the dynamic Landau’s system in Eqs. (15) and (16), there are four fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, ..., 4 [Eqs. (14) and (16), Table II, and Figs. 5 and 7–10]. Solution **r**_{1}(ω_{1}, **e**_{1}) with *C*_{1} ≠ 0 describes the Landau–Darrieus instability. For solution **r**_{2}(ω_{2}, **e**_{2}), the integration constant is *C*_{2} = 0. Solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed fields of the velocity, the pressure, and the interface at any *C*_{3}. Solution $r4o\omega 4o,e4o$ is neutrally stable and may lead to a power-law dynamics near the equilibrium seeding LDI.

For solutions **r**_{i} with *i* = 1, 2, 3, the interface velocity in the laboratory frame of reference is constant, $V\u0303=V\u03030$, as discussed in the foregoing. For solution $r4o$with *C*_{4} ≠ 0, the interface velocity $V\u0303=V\u03030+v\u0303$ is slightly shifted from its steady value $V\u03030$, $v\u0303\u226aV\u03030$, and the shift is stationary, $v\u0303\u0307=0$. For zero integration constant *C*_{4} = 0, the interface velocity is constant.

## IV. PROPERTIES OF INTERFACIAL DYNAMICS IN THE CONSERVATIVE SYSTEM AND IN THE LANDAU’S SYSTEMS

In this section, we focus on the physical properties of the interfacial dynamics that are prescribed by fundamental solutions for the conservative system and the classic and dynamic Landau’s systems. We further find the new mechanisms of the interface stabilization due to the inertial effects, and we model the influence of energy fluctuations on the interface dynamics.

### A. Physical properties of fundamental solutions

#### 1. Stable, neutrally stable, and unstable standing waves

We compare the properties of fundamental solution **r**_{CD} for the conservative system, fundamental solution **r**_{LD} for the classic Landau’s system, and fundamental solution **r**_{DLD} for the dynamic Landau’s system, where **r**_{CD} = (**r**_{1} + **r**_{2})/2 in Eqs. (12), **r**_{LD} = **r**_{1} in Eqs. (14), and $rDLD=r4o$ in Eqs. (16). Solutions **r**_{CD}, **r**_{LD}, and **r**_{DLD} are standing waves; they have distinct properties (Tables III–V).

r_{CD} | $\omega CD=\xb1iR$ |

r_{LD} | $\omega LD=\u2212R+R3+R2\u2212RR+1$ |

r_{DLD} | $\omega DLD=0$ |

r_{CD} | $\omega CD=\xb1iR$ |

r_{LD} | $\omega LD=\u2212R+R3+R2\u2212RR+1$ |

r_{DLD} | $\omega DLD=0$ |

. | CD . | LD . |
---|---|---|

Conservation properties | Conserves mass, momentum, | Conserves mass and momentum |

and energy at the interface | and has zero perturbed mass flux at the interface | |

Interface velocity | Slight stable oscillations near | Constant value (postulate) |

the constant value | ||

Flow field | Potential velocity fields | Vortical structures in the light fluid |

Interfacial shear | Shear-free | Shear-free |

Formal properties | Non-degenerate; four fundamental | Degenerate; three fundamental |

solutions and four | solutions and four degrees of | |

degrees of freedom | freedom | |

Stability | Stable; stabilized by inertial | Unstable |

effects |

. | CD . | LD . |
---|---|---|

Conservation properties | Conserves mass, momentum, | Conserves mass and momentum |

and energy at the interface | and has zero perturbed mass flux at the interface | |

Interface velocity | Slight stable oscillations near | Constant value (postulate) |

the constant value | ||

Flow field | Potential velocity fields | Vortical structures in the light fluid |

Interfacial shear | Shear-free | Shear-free |

Formal properties | Non-degenerate; four fundamental | Degenerate; three fundamental |

solutions and four | solutions and four degrees of | |

degrees of freedom | freedom | |

Stability | Stable; stabilized by inertial | Unstable |

effects |

. | LD . | DLD . |
---|---|---|

Conservation properties | Conserves mass and momentum | Conserves mass and momentum |

and has the null perturbed | and has a constant perturbed | |

interfacial mass flux | mass flux at the interface | |

Interface velocity | Constant value (postulate) | May be slightly shifted near the |

constant value | ||

Flow field | The vortical field length-scale | The vortical field length-scale is |

depends on the density ratio and | independent of the density ratio | |

is large and finite when compared | and is infinite when compared to | |

to the perturbation wavelength | the perturbation wavelength | |

Interfacial shear | Shear-free | Shear-free |

Formal properties | Degenerate; three fundamental | Non-degenerate; four fundamental |

solutions and four degrees of | solutions and four degrees | |

freedom | of freedom | |

Stability | Unstable | Neutrally stable |

. | LD . | DLD . |
---|---|---|

Conservation properties | Conserves mass and momentum | Conserves mass and momentum |

and has the null perturbed | and has a constant perturbed | |

interfacial mass flux | mass flux at the interface | |

Interface velocity | Constant value (postulate) | May be slightly shifted near the |

constant value | ||

Flow field | The vortical field length-scale | The vortical field length-scale is |

depends on the density ratio and | independent of the density ratio | |

is large and finite when compared | and is infinite when compared to | |

to the perturbation wavelength | the perturbation wavelength | |

Interfacial shear | Shear-free | Shear-free |

Formal properties | Degenerate; three fundamental | Non-degenerate; four fundamental |

solutions and four degrees of | solutions and four degrees | |

freedom | of freedom | |

Stability | Unstable | Neutrally stable |

The dynamics **r**_{CD} is stable, the dynamics **r**_{LD} is unstable, and the dynamics **r**_{DLD} is neutrally stable. For solution **r**_{CD}, the velocity fields are potential; for solutions **r**_{LD} and **r**_{DLD}, the potential and vortical components of the velocities and the interface perturbation are strongly coupled. There is zero shear at the interface for each of the dynamics **r**_{CD}, **r**_{LD}, and **r**_{DLD}. Solution **r**_{CD} conserves mass, momentum, and energy at the interface. Solution **r**_{LDG} conserves mass and momentum and has zero perturbed mass flux at the interface. Solution **r**_{DLD} conserves mass and momentum and has a constant perturbed mass flux at the interface. The dynamics **r**_{CD} is non-degenerate (four solutions and four degrees of freedom). The dynamics **r**_{LD} is degenerate (three solutions and four degrees of freedom). The dynamics **r**_{DLD} is non-degenerate (four solutions and four degrees of freedom) (Tables III–V).

For solutions **r**_{DLD} and **r**_{LD}, the potential and vortical components of the velocities are strongly coupled with the interface perturbations. Yet, for solution **r**_{LD}, the length-scale of the vortical field $\lambda \u0303=2\pi /k\u0303$ depends on the density ratio, whereas for solution **r**_{DLD}, the length-scale of the vortical field is independent of the density ratio and is infinitely large ($\lambda \u0303=\u221e$) (Tables III–V).

#### 2. Other unstable solutions

The conservative system, the classic Landau system, and the dynamic Landau system have a fundamental solution in common, **r**_{3}(ω_{3}, **e**_{3}) in Eqs. (12), (14), and (16), which can be implemented for any integration constant *C*_{3}. This solution is the fastest unstable solution, as suggested by its eigenvalue ω_{3}. This solution has the null perturbation fields of the velocity, pressure, and interface, as suggested by its eigenvector **e**_{3} (Fig. 5).

Note that this solution is in agreement with the classic works (Landau and Lifshitz, 1987) and can be derived from the governing equations in the framework (Landau, 1944). The existence of this solution may imply that the velocity of the inertial frame of reference **V**_{0} is a free quantity, and it may differ from the velocity of the planar steady interface $V\u03030$. It may also imply that the interfacial dynamics may have flow fields that are uniform and are indistinguishable from those of the equilibrium state. While the solution is formally unstable, its dynamics may appear as being “fully stabilized.” This apparent stabilization is achieved in the case of ideal incompressible fluids and at any density ratio. In realistic fluids, various physical effects may influence the null perturbation fields of the velocity, pressure, and interface associated with this solution, to be studied in the future research.

#### 3. Other stable solutions

Our analysis identifies some new properties of stable fundamental solutions, which are the stable solutions that decay exponentially with time due to their real negative eigenvalues. These include solution **r**_{4} in Eqs. (12) for the conservative system and solution **r**_{2} in Eqs. (14) and (16) for the classic and dynamic Landau’s systems. As discussed in the foregoing, these stable solutions require a null integration constant in order to obey the boundary conditions at the outside boundaries of the domain at any times.

A somewhat similar requirement was discussed in the study of the LDI by means of the Laplace transform (Zeldovich, 1944). The study interpreted a non-zero vortical component of the velocity field away from the interface as vortical structures that are produced at the interface and are moved to the bulk. The study required a null integration constant for the corresponding solution because these vortical structures have no influence on the development of LDI (Zeldovich, 1944). In our analysis, we require the null integration constant for the exponentially decaying stable solutions to obey the boundary conditions at the outside boundaries of the domain at any time for any density ratio. We further note that for stable solutions decaying exponentially with time, the requirement of zero integration constant directly follows from the solution structure. Thus, it should hold in realistic fluids, especially in the case of fluids with contrasting densities, when the contribution of the vortical field to the dynamics is large.

### B. Mechanisms of the interface stabilization and destabilization

#### 1. Interface stabilization for the conservative dynamics

To better understand the mechanism of the interface stabilization of the conservative dynamics, we consider the interface velocity. In the laboratory reference frame, the interface velocity is $V\u0303=V\u03030+v\u0303$, with $v\u0303\u22c5n0=\u2212u\u22c5n0+\theta \u0307\theta =0$.

For the conservative dynamics **r**_{CD}, the values are $u\u22c5n0\theta =0\u223ce\xb1iRt,\u2009\theta \u0307\theta =0\u223ce\xb1iRt$, leading to $u\u22c5n0+\theta \u0307\theta =0\u223ce\xb1iRt$ and $v\u0303\u22c5n0\u223ce\xb1iRt$. Thus, the interface velocity experiences small stable oscillations near the steady value $V\u03030,\u2009V\u0303\u2212V\u03030\u22c5n0\u223ce\xb1iR\u2009t$, as discussed in the foregoing.

This suggests the inertial effect as the stabilization mechanism of the conservative dynamics **r**_{CD}. When the interface is slightly perturbed, the parcels of the heavy fluid and the light fluid follow the interface perturbation causing the change of momentum and energy of the fluid system. Yet, the dynamics is inertial. To conserve the momentum and energy, the interface should slightly change its velocity. This causes the reactive force occurs and to stabilize the dynamics. The resulting flow in the conservative dynamics is a superposition of two motions: the background motion of the fluids following the interface, whose velocity slightly oscillates near the steady value, and the stable oscillations of the interface perturbations (Fig. 11).

#### 2. Interface destabilization for the classic Landau’s dynamics

For the classic Landau–Darrieus dynamics **r**_{LD}, the interface velocity is constant, $V\u0303=V\u03030$, as it is postulated by the boundary condition $u\u22c5n0=0$, leading to $u\u22c5n0+\theta \u0307\theta =0=0,\u2009v\u0303=0$. The constancy of the interface velocity prevents the reactive force to occur. The dynamics is destabilized. The LD unstable flow is a superposition of two motions—the background motion of the fluids following the interface with the constant velocity and the growth of the interface perturbations (Fig. 11).

Since, in the classic LD dynamics, mass and momentum are conserved, for capturing the destabilization mechanism, we consider how the condition $u\u22c5n0=0$ may influence the energy transport. Remarkably, for ideal incompressible fluids, the solution **r**_{LD} is incompatible with the condition for energy balance at the perturbed interface (Abarzhi *et al.,* 2019). Indeed, by substituting $u\u22c5n=0\u2009\u2009jn=0$ in the condition for the perturbed energy balance $Jnw+J\u22c5j/\rho 2=0$, one obtains $Jnw+J\u22c5j/\rho 2=Jnw=0$ and, with [*J*_{n}] = 0, reduces it further to [$w$] = 0 (Abarzhi *et al.,* 2015).

The enthalpy perturbations are $w$ = *p*/ρ + *δe*, where *δe* are the fluctuations of internal energy (in physics sense). In ideal incompressible fluids, free from energy sources, these fluctuations are zero, *δe* = 0, because $\u0117=0,\u2009\u2207e=0$. Thus, with $w$ = *p*/ρ and with ρ_{h} ≠ ρ_{l}, the condition for energy balance, [$w$] = 0, contradicts the condition for momentum balance, [*p*] = 0. We see that for ideal incompressible fluids, the classic Landau’s dynamics, while conserving the energy (and the physics enthalpy) to the zero order, [*W*_{0}] = 0, requires energy imbalance at the interface to the first order. This imbalance is induced by the first order enthalpy perturbations, [$w$] = [*p*/ρ]. This is the work done by the fluid when a parcel of fluid of a unit mass expands its volume from $\rho h\u22121$ to $\rho l\u22121$ under pressure *p*.

In realistic fluids, this imbalance of the perturbed energy can be induced by energy fluctuations. The effect can be self-consistently derived from entropy conditions accounting for chemical reactions. In ideal fluids, to evaluate the effect of the imbalance of perturbed energy on the interface stability, we may introduce an artificial energy flux and study a transition from stable to unstable dynamics with an increase in its strength (Abarzhi *et al.,* 2015).

### C. Conservative system with effect of energy fluctuations

The theoretical framework of Landau (1944) describes the evolution of a phase boundary: On each side of the interface, the fluids are ideal and incompressible. They are transformed into one another at the interface. The process is accompanied by a mass flow across the interface. At the planar interface, the fluxes of mass, momentum, and energy are fully balanced. For the perturbed interface, in some circumstances, the interface becomes unstable. The instability is accompanied by the appearance of a vortical field and requires the constancy of the interface velocity, the complete balance of the perturbed momentum and, for ideal incompressible fluids, an imbalance of perturbed energy, or energy fluctuations by the perturbed interface. Here, we model the effect of energy fluctuations on the interface stability and the flow fields’ structure (Abarzhi *et al.,* 2015).

Mathematically, it is straightforward to verify that adding a stationary energy source at the moving interface in the governing equations, such as $j\u0303\u22c5nW+j\u03032/2\rho 2\u2212Q\u0303=0$, keeps the boundary conditions and the interfacial dynamics the same, except for the modification of the zero-order enthalpy as $W0\u2192W\u03030$, assuming that the fluids are ideal and incompressible. To influence the interface stability, the energy imbalance, or energy fluctuations, should be produced by the perturbed interface.

For realistic fluids, the flux of energy fluctuations occurs due to the difference in the thermodynamic properties of the fluids (Kadanoff, 2000; Landau and Lifshitz, 1987). It should be derived self-consistently from the condition for entropy conservation (Landau and Lifshitz, 1987). For ideal incompressible fluids, we model the effect of energy fluctuation by introducing an artificial energy flux produced by the perturbed interface. We modify the boundary conditions as (Abarzhi *et al.,* 2015)

This model [Eqs. (17)] can be used only for estimates of the effect of the perturbed energy imbalance, or the energy fluctuations, on the interface stability. A self-consistent study is addressed to the future.

#### 1. Fundamental solutions

For the conservative system with energy fluctuations, with $Qhl=qhlVh2$ in Eqs. (17), the matrix Μ is $M=M\u0303$,

Its rank is 4. It has four eigenvalues ω_{i} and four corresponding eigenvectors **e**_{i}. With *q*_{h} − *q*_{l} = *q*, its determinant is $detM\u0303=iR\u221212\omega \u2212R\omega +R\omega 2+R+iq\omega \u2212RR+1\omega 2+2R\omega \u2212RR\u22121$. The parameter *q*, *q* > 0, describes the strength of energy fluctuations (the perturbed energy imbalance).

With $detM\u0303=RdetM\u2212qR/R\u22121detL$ and $detM\u0303=RdetM\u2212q/\omega R/R\u22121detL$, the condition $detM\u0303=0$ leads to (*R* − 1)*detM* = *qdetL* and $R\u22121detM=q/\omega detL\u0303$. Similarly to Abarzhi *et al.* (2015), we scale the fluctuations strength as $q=fR\u221212R/R+1$, where *f* ≥ 0 is a constant, and reduce equation $detM\u0303=0$ to

The dependence of fundamental solutions **r**_{i} = **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 4 on the density ratio *R* and the fluctuation strength *f* is cumbersome and not presented here. Solution **r**_{3} = **r**_{3}(ω_{3}, **e**_{3}) with ω_{3} = *R* and **e**_{3} = (0, *i*, 0, 1) is the same as for the conservative system, the classic Landau’s system, and the dynamic Landau’s system. Figure 12 illustrates the dependence of the eigenvalues on the fluctuation strength and the dependence of the eigenvalues on the density ratio for weak and strong fluctuations.

#### 2. Stability of fundamental solutions

In the limiting cases of weak (*f* → 0) and strong (*f* → ∞) fluctuations for fluids with very different densities (*R* → ∞) and with very similar densities (*R* → 1^{+}) in Eqs. (18), the properties of the dynamics are the following (Table VI and Fig. 12).

. | f ≡ 0
. | f > 0
. | 0 < f < f_{cr}
. | f_{cr} < f
. |
---|---|---|---|---|

r_{1} | ||||

Re[ω] | 0 | <0 | <0 | >0 |

C | * | 0 | 0 | * |

r_{2} | ||||

Re[ω] | 0 | <0 | <0 | <0 |

C | * | 0 | 0 | 0 |

r_{3} | ||||

Re[ω] | >0 | >0 | >0 | >0 |

C | * | * | * | * |

r_{4} | ||||

Re[ω] | <0 | <0 | <0 | <0 |

C | 0 | 0 | 0 | 0 |

. | f ≡ 0
. | f > 0
. | 0 < f < f_{cr}
. | f_{cr} < f
. |
---|---|---|---|---|

r_{1} | ||||

Re[ω] | 0 | <0 | <0 | >0 |

C | * | 0 | 0 | * |

r_{2} | ||||

Re[ω] | 0 | <0 | <0 | <0 |

C | * | 0 | 0 | 0 |

r_{3} | ||||

Re[ω] | >0 | >0 | >0 | >0 |

C | * | * | * | * |

r_{4} | ||||

Re[ω] | <0 | <0 | <0 | <0 |

C | 0 | 0 | 0 | 0 |

##### a. Weak fluctuations.

For weak fluctuations and very different densities,

In Eqs. (18c), solutions **r**_{i} with *i* = 1, 2 are stable and describe the decaying oscillations of the flow fields. Solution **r**_{3} is the same as in the foregoing cases. Solution **r**_{4} is stable (Fig. 12).

For weak fluctuations and similar densities,

##### b. Strong fluctuations.

For strong fluctuations and very different densities,

In Eqs. (18e), the fundamental solutions are similar to the corresponding solutions in the classic Landau’s system. Solution **r**_{1} is unstable and corresponds to the Landau–Darrieus instability. Solution **r**_{2} is formally stable. Solution **r**_{3} is the same as before. Solution **r**_{4} is stable (Fig. 12).

For strong fluctuations and similar densities,

#### 3. Flow fields of fundamental solutions

##### a. Weak fluctuations.

In Eqs. (18), in the case of weak fluctuations and very different densities, solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2 describe slightly decaying oscillations of the flow fields (Fig. 12). At the same time, the condition Re [ω_{1(2)}] < 0 leads to an increase in the vortical velocity component far from the interface; thus, the values of integration constants $C12,C1=C2*$ should be zero *C*_{1(2)} = 0 in order to satisfy the boundary conditions at the outside boundaries at all time. Solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed flow fields in the entire domain at any time for any *C*_{3} as before. For solution **r**_{4}(ω_{4}, **e**_{4}), due to Re[ω_{4}] < 0, the integration constant must be zero, *C*_{4} = 0, in order to satisfy the boundary conditions at the outside boundaries. In Eqs. (18), in the case of weak fluctuations and very similar densities, for solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 4, the corresponding integration constants *C*_{i} should be zero, *C*_{i} = 0, in order to satisfy all the boundary conditions. As before, solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed flow fields at any time for any *C*_{3}.

##### b. Strong fluctuations.

In Eqs. (18), in the case of strong fluctuations and very different densities, fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 3 are similar to the corresponding solutions in the classic Landau’s system. Specifically, (1) for solution **r**_{1}(ω_{1}, **e**_{1}), the interface perturbations are strongly coupled with the vortical and potential components of the velocity fields. This solution is unstable. (2) For solution **r**_{2}(ω_{2}, **e**_{2}), the interface perturbation and vortical and potential velocity components are also strongly coupled. While this solution is formally stable, its velocity field has a vortical structure that decays exponentially in time yet grows away from the interface because of Re[ω_{2}] < 0. For solution **r**_{2}(ω_{2}, **e**_{2}) to satisfy the boundary condition **u**_{l}|_{z→+∞} = 0 at all times, the integration constant must be zero *C*_{2} = 0. (3) Solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed flow fields in the entire domain at any time and for any *C*_{3}, as before. (4) For the formally stable solution **r**_{4}(ω_{4}, **e**_{4}) due to Re[ω_{4}] < 0, the integration constant must be *C*_{4} = 0 in order to satisfy the boundary condition at the outside boundaries of the domain at any time. In Eqs. (18), in the case of strong fluctuations and similar densities, the solutions are similar to those with very different densities. These similar features include the development of LD-type of instability for solution **r**_{1}(ω_{1}, **e**_{1}), the choice of *C*_{2} = 0 for solution **r**_{2}(ω_{2}, **e**_{2}) with Re[ω_{2}] < 0, the unperturbed flow fields for solution **r**_{3}(ω_{3}, **e**_{3}), and the choice of *C*_{4} = 0 for solution **r**_{4}(ω_{4}, **e**_{4}) with Re[ω_{4}] < 0.

#### 4. Interface velocity for fundamental solutions

The properties of fundamental solutions for the conservative dynamics with the effect of energy fluctuations indicate that for *f* > 0, the interface velocity in the laboratory frame of reference is constant, $V\u0303=V\u03030$, for these solutions in Eqs. (18) in the case of either strong *f* → ∞ or weak *f* → 0^{+} fluctuations. Indeed, solution **r**_{3} has the null perturbed fields of the velocity and the interface at any integration constant *C*_{3}. For stable solutions **r**_{2(4)}, the integration constants should be zero to satisfy the boundary conditions at the outside boundaries of the domain, *C*_{2(4)} = 0. For solution **r**_{1} the integration constant should be zero in the case of weak fluctuations, *f* → 0^{+}, *f* ≠ 0, in order to satisfy the boundary conditions at the outside boundaries of the domain; in the case of strong fluctuations *f* → ∞, this solution approaches the classic LD solution **r**_{1} → **r**_{LD} and has a constant interface velocity $V\u0303=V\u03030$, as discussed before. This result holds even when the boundary conditions away from the interface are slightly perturbed, since, for these solutions, the vortical component of the velocity grows exponentially far from the interface.

In the limit of zero fluctuations, *f* ≡ 0, the interface velocity experiences stable oscillations near the steady value $V\u03030$. These oscillations enable the inertial stabilization mechanism of the conservative dynamics to occur, as found in the foregoing.

#### 5. Formal properties of the conservative system with effect of energy fluctuations

To study the formal properties of the dynamic Landau’s system, we find for matrix $M=M\u0303$ the associated matrices $S=SM\u0303$ and $P=PM\u0303$,

Equations $detPM\u0303\u22121SM\u0303\u2212\omega I=0$ and $detM\u0303=0$ yield the same eigenvalues $\omega =\omega i,\u2009i=1,...,4$. The conservative system with the effect of energy fluctuations is non-degenerate and has four fundamental solutions for four degrees of freedom.

#### 6. Summary of properties of fundamental solutions for the conservative system with effect of energy fluctuations

##### a. Outline of properties of fundamental solutions.

For strong and weak fluctuations and for fluids with similar and contrasting densities, there are four fundamental solutions **r**_{i}(ω_{i}, **e**_{i}) with *i* = 1, 2, 3, 4 [Eqs. (18), Table VI, and Fig. 12]. Solution **r**_{3}(ω_{3}, **e**_{3}) corresponds to the unperturbed flow fields for any *C*_{3} and any *f*. In the case of strong fluctuations, *f* → ∞, solution **r**_{1}(ω_{1}, **e**_{1}) describes (for *C*_{1} ≠ 0) unstable dynamics of LD-type, whereas for solutions **r**_{2}(ω_{2}, **e**_{2}) and **r**_{4}(ω_{4}, **e**_{4}), the constants are *C*_{2(4)} = 0. In the case of weak fluctuations, *f* → 0, for solutions **r**_{1}(ω_{1}, **e**_{1}), **r**_{2}(ω_{2}, **e**_{2}), and **r**_{4}(ω_{4}, **e**_{4}), the constants are *C*_{1} = *C*_{2} = *C*_{4} = 0. For *f* > 0, the interface velocity in the laboratory frame of reference is constant, $V\u0303=V\u03030$, for these solutions.

At *f* ≡ 0, solutions **r**_{1}(ω_{1}, **e**_{1}) and **r**_{2}(ω_{2}, **e**_{2}) with non-zero integration constants *C*_{1}, *C*_{2} describe stable oscillations of the flow fields for the conservative dynamics, accompanied by stable oscillations of the interface velocity near the steady value $V\u03030$.

##### b. Transition from the stable to the unstable dynamics.

Transition from the stable to unstable dynamics is illustrated in Fig. 12 showing eigenvalues ω_{i}(*R*, *f*) with *i* = 1, 2, 3, 4 in Eqs. (17) and (18) for a fixed *R* as a function of *f* as well as in cases of small and large *f*. For very small *f*, the values Im[ω_{1(2)}] and ω_{3(4)} are indistinguishable from the corresponding values in conservative dynamics. For very large *f*, eigenvalues ω_{i} with *i* = 1, 2, 3 are indistinguishable from the corresponding values in the classic Landau’s system.

The properties of the flow fields of fundamental solutions for the conservative system with energy fluctuations indicate that when the value *f* changes from 0 to ∞ in the system [Eqs. (17) and (18)], an abrupt transition can occur to the unstable dynamics. The critical value *f* = *f*_{cr}, at which the transition is expected, can be estimated from condition ω = 0 in Eq. (18), leading to *f*_{cr} = (*R* + 1)/(*R* − 1), *q* = *q*_{cr} = *R*(*R* − 1), and $Qh\u2212Qlcr=RR\u22121Vh2$.

The flow dynamics can behave as follows: At *f* ≡ 0, the flow fields experience stable oscillations near the uniform. For 0 < *f* < *f*_{cr}, the flow fields are uniform. For *f* > *f*_{cr}, the dynamics becomes unstable. For large fluctuations, *f* → ∞, the unstable dynamics fully recreates the properties of the Landau–Darrieus instability.

In realistic fluids, other mechanisms are possible that account for the influences of entropy flux and the thermodynamic properties of the fluids on energy fluctuations produced by the perturbed interface (Kadanoff, 2000; Landau and Lifshitz, 1987) as well as for nonlinear effects (Sivashinsky, 1983).

## V. COMPARISON WITH TRADITIONAL STUDIES OF COMBUSTIBLE SYSTEMS

In this section, we briefly discuss the connection of our work to the studies in the theory, simulations, and experiments of combustion. We propose new experiments and experimental diagnostics on the basis of our analysis.

### A. Landau theoretical framework and traditional theories, simulations, and experiments of combustion

While the theoretical framework (Landau, 1994) describes the dynamics of a phase boundary broadly defined, LDI is traditionally associated with the theory of combustion. Landau (1944) considered the evolution of a flame in premixed combustion by treating the unburned gas as a heavy fluid, the burned gas as a light fluid, both fluids as ideal and incompressible, and the interface as a discontinuity.

A premixed flame is a self-sustained wave of an exothermic reaction in a gaseous mixture (Sivashinsky, 1983; Williams, 1965; and Zeldovich, 1944). The wave is initiated by a spark, propagates through the mixture, and produces combustion products. A freely propagating flame is a multi-scale, multi-parameter, multi-physics, and multi-chemistry system. Its evolution is a long-standing problem (Peters, 2000). For planar and steady flames, the interface dynamics is well captured by rigorous theories (Mayo *et al.,* 1990; Frank-Kamenetsky, 1969; Zeldovich and Frank-Kamenetsky, 1938; and Mallard and Le Chatelier, 1883). For unsteady and curved flames, the flame interaction with the hydrodynamic flow that it produces is a challenging problem (Class *et al.,* 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Liberman *et al.*, 1994; Sivashinsky, 1983; Matalon and Matkowsky, 1982; and Williams, 1965).

Landau’s result on the unconditional instability of the flame front contradicts some experiments because stable flames exist and are observed in a laboratory (Clavin, 2000; Clanet and Searby, 1998; Buckmaster, 1979; and Markstein, 1951). Several approaches were developed to explain these observations. Theoretical research was focused on bridging the gap between the idealistic framework and realistic conditions and on augmenting the theory (Landau, 1944) with physical effects typical for flames in premixed combustion (Clavin and Searby, 2016; Class *et al.,* 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Liberman *et al.,* 1994; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Williams, 1971; Sivashinsky, 1983; Matkowsky and Sivashinsky, 1979; Istratov and Librovich, 1966; and Zeldovich, 1944). These effects included the finite thickness of the flame and the finite curvature of the interface (Markstein, 1951), the molecular diffusion and the thermal expansion (Williams, 1971), the dynamics of a reaction zone inside a thin flame (Matalon and Matkowsky, 1982), and the kinematic viscosity (Frankel and Sivashinsky, 1982) as well as the assumptions of weak compressibility and a gradual change in the flow quantities across the interface (Clavin, 2000; Williams, 1971). These fundamental works have found that various physical effects can stabilize the development of LDI at small scales (Sivashinsky, 1983). Peters (2000) and Clavin (2000) summarized these studies, reconciled the theory and the experiment, and found that in premixed combustion, the flame stability depends on the interplay of compressibility, mass diffusion, and heat diffusion. Class *et al.* (2003a; 2003b), Klimenko and Class (2000), and Sivashinsky (1983) analyzed the stability of a uniformly propagating planar flame as a solution of the unified model that accurately captured Landau–Darrieus, cellular, and pulsating instabilities including a high-wavenumber cutoff for each.

Our results agree with Landau’s results (Landau, 1944) and thus substantiate that the theories of premixed combustion are accurate and complete, within the range of their applicability (Class *et al.,* 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Sivashinsky, 1983; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Sivashinsky, 1980; Matkowsky and Sivashinsky, 1979; Williams, 1971; Markstein, 1951). Particularly, by applying proper transformations of variables, one can directly match our results and the classical results (Landau and Lifshitz, 1987). This includes the governing equations, the boundary conditions, the structure of the solution, the characteristic length scale of the vortical field, and the instability growth-rate. Furthermore, the flow fields in Figs. 8 and 9 can be obtained with some effort from classical works (Peters, 2000; Landau and Lifshitz, 1987; Sivashinsky, 1983; Landau, 1944; Zeldovich, 1944).

In experiments, the implementation and diagnostics of Landau–Darrieus instability is an extreme challenge. It is remarkable that the first experimental study of LDI (Clanet and Searby, 1998) was reported nearly 60 years after the theoretical work was published (Landau, 1944). These experiments have measured the LDI growth-rate in the compressible reactive fluids (gases) for a relative thick interface (Markstein parameter 4.0) in the presence of acceleration (gravity) under the influence of acoustic parametric forcing. More recently, the observations of Landau–Darrieus instability have also been reported for a Bunsen burner and for an inverted V-flame (Clavin, 2000). These experiments have observed that under the influence of Landau–Darrieus instability, the amplitude of the flame wrinkles increases; the instability growth-rate agrees with the theory accounting for the finite flame thickness and non-ideal effects (Markstein parameters 4.3, 4.6, 5.0) (Clavin, 2000). A direct quantitative comparison of these experimental results and our theory is a challenge since the experimental conditions strongly deviate from our theoretical approximation (Clavin, 2000; Clanet and Searby, 1998).

For numerical modeling of Landau–Darrieus instability of combustion fronts in premixed gases, the simulations apply the so-called low Mach number approximation and consider non-ideal fluids with simplified chemical kinetics and gradually changing flow fields (Schlimpert *et al.,* 2016; Miglani *et al.,* 2014; Clavin, 2000; and Denet and Haldenwang, 1995). Particularly, simulations (Denet and Haldenwang, 1995) adopt the isobaric approximation valid for low Mach number flows, apply momentum–pressure formulation for variable density flows, and solve numerically the equations for a premixed flame in a continuous non-ideal fluid. The simulations obtain curved flames with large amplitudes and provide measurements of the growth-rate of LDI (Denet and Haldenwang, 1995). This dispersion relation is consistent with theory (Clavin, 2000; Pelce and Clavin, 1982) and detects the stabilization of short wavelengths.

In addition to the quantitative results in a broad parameter regime, the simulations (Denet and Haldenwang, 1995) also present the flow fields, including temperature, velocity streamlines, vorticity lines, and the interfacial shear. The presence of the vortical field in the bulk of the light fluid is consistent with traditional theories of combustion and with our results for LDI (Fig. 8) (Clavin, 2000; Williams, 1971; Landau, 1944). The simulations also report the presence of shear at the interface with a finite thickness (Denet and Haldenwang, 1995), whereas for Landau’s solution, the velocity field is shear free at the discontinuous interface [Eqs. (2) and (14) and Fig. 8] (Landau and Lifshitz, 1987). The absence of interfacial shear at the discontinuous interface with interfacial mass flux is a rigorous theoretical result for ideal fluids [Eqs. (2)] (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018; and Landau and Lifshitz, 1987). Further investigations are required to better understand the processes of formation of shear at the interface and the vortical field in the bulk in real fluids.

Numerical simulations of premixed combustion, while employing realistic conditions of weakly compressible non-ideal fluids, usually presume the constancy of the interface velocity (Schlimpert *et al.,* 2016; Clavin, 2000; Denet and Haldenwang, 1995). The direct observation of the inertial stabilization mechanism, which leads to the stable oscillations of the interface velocity near the constant value, and which our analysis finds for ideal incompressible fluids, may be a challenge to detect in these simulations. Theories of premixed combustion also often presume the constancy of the interface velocity (Class *et al.,* 2003a; 2003b; Clavin, 2000; and Peters, 2000). For instance, this assumption is applied in the studies of the influence of temperature dependence of diffusivities on the dynamics of flame fronts. This analysis has derived the 4 × 4 matrix, has found the eigenvalues for the matrix, and has identified the instability growth-rate (Clavin, 2000). A derivation of eigenvectors for the 4 × 4 matrix may further enable the study of structure of flow fields, as in Eqs. (1)–(18).

### B. New experiments and experimental diagnostics of combustible systems

The existing experiments and simulations of the interfacial dynamics are focused on the evolution of the perturbation amplitude (Class *et al.,* 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Clanet and Searby, 1998; Sivashinsky, 1983; Matalon and Matkowsky, 1982; and Matkowsky and Sivashinsky, 1979). Our theoretical analysis, in addition to identifying the evolution of the perturbation amplitude, finds that the flow is highly sensitive to the interfacial boundary conditions and that the properties of interfacial transports at microscopic scales can be linked to the flow fields at macroscopic scales. Note that Figs. 3–6 and 8–10 present the perturbations of the flow fields. In order to obtain the flow fields, the corresponding zero-order component should be added to the perturbations.

Figure 13 illustrates the velocity streamlines **S**_{h(l)} of the heavy and light fluids for the conservative dynamics and the classic Landau’s dynamics at some density ratio and some instance of time. Each plot has its own range of values. The streamlines are defined as $dShl/dt\xd7vhl$, where the heavy (light) fluid velocity is $vhl=Vhl+uhl$. The integration constants for solutions **r**_{CD} and **r**_{LD} are chosen such that the amplitude of the interface perturbation is the same for the conservative dynamics and the classic Landau’s dynamics and equals to $z\xaf=1/3$. Figure 13 illustrates slight deviations of the velocity streamlines from the straight lines for both conservative dynamics and classic Landau’s dynamics. For the conservative dynamics, the deviations of the velocity streamlines are clearly seen in the heavy fluid; they are also present in the light fluid though they are somewhat less pronounced due to the large value of the uniform velocity of the light fluid. For the classic Landau’s dynamics, the deviations of the velocity streamlines are clearly seen in both fluids. For a given value of the perturbation amplitude, the fluid motion is more intense in the classic Landau’s dynamics (due to the presence of the vortical components of the velocity of the light fluid) when compared to that in the conservative dynamics. The results in Figs. 4, 8, and 13 are in consistency with one another.

Note that due to the small values of the perturbed velocities, $uhl\u226aVhl$, the properties of the flow fields in the fluid bulk may be a challenge to precisely diagnose. For addressing these challenges, advanced experimental technologies may be applied (Orlov *et al.,* 2010). New experiments with new diagnostics can be conducted on the basis of our results to better understand the interfacial dynamics. Here, we outline some parameters of a new experiment.

According to our results, in the experiment, the interface velocity is sub-sonic, the dynamics is observed from a far field, and its characteristic length scale 1/*k* = *λ*/2π and time scale *τ* = 1/*kV*_{h} = *λ*/2*πV*_{h} are greater than those that are induced by the processes of dissipation and diffusion and by a finite thickness of the interface. These scales are set by the initial conditions. To select the density ratio for the experimental fluid system, we note that, in the classical Landau’s solution for a given *λ*, the length-scale $\lambda \u0303$ of the vortical field is a non-monotone function of the density ratio *R* with $\lambda \u0303/\lambda =k/k\u0303=R1+R/\u2212R+\u2212R+R2+R3$. The minimum value of $\lambda \u0303/\lambda $ is achieved at $R=2+5\u22484.24$. Thus, the vortical field can be easier to diagnose when the fluids density ratio is $R=2+5\u22484.24$ and the ratio of uniform velocities is $Vl/Vh=2+5\u22484.24$.

In this case, the Landau solution has the following characteristics:

or, numerically, *R* = 4.24, $\lambda \u0303/\lambda =4.24$, ω = 1.00, and **e**_{1} = (−0.62 *i*, 1.6 *i*, 0.62 *i*, 1). The choice of the integration constant for the solution is suggested by the requirement that the perturbations of the flow fields are small. By evaluating the amplitude of the interface perturbation for which linear theory is valid as $z\xaf=1/3$, we estimate the integration constant as $C=5.4\xd710\u22121$. By implementing in experiments the fluid density $R=2+5\u22484.24$ and the corresponding initial conditions, one can further diagnose the flow fields in the vicinity of the interface and in the bulk and compare with the theoretical solutions.

The experimental outcome may include the following scenarios: (1) The interface is stable; the perturbed flow fields are the same as those of the conservative solution; they slightly oscillate near the equilibrium. This can be interpreted as the absence of energy fluctuations and/or as their weak effect. (2) The interface is unstable; the flow fields are the same as in the classical Landau’s solution. This can be interpreted as the existence of LDI and the applicability of Landau’s framework for the experimental system. (3) The interface is unstable; the flow fields deviate qualitatively from those of the classical LDI and from those of the conservative dynamics with strong energy fluctuations. This can be interpreted as that the approximation of ideal incompressible fluids is not a valid approximation for the fluid system. It may also suggest that the fluid system experiences a hydrodynamic instability that is different from LDI. If the growth-rate of this instability agrees with the results of the seminal theories of combustion, it may give grounds to name those instabilities in honor of those who studied its properties (Class *et al.,* 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Clanet and Searby, 1998; Liberman *et al.,* 1994, Sivashinsky, 1983; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Matkowsky and Sivashinsky, 1979; Williams, 1971; and Markstein, 1951).

To even more accurately outline the regime of experimental parameters, our analysis can be expanded to include compressibility and non-ideal effects. We address these studies to the future.

## VI. CHEMICALLY REACTIVE SYSTEMS

In reactive fluids, it is generally well understood that chemically reactive systems can be hydro-dynamically unstable (Ilyin *et al.*, 2019; Zhakhovsky *et al.,* 2014; Mayo *et al.,* 1990). Yet, it is a challenge to construct a model experiment or a simulation that clearly displays a significant chemical reaction instability at a simple interface. We can use atomistic simulations to study the energy transport in a reactive system and the effect of chemical reaction on the interface stability. Here, we briefly outline the properties of a simple chemically reactive system. By applying reactive molecular dynamics simulations, we find the following: (i) Simple chemical reactions can be accompanied by complex processes with transient states (Ilyin *et al.*, 2019). (ii) Energy fluctuations due to chemical reactions can destabilize the interface, leading to chemistry-induced instability. (iii) Energy scatter (via rotational and vibrational degrees of freedom rather than translational degrees of freedom) may hinder the instability development. Our results are illustrated in Figs. 14–17.

### A. Reactive molecular dynamics

In our reactive molecular dynamics simulations, we find that reactions that are simple at first glance may exhibit complex transient behavior (Ilyin *et al.*, 2019). As a particular example, we consider radical initiated hydrogen peroxide decomposition to water and oxygen. We start with the reactants hydrogen peroxide (HOOH) and hydroxyl radical (OH) and allow them to react under various temperatures and pressures to form a series of intermediate species (OH and HOO in steady-state as well as occasional OH–OH and HOO–HOOH complexes) to products (HOH and O_{2}). Here, we use standard notations O for oxygen and H for hydrogen.

### B. RMD2Kin method and computational setup

We analyze this system by applying the RMD2Kin method (Ilyin *et al.*, 2019). The molecular dynamics is run with the reactive force field ReaxFF, which allows electric charges to flow freely and permits continuous bond breaking and formation (van Duin *et al.,* 2001). Molecules are detected as they form and break bonds, and the chemical reactions are analyzed directly with no preconceived ideas about the chemistry. We follow the products back in time to find the specific reaction that formed them. We then statistically analyze all reactions at all time steps and at various temperatures to extract reaction rates and barriers using the Arrhenius relation (Ilyin *et al.*, 2019).

The computational setup is as follows: The system is enclosed in a flat periodic box 20 × 2 × 100 nm^{3} in size. At the start of the simulation, the box is filled with 100 000 peroxide molecules. After equilibration and heating, reactions are initiated by introducing radicals in a specific region [Fig. 15(a)] or by thermal shock in a specific region (Fig. 17). The system is then allowed to react for tens to hundreds of picoseconds.

From the reaction trajectories, the observed dynamics are quite complex. Many reactions can be described as two bodies colliding to form two products (OH + HOOH → H_{2}O + HOO) or three products (HOO + HOOH → O_{2} + H_{2}O + OH). For uni-molecular reactions (HOOH → OH + OH), third bodies play the role in energizing the reacting molecules without being modified. However, we also find additional interactions with no net reaction that are essential to accurately describe the reactive process. For instance, two OH radicals may end up in a rotationally bound state that temporarily reduces their concentration and affects the overall reaction rate. As another example, HOO and HOOH can transfer the H back and forth, which effectively removes some HOO radicals from interactions and also has to be accounted for.

We find five dominant species and seven dominant reactions. The five dominant species are HOOH, H_{2}O, O_{2}, HOO, and OH. The seven dominant reactions are (1) HOOH + OH → HOO + H_{2}O, (2) HOO + HOOH → O_{2} + H_{2}O + OH, (3) HOO + OH → O_{2} + H_{2}O, (4) HOO + HOO → O_{2} + HOOH, (5) HOO + H_{2}O → OH + HOOH, (6) HOOH + M (a third body) → OH + OH, and (7) OH + OH <=> HO–OH (see Fig. 14). We find that each reaction transfers the translational motion of the reactants into rotational and vibrational energy of the products. Thus, the kinetic energy of the reactants is effectively absorbed into molecular degrees of freedom.

### C. Chemistry-induced instability

To study the effect of energy fluctuations on the interface dynamics, we carry out atomistic simulations aimed specifically at searching for chemically induced instabilities such as LDI. We use two ways of initiating the reactions: introduction of radicals (Figs. 15 and 16) and thermal shock (Fig. 17).

When the radicals are introduced, as in Figs. 15 and 16, then, at the interface between a region with HOOH seeded with OH and region of pure HOOH, we observe the initiation of a local plume developing into HOOH. While the plume starts to grow, it, however, is hindered before developing a clear instability (Figs. 15 and 16). When the thermal shock is introduced, as in Fig. 17, then, after some growth, small perturbations along the thermal shock interface are hindered as well (Fig. 17). We observe that at late times, both methods produce a stable front propagation.

In Fig. 15, energetic radicals are propelled forward and initiate the formation of a plume at 10 ps; however, the plume fails to propagate at 12 ps due to energy scatter into the bulk, and by 24 ps, the plume has effectively disappeared into the rest of the flame front. Similarly, in Fig. 16, a perturbation is formed around 1 ps and grows over the next 16 ps. However, it does not propagate forward, and by 24 ps, the rest of the front “catches up” and the perturbation disappears (Fig. 16). A thermal shock in the same region also does not lead to instability. The temperature map (Fig. 17) shows the location of high-energy molecules, and the temperature front propagates in a stable manner.

In all systems studied, small local perturbations appear but then quickly disappear into the rest of the reaction front. As a result, the front is not perfectly flat but propagates in a stable way. Any initial forward momentum of the front is rapidly scattered, and any initial perturbations that may have led to an LDI-like instability completely disappear.

We hypothesize two main reasons for the hindrance of perturbations. In order to develop into an instability, a perturbation needs (1) forward momentum and (2) exponential growth. First, the momentum that could potentially cause formation and propagation of a plume comes from the individual radicals. These radicals are smaller than the bulk HOOH molecules, and thus, any initial forward momentum is rapidly scattered after a few collisions. Second, the reactions in this particular chemistry are such that each reaction consumes and produces one radical. There is no multiplication of the number of radicals or chain reaction. In fact, studies of similar chemistry have shown that the number of radicals remains constant throughout (Ilyin *et al.*, 2019). Therefore, the perturbation does not experience exponential growth, and the evolution of the interface is stabilized.

These issues may be solved by (i) introducing outside forces that may provide the necessary forward momentum, (ii) selecting chemistry that yields more easily a radical chain reaction, (iii) including fast propagation of high energy atoms and radicals, and (iv) selecting the proper initial conditions for the flow fields. For instance, we might see more radical multiplication that could lead to the exponential growth of the plume by changing the reactants or parameters of the system to an extreme, e.g., by replacing HOOH with analogs that produce two or more radicals per each consumed radical, such as HOOOH or HOOOOH, or by using a different chemistry altogether. To further disturb the flat interface, the reactive system may be augmented with the fast propagation of high energy atoms and radicals, which, in turn, may initiate the reaction and the formation of hotspot ahead of the interface.

Emphasize that in our reactive molecular dynamics simulations, the energy is effectively scattered via rotational or vibrational degrees of freedom (rather than diffused). This process is substantially quicker than the energy dissipation (diffusion) via translational degrees of freedom, which require multiple collisions. Our theoretical analysis at continuous scales suggests the importance of the large-scale vortical field in the bulk of the reacted material (in order for the instability to occur). Its implementation in reactive molecular dynamics simulations is challenging, and we address the resolving of this challenge to the future research.

Hence, atomistic simulations of combustible systems find that even a simple chemical reaction can be accompanied by complex processes with transient states strongly influencing energy transport. Energy fluctuations by chemical reactions can potentially lead to chemistry-induced instabilities. The hindering of the instability can occur via energy scatter locally into molecular degrees of freedom (rather than the longer-range energy diffusion due to multiple collisions). Further investigations are required to better understand the properties of chemistry-induced instabilities at atomic and continuous scales. Our analytical and numerical methods can serve for accurate modeling of reactive systems and their hydrodynamic stability and for expanding the theory of combustion beyond the diffusive approximation.

## VII. INTERFACE DYNAMICS INFLUENCED BY ACCELERATION

The inertial stabilization mechanism of conservative dynamics suggests that for ideal fluids, the interface velocity experiences small stable oscillations near the constant value. In realistic fluids, these free oscillations can be damped by dissipation, diffusion, and other effects and can lead to the constancy of the interface velocity. A question thus appears whether the inertial stabilization mechanism can be directly observed in experiments and simulations. To answer this question, we analyze the interface dynamics influenced by acceleration. Here, we briefly outline the properties of accelerated conservative dynamics and compare them with those of inertial conservative dynamics as well as with accelerated Landau–Darrieus and Rayleigh–Taylor dynamics. Details can be found in Ilyin *et al.* (2018).

Under the influence of acceleration **g** directed from the heavy to the light fluid, **g** = (0, 0, 1), pressure is modified as *P* → *P* + ρ *g z*. Accounting for the structure of the solution [Eqs. (10) and (11)], the pressure perturbations are transformed to $phl=\u2212\rho hl\Phi \u0307hl+Vhl\u2202\Phi hl/\u2202z\u2212gz$. By introducing the dimensionless acceleration value $G=g/kVh2$ with *G* > 0, one can further find, for accelerated conservative dynamics, the matrix *M*_{G} and the corresponding fundamental solutions. While matrices *M*_{G} and *M* appear similar [except for element 33, which is *G*(*R* − 1) in matrix *M*_{G} and is 0 in matrix M], a number of properties of accelerated conservative dynamics deviate substantially from those of inertial conservative dynamics. Particularly, accelerated conservative dynamics exhibits a new fluid instability (Ilyin *et al.,* 2018).

Per condition *detM*_{G} = 0, the fundamental solutions for accelerated conservative dynamics are

where *G*_{cr} = *R*(*R* − 1)/(*R* + 1) and asterisks mark functions of *R*, *G*. Fundamental solutions **r**_{1(2)} have potential velocity fields in the fluid bulk for any *G* ≥ 0 and are transformed to the corresponding solutions in Eqs. (12) at *G* = 0. Fundamental solutions **r**_{3(4)} are the same as the corresponding solutions in Eqs. (12) for any acceleration value *G* (Figs. 5 and 6) (Ilyin *et al.,* 2018).

For small acceleration values *G* < *G*_{cr}, solutions **r**_{1}, **r**_{2} describe stable traveling waves whose interference results in stably oscillating standing waves, similarly to Figs. 3 and 4. For large acceleration values *G* > *G*_{cr}, solutions **r**_{1(2)} describe a standing wave with an exponentially growing (decaying) amplitude. The stable and unstable regions in the (*G*, *R*) plane are separated by the curve *G* − *R*(*R* − 1)/(*R* + 1) = 0 and the curve $g\u2212kVh2\rho h/\rho l\rho h\u2212\rho l/\rho h+\rho l=0$ in dimensional form. This suggests that the stability of accelerated conservative dynamics is defined by a balance of inertial and acceleration effects, that is, by reactive force and gravity. For *G* < *G*_{cr}, the inertial effect ∼$kVh2$ dominates and solution **r**_{CDG} = (**r**_{1} + **r**_{2})/2 is stable. In this regime, the interface is stable and the interface velocity experiences stable oscillations near the constant value, similarly to the case of inertial conservative dynamics [Eqs. (12)]. For *G* > *G*_{cr}, acceleration effect ∼*g* dominates and solution **r**_{CDG} = **r**_{1} is unstable. In this regime, the interface is unstable and the interface velocity increases with time (Ilyin *et al.,* 2018).

The properties of this new fluid instability deviate substantially from those of the accelerated Landau–Darrieus instability and Rayleigh–Taylor instability (Ilyin *et al.,* 2018). Particularly, the new fluid instability develops for acceleration values greater than a threshold, *G* > *G*_{cr}, whereas the accelerated Landau–Darrieus and Rayleigh–Taylor instabilities develop for any acceleration value *G* > 0. For accelerated conservative dynamics, the velocity field is potential in the bulk, similarly to Rayleigh–Taylor unstable dynamics and unlike accelerated Landau–Darrieus dynamics having the coupled potential and vortical components of the velocity field. The accelerated conservative dynamics is shear-free at the interface, similarly accelerated Landau–Darrieus dynamics and unlike Rayleigh–Taylor unstable dynamics exhibiting the interfacial shear. The interface velocity is time-dependent for accelerated conservative dynamics, is constant for accelerated Landau–Darrieus dynamics, and is zero for Rayleigh–Taylor unstable dynamics. At a given density ratio *R* > 1 and at large enough acceleration values *G* > *G*^{*}, where *G*^{*} = (*R*^{2} − 1)/4 with *G*^{*} > *G*_{cr}, the new fluid instability has the largest growth-rate when compared to the accelerated Landau–Darrieus instability and Rayleigh–Taylor instability (Ilyin *et al.,* 2018).

We emphasize that the velocity fields have the same qualitative features for the accelerated and inertial fluid systems. For conservative dynamics, the velocity field is potential in the bulk and is shear-free at the interface. For Landau’s dynamics, the velocity field has the coupled potential and vortical components in the bulk and is shear-free at the interface. Pressure fields are, however, distinct in the inertial and accelerated cases (Figs. 4 and 8) (Ilyin *et al.,* 2018). In the presence of acceleration, the perturbed pressure fields become asymmetric in the light and heavy fluids. This anisotropy leads to the formation of the structure of bubbles and spikes at the interface and is present in the accelerated conservative dynamics and accelerated Landau–Darrieus and Rayleigh–Taylor dynamics (Ilyin *et al.,* 2018).

## VIII. MULTIPHASE FLOWS

Our theoretical framework is developed to study interface dynamics in a multiphase flow. The interface is broadly defined as a phase boundary. It can be an interface between two distinct fluids, an interface between distinct phases of the same fluid, or an interface at which fluids experience a phase transition and change their thermodynamic properties (Abarzhi *et al.,* 2019). In addition to combustible systems, some examples of multiphase flows with interfacial mass flux may include, e.g., explosion of type-Ia supernova, convection in stellar interiors, indirect-drive inertial confinement fusion, and transportation of liquefied natural gas (Abarzhi *et al.,* 2019; Karthik *et al.,* 2018; Arjomandnia *et al.,* 2014; Hurricane *et al.,* 2014; Nepomnyashchy *et al.,* 2012; Bo *et al.,* 2011; Gorokhovski and Herrmann, 2008; Prosperetti, 2017; Bell *et al.,* 2004a; 2003b; Stein and Norlund, 2000; Clavin, 2000; Peters, 2000; and Sivashinsky, 1977).

Our theoretical framework can be employed to study multiphase flows with sharply changing flow fields in nature and technology (Abarzhi *et al.,* 2019; Karthik *et al.,* 2018; Arjomandnia *et al.,* 2014; Hurricane *et al.,* 2014; Nepomnyashchy *et al.,* 2012; Bo *et al.,* 2011; Prosperetti, 2017, Bell *et al.,* 2004a; 2003b; and Peters, 2000). Our analysis finds strong sensitivity of the interface dynamics to the interfacial boundary conditions, directly links microscopic processes at the interface to macroscopic properties of the flow fields in the bulk, and elaborates extensive theoretical benchmarks for experiments and simulations Eqs. (1)–(18) and Figs. 1–12. To the best of the authors’ knowledge, these benchmarks were not discussed before. They can serve for advancement of methods of experimental diagnostics and numerical modeling of multi-phase flows in a broad parameter regime and can expand the capabilities of accurate methods of numerical modeling of the interface dynamics, including the level set and the in-phase fluid method (Karthik *et al.,* 2018; Bo *et al.,* 2011; Bell *et al.,* 2004a; 2003b; Badalassi *et al.,* 2003; Osher and Fedkiw, 2002; Peters, 2000; and Anderson *et al.,* 1998).

Our theoretical approach accurately describes the formal properties of the interface dynamics, finds the degeneracy of the classic Landau’s solution, and connects the lifting of this degeneracy to a singular, possibly—self-similar, character of the dynamics. Our theory finds that for the conservative dynamics, the interface between two fluids, which may mix at small scales due to the interfacial mass flux (i.e., are miscible at molecular scales), can be stable at global scale. One possible application of our results in multiphase geophysical flows can be a qualitative explanation of confluence of rivers—when waters in two rivers meet and not mix at the global scale (Fig. 18).

In order to study realistic multiphase flows, including premixed combustion, accounting for compressibility and non-ideal effects, our analysis can be augmented with the conditions for entropy conservation, the corresponding boundary conditions, and the equation of state. The formalism developed in this work provides the principal opportunity to systematically study the influence of realistic properties of the fluids on the quantitative and qualitative properties of interface dynamics and the structure of flow fields. We address these studies to the future.

## IX. DISCUSSION

We have studied the dynamics of the interface separating ideal fluids with different densities and having interfacial mass flux [Eqs. (1)–(3)]. We have developed and applied the general matrix method to solve the boundary value problem for the system of governing equations linearized near the equilibrium state of the uniform flow fields [Eqs. (1)–(11)]. We have formally derived the fundamental solutions, have found their eigenvalues and eigenvectors, and have directly linked the interface stability to the structure of the flow fields (Figs. 1–13 and Tables I–VI). Our theoretical framework is consistent with and generalizes the classic approach (Landau and Lifshitz, 1987) and assumes sharp changes in the flow quantities at the interface and the negligible effects of dissipation, diffusion, compressibility, surface tension, and the interface thickness. We have identified the extreme sensitivity of the interfacial dynamics to the interfacial boundary conditions [Eqs. (12)–(18)]. We have found the new properties of the interface dynamics for the conservative and the Landau systems and the new mechanisms of the interface stabilization and destabilization [Eqs. (12)–(18) and Figs. 1–13].

Particularly, the conservative dynamics is stable, and it is stabilized by the reactive force and the inertial effects. The stable solutions are traveling waves; their interference results in the stably oscillating standing waves. The stable flow is a superposition of two motions—the background motion of the fluids following the interface, whose velocity slightly oscillates near the steady value, and the stable oscillations of the interface perturbations [Eqs. (12) and Figs. 1, 2, 4, and 11]. For the classic Landau’s dynamics, the interface velocity is constant, the reactive force is absent, and the dynamics is unstable. The unstable solution is the standing wave with the increasing amplitude. The LD unstable flow is a superposition of two motions—the background motion of the fluids following the interface with the constant velocity and the growth of the interface perturbations [Eqs. (14) and Figs. 1, 7, 8, and 11]. For the stable conservative dynamics, the velocity fields are potential in the fluids’ bulk (Fig. 4). For the classic unstable LD dynamics, the potential and vortical components of the velocity fields are strongly coupled with the interface perturbation (Fig. 8). In either case, there is the null velocity shear at the interface. In either case, the velocity streamlines slightly deviate from straight lines, though for the unstable Landau’s dynamics, the deviations are stronger when compared to those in the conservative dynamics (Fig. 13).

Our results clearly indicate that the problem of the Landau–Darrieus instability remains a source of inspiration for theory research. Landau (1944) analyzed the interface stability by applying compatibility conditions and identifying two eigenvalues of the dynamics (Landau, 1944). Over the past 70 years, the theoretical framework of Landau (1944) has been extended to incorporate the properties of realistic environments (Class *et al.,* 2003a; 2003b; Peters, 2000; Clavin, 2000; Sivashinsky, 1983; Matalon and Matkowsky, 1982; Williams, 1971; and Markstein, 1951). Our work focuses on general conditions of the stability of the interface with the interfacial mass flux, considers incompressible ideal fluids from a far field, and treats the interface as a discontinuity [Eqs. (1)–(16)]. We have developed and employed the general matrix method and have found the formal fundamental solutions for the linearized dynamics, which directly link the interface stability to the structure of the flow fields’ (Figs. 1–13). To our knowledge, this rigorous consideration has not been performed before.

By applying our general matrix method, we have found that the conservative system is non-degenerate and has four fundamental solutions associated with four independent degrees of freedom [Eqs. (12), Tables I and III, and Figs. 2–6]. These include the two stable traveling waves, whose interference results in two stably oscillating standing waves, and the two other solutions—the unstable solution with the null perturbation fields of the velocity and pressure and the stable solution with the null integration constant. We have further found that the classic Landau’s system is degenerate and has only three fundamental solutions [Eqs. (14), Tables II–IV, and Figs. 5 and 7–9]. These include the unstable standing wave corresponding to LDI, the stable solution with the null integration constant, and the unstable solution with the null perturbation fields of the velocity and pressure, which is the same as for the conservative dynamics. The degeneracy is associated with the static nature of the special postulated condition at the interface (Landau, 1944) and can be lifted with the use of a more general condition.

When the degeneracy is lifted by modifying the postulated condition of zero perturbed mass flux by a constant mass flux across the interface [Eqs. (15)], the fourth solution appears [Eqs. (16), Tables II and V, and Fig. 10]. This solution is neutrally stable. It has a null eigenvalue, and its vortical field can seed the LDI. The null eigenvalue implies that the perturbations of the flow fields can be scale-invariant power-law functions of time rather than scale-dependent exponential functions. Furthermore, by noting that the interface is a phase boundary at which one fluid is transformed into another, this result is in agreement with the linear response theory of phase transitions, where thermodynamic fluctuations are known to be power-laws rather than exponentials (Landau and Lifshitz, 1987; Kadanoff *et al.,* 1967; Kadanoff, 2000).

We have identified the new stabilization mechanism of the conservative dynamics. It is due to the inertial effects that cause the reactive force to occur and, by influencing the interface velocity, to balance the change of momentum produced by the interface perturbation. Note that this inertial stabilization mechanism leads to the new hydrodynamic instability of the accelerated interface that develops when the acceleration is directed from the heavy fluid to the light fluid and when the acceleration magnitude exceeds a threshold (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018). Both in the inertial and accelerated conservative dynamics, the flow has potential velocity fields in the fluid bulk and is shear free at the interface (Abarzhi *et al.,* 2019; Ilyin *et al.,* 2018).

In LDI, the inertial stabilization mechanism is absent due to the perfect constancy of the interface velocity. Furthermore, in ideal incompressible fluids, the development of the LDI is associated with an imbalance of perturbed energy due to the energy fluctuations by the perturbed front. To model the effect of energy of fluctuations (i.e., the imbalance of the perturbed energy), we have introduced an artificial energy flux and have found that the LDI may develop when the fluctuations are strong. The structure of the flow fields of the solutions for the conservative dynamics with the effect of energy fluctuations suggests a possibility of an abrupt transition to unstable dynamics with the increase in fluctuation strength over a critical value (Fig. 12). In realistic systems, the energy imbalance can be induced by chemical reactions and by difference in the thermodynamic properties of the fluids (Abarzhi *et al.,* 2019; Ilyin *et al.*, 2019).

To our knowledge, these properties of the conservative dynamics and the Landau dynamics have not been discussed in detail in other studies (Abarzhi *et al.,* 2019). Traditionally, it is believed that an interface separating nearly ideal incompressible fluids and having an interfacial mass flux is a subject to the LDI at large scales and that in a realistic environment, the LDI is a challenge to implement because the effects of dissipation, diffusion, and finite interface thickness stabilize the small scales (Peters, 2000). Our analysis is in agreement with these results: For fluids with similar densities, the contribution of vortical fields to the dynamics is small, and the fluxes induced in realistic fluids by the stabilizing effects may dominate the dynamics and define the interface stability (Abarzhi *et al.,* 2015). For fluids with very different densities, a more careful consideration is required. Our analysis yields qualitative and quantitative characteristics of the dynamics that have not been measured before and that require new diagnostics (Figs. 1–13) (Abarzhi *et al.,* 2019; 2015).

One such experiment can be a study of the dynamics of fluids with very different densities, with diagnostics of the flow fields near the interface and in the fluids’ bulk, and with the measurements of the interface evolution, including the interface velocity as a whole and the interface perturbation growth-rate. By comparing the observations with our benchmarks, one can further identify the fundamentals of the interfacial dynamics in realistic environments and elaborate new approaches for the flow control (Figs. 1–13). Several questions may frame these studies: Can the LDI unconditionally develop at the large scales? How can the dynamics be stabilized–by inertial effect, by dissipation and diffusion, or by their combination? How strong should energy fluctuations be to destabilize the flow? Can these fluctuations be induced by chemical reactions? How are the properties of chemistry-induced instabilities compared to those of the LDI?

Our studies of chemically reactive flows by means of reactive molecular dynamics find the following: Even a simple chemical reaction can be accompanied by complex transient processes with strong energy fluctuations. Energy fluctuations by chemical reactions can lead to chemistry-induced instabilities. The hindering of the instability can occur via the energy scatter. Further investigations are required to better understand the properties of chemistry-induced instabilities at atomic and continuous scales.

The existing experimental and numerical studies of the interface stability are focused on the growth of the perturbation amplitude (Abarzhi *et al.,* 2019). Our analysis derives the amplitude growth-rate and finds that the dynamics is highly sensitive to interfacial boundary conditions. According to our theory, by measuring, at macroscopic scales, the flow fields in the bulk and at the interface, one can capture the properties of the interfacial transport at microscopic scales (Figs. 1–13). This information is particularly important in systems where experimental parameters can be difficult to control and where data are a challenge to obtain. They include, for instance, supernovae in astrophysics, laser ablation in plasma physics, plumes in the solar atmosphere, as well as convection in planetary interiors, and industrial applications (Abarzhi *et al.,* 2013). By directly connecting the macroscopic observations far from the interface to the microscopic properties at the interface, one can further get a better understanding of connection between the Lagrangian and Eulerian dynamics away from equilibrium and potentially achieve a better control of complex transport processes involving phase transitions (Abarzhi *et al.,* 2019).

Our approach can be employed in a broad range of problems with interfacial dynamics. Particularly, our theoretical framework can be expanded to other multiphase flows, such as ablative Rayleigh–Taylor instabilities, detonation to deflagration transition, convection in stellar and planetary interiors, and stability of the liquid–liquid and liquid–vapor interface (Abarzhi and Goddard, 2019; Abarzhi *et al.,* 2019; Ilyin *et al.*, 2019; 2018; Remington *et al.,* 2019; Zhakhovsky *et al.*, 2019; 2014; Schlimpert *et al.,* 2016; Hurricane *et al.,* 2014; Anisimov *et al.,* 2013; Abarzhi, 2010; Gorokhovski and Herrmann, 2008; Meshkov, 2006; Bell *et al.,* 2004a; 2003b; Class *et al.,* 2003a; 2003b; Stein and Norlund, 2000; Peters, 2000; Clavin, 2000; Klimenko and Class, 2000; Piriz *et al.,* 1997; Sanz, 1994; Prosperetti and Plesset, 1984; Sivashinsky, 1983; and Williams, 1965). When augmented with linear response theory, it can also serve to analyze dynamic processes in fluids experiencing phase transitions and super-critical fluids (Kadanoff, 2000; Landau and Lifshitz, 1987; and Kadanoff *et al.,* 1967). Our approach can self-consistently account for the effects of compressibility, viscosity, heat conduction, mass ablation, and mass diffusion and thus enable a unified treatment of the interfacial dynamics in realistic environments (Landau and Lifshitz, 1987). We address these studies to the future.

## X. CONCLUSION

To conclude, we have analyzed from a far field the dynamics of the interface separating incompressible ideal fluids with different densities and with mass flux across it. We have found fundamental solutions of the dynamics and directly linked the interface stability to the flow fields’ structure, associating rigorous mathematical attributes to physical observables. Our approach is consistent with the spirit of Landau’s approach and identifies the new properties of the interface evolution and the new mechanisms of the interface stabilization and destabilization. Our molecular dynamics simulations illustrate the complexity of chemically reactive systems and suggest possible chemistry-induced instabilities.

## AUTHORS’ CONTRIBUTIONS

The authors contributions to this work are as follows: D.V.I.—analysis, methodology, and investigations for the general theoretical framework and the reactive molecular dynamics simulations, writing—original draft, and figures’ preparation; W.A.G.—analysis, methodology, and investigations for the reactive molecular dynamics simulations and chemistry-induced instabilities, student’s supervision, and writing—original draft; S.I.A.—conceptualization, analysis, methodology, and investigations for the general theoretical framework, writing—original draft, and writing—review and editing.

## DATA AVAILABILITY

The methods, the results, and the data presented in this work are freely available to the readers and on request from the authors.

## ACKNOWLEDGMENTS

This work is dedicated to the memory of Leo P. Kadanoff.

The authors thank for support the University of Western Australia (Australia), Award/Contract Number 10101047 [S.I.A.]; the National Science Foundation (USA), Award/Contract Number 1404449 [S.I.A.]; the Japan Society for the Promotion of Science (Japan), Award/Contract Number S11140 [S.I.A.]; the Summer Undergraduate Research Fellowship Program at California Institute of Technology (USA), the Toni and Bob Perpall Fellowship [D.V.I.]; the Department of Energy (USA), Award/Contract Number DE-SC0017710 [W.A.G.]. S.I.A. expresses her gratitude to the late Professor Leo P. Kadanoff and the late Professor Sergei I. Anisimov for discussions on singularities and stabilization mechanisms and thanks Professor Yasuhide Fukumoto for discussions on fluid instabilities. The authors express their gratitude to Professor Sivashinsky for comments and for discussions of self-similar dynamics.

## REFERENCES

*et*al., “

*et*al., “

*et*al., “

*et*al., “

*et*al., “