The dynamics of an ideal fluid is governed by the laws of conservation of mass, momentum and energy: ∂ ρ ∂ t + ∂ ∂ x i ρ v i = 0 , ∂ ∂ t ρ v i + ∂ ∂ x j ρ v i v j + ∂ P ∂ x i = 0 , ∂ E ∂ t + ∂ ∂ x i E + P v i = 0 , (1) where i = 1, 2, 3, (x ₁, x ₂, x ₃) = (x, y, z) are Cartesian coordinates, t is time, ρ is its density, v its velocity, P is the pressure field and E = ρ (e + v ²/2) is the energy density field, in which e is the specific internal energy [40]. The governing equations Eq. 1 are augmented with the closure equation of state associating pressure P and internal energy e [17]. We presume the equation of state has the form P = sρe with some constant s, describing an ideal fluid behaving similarly to an ideal gas [17].

There is zero mass flux across the interface, and the fluxes of mass, momentum and energy are conserved at that interface. Consequently, the boundary conditions at the interface are v ⋅ n = 0 , P = 0 , v ⋅ τ = a r b i t r a r y , W = a r b i t r a r y , (2) where […] is used to denote the difference in function values across the interface; the unit normal and tangent vectors at the interface are n and τ, specifically, n = ∇θ/| ∇ ∇θ| and n ⋅ τ = 0, where θ = θ (x, y, z, t) is such that θ = 0 at the interface, θ > 0 in the bulk of the heavy fluid and θ < 0 in the bulk of the light fluid; W is the specific enthalpy W = e + P/ρ.

The heavier fluid is assumed to sit above the lighter and both are subject to a body force and a time-dependent acceleration, which is directed from the heavy fluid to the light. This acceleration influences the pressure field. The acceleration direction is chosen to be that of the z-coordinate. The acceleration g = (0, 0,−g) is such that g = Gt ^a , where a is the acceleration’s power-law exponent and G > 0 is its strength. We note that G has dimensions ms ^−(a+2) and a is dimensionless. The time t > t ₀ > 0 where t ₀ is some initial instant of time.

It is assumed that the upper and lower boundaries of the domain do not influence the dynamics and also that there is an absence of mass sources. The upper and lower boundary conditions are lim z → ∞ v h = 0 , lim z → − ∞ v l = 0 . (3) where h indicates the heavy fluid and l indicates the light fluid. The governing equations are subject to a set of initial conditions which include the perturbations of the interface and/or the flow fields. The initial conditions define the characteristic length scale, time scale, and symmetry of the resulting flow.

The large-scale coherent structures in the destabilised flow are periodic arrays of bubbles and spikes in the plane normal to the acceleration direction as is set by the initial conditions. As a specific example, we consider the flow which has square symmetry [8, 41] in the plane normal to the acceleration, with wavelength λ and wavevector k = 2π/λ.

The group theory approach [6, 8, 19, 40, 41, 48, 49] accurately describes both the linear and nonlinear dynamics of the unstable interface and can be used to identify similarities and differences in the dynamics of these flows. In particularly, for linear RTI and RMI with time-varying acceleration, the linear dynamics is single-scale and is defined by the spatial period of the coherent structure, whereas the process of formation of bubbles and spikes is set by the initial conditions. For nonlinear RTI with constant and time-varying accelerations and for RMI with impulsive and time-varying accelerations, there are families of regular asymptotic solutions with the number of relevant parameters defined by the flow symmetry; the dynamics is multi-scale and characterised by the contributions of two macroscopic length scales, these being the spatial period and the amplitude of the interfacial coherent structure. Furthermore, non-equilibrium dynamics is essentially interfacial, consisting of intense fluid motion near the interface and insignificant fluid motion away from it.

Here we consider dynamics of ideal fluids, which are inviscid and are free from thermal heat flux, Eq. 1 [17]. Their velocity field is irrotational in the bulk [17, 50]. The latter can be illustrated by the analysis of slight perturbations of the fluid system in Eq. 1 near the state of hydrostatic equilibrium, similarly to [50]. Particularly: 1) by introducing a small velocity field to an initially stationary system, 2) by representing the velocity field as a superposition of a gradient of a scalar potential and a curl of a vortical potential, and 3) by considering the associated changes in the fields of density and pressure, one can linearize the system of the governing equations in the bulk Eq. 1 and the equation of state near the equilibrium state. By finding the fundamental solutions for the linear system, including their eigenvalues and eigenvectors, one can further identify the structure of the perturbation waves [50].

According to the analysis, for both compressible and incompressible ideal inviscid fluids, the vortical component of the velocity is zero in the bulk [50]. This suggests that for ideal inviscid fluids Eq. (1), the velocity field is potential in the bulk, and the vortical structures can emerge only at the interface, due to discontinuities of the tangential component of velocity and enthalpy Eq. 2. In the limit of a small viscosity, the small-scale vortical structures at the interface may be influenced by the viscosity, and the vortical structures may also appear in the bulk. However, given the interfacial character of RT/RM dynamics, the influence of the viscosity on the bulk motion may also be negligible. In this work we consider ideal incompressible fluids and their three-dimensional dynamics. We focus on RTI/RMI driven by variable accelerations with inverse quadratic power-law time-dependence and on RT-to-RM transition. We address to the future the detailed analysis of compressible and viscous effects on RT/RM dynamics.

We emphasize that even for ideal incompressible fluids, with e ̇ = 0 and ∇e = 0, the system (Eqs 1–3) is extremely challenging, requiring the solution of a system of nonlinear partial differential equations in a bulk, the solution of the boundary value problem at the nonlinear unstable interface, in addition to the solution of the relevant initial value problem, which is ill-posed [6, 8, 17, 19]. It is remarkably that even with this amount of complexity, the dynamics nevertheless has features of universality and order and can thus be approached using group theory [6, 19].

We focus on bubbles and spikes moving in the z-direction and for convenience work in a frame of reference which moves with velocity v (t) in the z-direction, where v (t) = ∂z ₀/∂t and z ₀ is the position of the bubble or spike in the laboratory frame. We define the function θ (x, y, z, t) as θ (x, y, z, t) = z* (x, y, t)−z, where z = z* (x, y, t) describes the interface.

The fluids are assumed to be incompressible and stratification and density variation are negligible. We focus on large-scale coherent dynamics having length scale λ, and presume that the length scale of shear-driven interfacial vortical structures is comparatively small, ≪ λ. The fluid motion in the bulk is potential. Consequently, the velocity of the heavy fluid v _h = ∇Φ _h and that of the light fluid v _l = ∇Φ _l and the equation of conservation of mass means that Laplace’s equation ΔΦ _h(l) = 0 is satisfied in the respective bulks, that is, θ > 0 (θ < 0).

In the moving frame of reference, the interface conditions are ρ h ∇ Φ h ⋅ n ̂ + θ ̇ | ∇ θ | = 0 = ρ l ∇ Φ l ⋅ n ̂ + θ ̇ | ∇ θ | , ρ h ∂ Φ h ∂ t + | ∇ Φ h | 2 2 + g t + d v d t z = ρ l ∂ Φ l ∂ t + | ∇ Φ l | 2 2 + g t + d v d t z , (4) and the vertical far-field boundary conditions are ∂ Φ h ∂ z | z → ∞ = − v t , ∂ Φ l ∂ z | z → − ∞ = − v t . (5)

The periodic nature of the large-scale coherent structure can be address by means of group theory [6]. Upon identification of a group which enables structurally stable dynamics, a specific Fourier series (an irreducible representation of that group) can be employed to solve the nonlinear boundary value problem (Eqs 4, 5) For three-dimensional flow with square symmetry the potentials are [41]. Φ h x , y , z , t = ∑ m , n = 0 ∞ Φ m n t cos m k x cos n k y e − α m n k z α m n k + z + f h t , Φ l x , y , z , t = ∑ m , n = 0 ∞ Φ ̃ m n t cos m k x cos n k y e α m n k z α m n k − z + f l t , (6) where α m n = m 2 + n 2 , m and n are integers, k = 2π/λ is the wavenumber, Φ _mn and Φ ̃ m n are the Fourier amplitudes for the heavy and light fluids respectively, with Φ₀₀ ≡ 0, Φ ̃ 00 ≡ 0 , and f _h (t) and f _l (t) are time-dependent functions. Symmetry requires that Φ _mn = Φ _nm and Φ ̃ m n = Φ ̃ n m .

In the vicinity of the tip of a bubble or spike we make expansions in terms of the spatial coordinates. This reduces the set of governing equations to a system of oridonary differential equations (ODEs) in terms of the interface variables and the moments [6–8, 19, 40, 41].

The local behaviour of the interfacial dynamics in the vicinity of the tip of a bubble or spike can be investigated by expanding the interface function in a power series about (x, y) = (0, 0), this being z * x , y , t = ∑ N = 1 ∞ ∑ i + j = 1 N ζ i j t x 2 i y 2 j , (7) where ζ _ij (t) = ζ _ji (t) due to symmetry, ζ (t) = ζ ₁₀(t) is the principal curvature at the tip, and N is the order of the approximation. To lowest order (that is, N = 1), the interface is z* (x, y, t) = ζ (t) (x ² + y ²). Note that ζ (t) < 0, v (t) > 0 for bubbles, and ζ (t) > 0, v (t) < 0 for spikes.

The Fourier series and interface functions are substituted into the boundary conditions at the interface and the these expressions are then expanded in power series, yielding an infinite system of ordinary differential equations for Φ _mn (t), Φ ̃ m n ( t ) and ζ _ij (t).

We introduce the following moments for the heavy and light fluids: M a , b , c t = ∑ m n Φ m n t m k a n k b α m n k c , M ̃ a , b , c t = ∑ m n Φ ̃ m n t m k a n k b α m n k c (8) These are weighted sums of the (infinite number of) Fourier amplitudes. We note that M _{a, b, c} = M _{b, a, c} and M _{a+2, b, c} + M _{a, b+2, c} = M _{a, b, c+2}, by symmetry, and similarly for M ̃ . At N = 1, the series is abbreviated to second order in x, y and first order in z, since z* (x, y, t) is quadratic in both x and y.

We also define a shear function at the interface. Traditionally, in viscous fluids, such as boundary layers, the velocity share is quantified with the shear rate, which is the rate of change of velocity at which one layer of fluid passes over an adjacent layer [51]. The velocity change is caused by, e.g., the fluid’s viscosity, and the velocity shear rate is estimated in the direction normal to the direction of motion [17]. This traditional definition has a limited applicability in our problem, because in RT/RM dynamics the shear is produced by the finite jump(s) of tangential component of velocity (and enthalpy) at the interface Eq. 2. In RT/RM dynamics, the rate of change of the tangential component of velocity (and the enthalpy) is infinite at the interface in the direction normal to the interface, since the jump of the tangential component of velocity (and the enthalpy) is finite at the interface, and since the interface is a contact discontinuity separating inviscid fluids (Eqs 1–3).

We notice however that in RT/RM dynamics the jump of the tangential velocity is non-uniform along the interface and it is a function changing continuously along the interface. The jump of the tangential component of velocity and, hence, the interfacial shear is exactly zero at the very tip of the bubble and the very tip of the spike, and it achieves its extreme values in-between the bubble and the spike. In experiments and simulations (with a finite viscosity) this is observed as the formation of shear-driven vortical structures on the side of the bubble and the side of the spike and in-between of (rather than at) the tip of the bubble and the tip of the spike [20, 23, 29, 30, 36, 45, 52–54].

Hence, in the RT/RM problem, we consider the jump of the tangential component of velocity across the interface in the direction normal to the interface as a function changing continuously along the interface. We define the shear Γ to be the spatial derivative of this function along the interface. Note that the shear Γ (hereafter—the shear function or the shear) depends on the problem parameters (including the acceleration parameters G, a, the initial growth-rate v ₀, the Atwood number A), and it also depends on time t and on the position along the interface and thus on the interface morphology, i.e., the curvature ζ (t) of the bubble/spike.

At N = 1 the shear function is Γ = Γ _x(y) with Γ _x(y) = ∂(v _x(y))/∂x (y), where [ v x ( y ) ] is the jump of the x(y) component of the velocity at the interface, with [ v x ] = x ( M ̃ 1 − M 1 ) and [ v y ] = y ( M ̃ 1 − M 1 ) , leading to the shear function Γ = M ̃ 1 − M 1 in the vicinity of the tip.

The boundary conditions at the interface become ζ ̇ = 4 M 1 ζ + M 2 2 , ζ ̇ = 4 M ̃ 1 ζ − M ̃ 2 2 , (9) 1 + A M ̇ 1 2 + ζ M ̇ 0 − M 1 2 2 − ζ g = 1 − A M 1 ̃ ̇ 2 − ζ M ̃ 0 ̇ − M ̃ 1 2 2 − ζ g . (10) M 1 − M ̃ 1 = a r b i t r a r y , (11) and the vertical far-field requirement becomes M 0 = − M ̃ 0 = − v t (12) where M ₀ = M _{0, 0, 0}, M ₁ = M _{2, 0, −1} and M ₂ = M _{2, 0, 0}, and the ratio of the densities of the fluids is parametrised by the Atwood number A = (ρ _h − ρ _l )/(ρ _h + ρ _l ), and 0 ≤ A ≤ 1.

Initial conditions for the dynamical system (Eqs 9–12) include the conditions for the initial curvature ζ (t ₀) and the initial velocity v (t ₀) at some initial instance of time t ₀. The initial velocity sets the initial growth-rate v ₀ = |v (t ₀)|.

Representation in terms of the moments M , M ̃ and interface variable ζ allows us to accommodate the nonlocal nature of RT/RM dynamics, and also enables us to study the interplay of the harmonics and to derive solutions in both the linear and the nonlinear regimes.

Equations 1–5 and Equations 9–12 are in the dimensional form, hence requiring us to identify the problem’s length and time scales. The length scale is defined by the wavelength λ with the corresponding wavevector k. In this work we focus of the case a =−2, that is g = G/t ². The dimension of the acceleration strength G is dim (G) = dim (1/k) = dim (λ) = m. The dimension of the initial growth-rate v ₀ is   dim (v ₀) = m/s. Hence, at a =−2, there are two natural time scales in the problem: τ _G = G/v ₀ and τ ₀ = 1/(kv ₀). Their interplay is depends on the parameter (Gk), with τ _G /τ ₀ ≪ (≫)1 for (Gk) ≪ (≫)1. Here we use τ ₀ for time-scale and consider the broad range of values (Gk) ∈ (0, ∞). We also presume that t ₀ ≫ (τ _G , τ ₀) in order to preempt the effect of the choice of time origin on the dynamics.

When t−t ₀ ≪ τ ₀, the perturbation amplitude is small and the interface is approximately flat, meaning that only first order harmonics need be retained in moments, that is, M ₀ = 2Φ₁₀, M ̃ 0 = 2 Φ ̃ 10 ; M _n = k ⁿ Φ₁₀, M ̃ n = k n Φ ̃ 10 , n = 1, 2. Continuity of the normal components of velocity and momentum at the fluid interface, that is Eqs 9 and 10, become d ζ ̂ d t + k 4 v = 0 , d v d t + A k 2 v 2 + 4 A G t 2 ζ ̂ = 0 (13) where ζ ̂ = ζ / k is a nondimensional curvature. These apply for any sign of kζ ₀ and v (t ₀)/v ₀, and can be combined to give d 2 ζ ̂ d t 2 = 2 A d ζ ̂ d t 2 + A k G t 2 ζ ̂ .

For Gk ≫ 1, the dynamics is acceleration-driven and consequently the relative contributions of the terms is | A k G ζ ̂ / t 2 | ≫ | ( 2 A d ζ ̂ / d t ) 2 | . The system becomes d 2 ζ ̂ d t 2 = A k G t 2 ζ ̂ , v = − 4 k d ζ ̂ d t . (14) The curvature equation has solution ζ ̂ = C 1 t τ p + C 2 t τ − p , p = A k G , (15) where C ₁ and C ₂ are integration constants which can be determined from the initial conditions ζ (t ₀) and v (t ₀), with |ζ (t ₀)/k|≪ 1 and |τkv (t ₀)|≪ 1. This dynamics is driven by the acceleration and is of RT type. [2, 8, 17, 49].

For kG ≪ 1, the dynamics is driven by the initial conditions and consequently the relative contributions of the terms is | A k G ζ ̂ / t 2 | ≪ | ( 2 A d ζ ̂ / d t ) 2 | . The system becomes d 2 ζ ̂ d t 2 = 2 A d ζ ̂ d t 2 , v = − 4 k d ζ ̂ d t . (16) The curvature equation has solution, with integration constants C ₁ and C ₂, ζ ̂ = − 1 2 A ln C 1 t + C 2 . This dynamics is driven by the initial growth-rate and is of RM type.

The very-early-time (t ∼ t ₀) dynamics yields ζ − ζ 0 k = − 1 4 t − t 0 τ 0 sgn v t 0 v 0 , v − v 0 = − A 2 v 0 t − t 0 τ 0 + 4 A G k v 0 t 0 2 ζ 0 k t − t 0 τ 0 . (17) suggesting that the formation of RT/RM bubbles and spikes are defined by the initial conditions [8, 49].

For t ≫ τ, we must retain higher order harmonics in the moments. Asymptotic solutions of the resulting set of equations can be derived. To leading order, these will have the following time-dependence: ζ k ∼ c o n s t , v ∼ 1 k τ 0 t τ 0 b , Φ n , Φ ̃ n ∼ 1 k τ 0 t τ 0 a 2 , M n , M ̃ n ∼ k n Φ n , k n Φ ̃ n , (18) where b is a constant to be determined. For N = 1, the first two harmonics are retained and we arrive at a one-parameter family of solutions. We choose the nondimensional curvature ζ ̂ = ζ / k to parametrise the family. The velocity is v t = ± 9 − 64 ζ ̂ 2 3 + 10 A ζ ̂ − 128 A ζ ̂ 3 + q 2 k t 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ , (19) where q = 3 + 10 A ζ ̂ − 128 A ζ ̂ 3 2 − 8 A k G ζ ̂ 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ , the shear function is Γ t = ± 3 3 + 10 A ζ ̂ − 128 A ζ ̂ 3 + q t 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ , (20) and the corresponding Fourier amplitudes are Φ 10 = − 2 + 8 ζ ̂ 3 + 8 ζ ̂ v , Φ 20 = 1 + 8 ζ ̂ 6 + 16 ζ ̂ v , Φ ̃ 10 = 2 − 8 ζ ̂ 3 − 8 ζ ̂ v , Φ ̃ 20 = − 1 − 8 ζ ̂ 6 − 16 ζ ̂ v . (21) In (Equations 19–21), the positive sign applies for bubbles and the negative sign applies for spikes.

In (Equations 19–21) the limit k G → 0 , v t → ± 9 − 64 ζ ̂ 2 3 + 10 A ζ ̂ − 128 A ζ ̂ 3 k t 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ , Γ t → ± 6 3 + 10 A ζ ̂ − 128 A ζ ̂ 3 t 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ (22) is consistent with RM dynamics induced by the acceleration g = Gt ^a with a →−2⁻, whereas the limit k G → ∞ , v t → ± 9 − 64 ζ ̂ 2 k t − 2 A k G ζ ̂ 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ , Γ t → ± 6 t − 2 A k G ζ ̂ 64 A ζ ̂ 2 + 9 A − 48 ζ ̂ (23) is consistent with RT dynamics induced by the acceleration g = Gt ^a with a →−2⁺ [8, 46].

To conveniently describe properties of special solutions for nonlinear bubbles, we introduce dimensionless variables ζ ˇ = − ζ k , υ ˇ = A k t 3 v t , Γ ˇ = A t 2 Γ t , and use superscript 0 to refer to the limit kG → 0 and superscript ∞ to refer to the limit kG → ∞.

To conveniently describe special solutions for nonlinear spikes, we introduce the dimensionless and negated variables ζ ̂ = ζ k , v ̂ = − A k t 3 v t , Γ ̂ = − A t 2 Γ t and use superscript 0 to refer to the limit kG → 0 and superscript ∞ to refer to the limit kG → ∞.

In classical RTI/RMI, RT bubbles and spikes move significantly faster than RM bubbles and spikes in either linear or nonlinear regime. Hence, experiments and simulations can differentiate between RT dynamics and RM dynamics by measuring the growth-rate of RT/RM bubbles/spikes [18, 20, 21, 27, 32, 33, 45, 52–54]. For variable acceleration g = G/t ² considered in the present work, RT/RM bubbles and spikes move slowly in the linear regime, and have nonlinear asymptotic solutions changing with time as ∼ 1 / t in the nonlinear regime. To accurately identify a pre-factor of a power-law, a substantially span of temporal and spatial scales is required, which is usually very challenging to achieve in experiments and simulations [45, 56]. An accurate implementation of the acceleration with the magnitude g = G/t ² is also challenging. Nevertheless, since our analysis finds that, expect for the time-dependence, the properties of RT/RM dynamics of bubbles/spikes are similar in many regards to those of classical RTI/RMI, we may compare our analysis with existing experiments and simulations. Our theoretical results are in agreement with available observations [18, 20, 21, 27, 32, 33, 45, 52–54].

Particularly, our theory accurately reproduces the interfaciality of RT/RM dynamics, which is observed in experiments and simulations. This includes the velocity field with intense motion of the fluids in a vicinity of the interface, with effectively no motion away from the interface, and with vortical structures produced by shear at the interface [20–22, 27, 45]. Our analysis further finds that in the nonlinear regime RT bubbles are curved and RM bubbles are flat. The flattening of RM bubble front and the curved shape of RT bubble are distinct feature of nonlinear RMI and RTI, respectively; they are observed in experiments and simulations in fluids and plasmas [18, 20, 21, 27, 32, 33, 45, 52–54]. Our analysis also finds that, at any density ratio and acceleration strength, nonlinear RT/RM spikes can achieve a finite curvature magnitude and can move with a time-dependence faster than that of nonlinear RT/RM bubbles, in agreement with observations [18, 20, 21, 27, 32, 33, 45, 52–54].