Rayleigh**–**Taylor instability (RTI) has critical importance for a broad range of plasma processes, from supernovae to fusion. In most instances, RTI is driven by variable acceleration, whereas the bulk of existing studies have considered only constant and impulsive acceleration. This work focuses on RTI driven by acceleration with a power-law time-dependence. We review the existing theoretical approaches, apply the group theory approach to solve this long-standing problem, and yield the unified framework for the scale-dependent dynamics of Rayleigh**–**Taylor (RT) bubbles and RT spikes. For the early-time linear dynamics, we provide the dependence of RTI evolution on the acceleration parameters and the initial conditions. For the late-time nonlinear dynamics, we find a continuous family of asymptotic solutions, directly link the interface velocity to the interface morphology and the interfacial shear, derive solutions for the regular bubbles and for the singular spikes, and study the stability of these solutions. The properties of special nonlinear solutions in the RT family are scrupulously described, including the critical, Taylor, Layzer-drag, and Atwood solutions. It is shown that the fastest Atwood bubble is regular and stable, and the fastest Atwood spike is singular and unstable. The essentially multi-scale and interfacial character of RT dynamics is demonstrated. The former can be understood by viewing the RT coherent structure of bubbles and spikes as a standing wave with the growing amplitude. The latter implies that RT flow has effectively no motion of the fluids away from the interface and has intense motion of the fluids near the interface, with shear-driven vortical structures appearing at the interface. Our theory agrees with available observations and elaborates extensive benchmarks for future research and for better understanding of RT driven phenomena in plasmas.

## I. INTRODUCTION

Rayleigh–Taylor instability (RTI) has critical importance for a broad range of processes in nature and technology at astrophysical and molecular scales.^{1–4} Examples of plasma processes governed by RTI include the fingering of interstellar media along the edges of black holes, the abundance of chemical elements in supernova remnants, the effect of plasma irregularities in the Earth ionosphere on the regional climate change, the formation of hot spots in inertial confinement fusion, and the efficiency of plasma thrusters.^{1–9} While realistic plasmas are necessarily electro-dynamic, with magnetic fields and charged particles involved, the grip and control of Rayleigh**–**Taylor (RT) plasma phenomena often require the better understanding of hydrodynamic aspects of non-equilibrium dynamics and the development of rigorous yet practical theory of RTI in neutral plasmas (fluids) at continuous (hydrodynamics) scales.^{1,2} In this work, we study the long-standing problem of Rayleigh**–**Taylor instability driven by acceleration varying as the power-law with time, and we employ the group theory approach to solve the boundary value problem for the scale-dependent RT dynamics in the early-time linear and late-time nonlinear regimes.^{1–4,10–12}

This paper is organized as follows: We start with the brief introduction in Sec. I. We provide the method in Sec. II, and we overview the governing equations (II A), the theoretical approaches (II B), and the group theory approach (II C), including the space groups, the expansions, and the dynamical system. We present the results for the scale-dependent Rayleigh**–**Taylor instability driven by variable acceleration in Sec. III. This includes the early-time linear dynamics (III A), the late-time nonlinear dynamics (III B), the focused analysis of the solutions for nonlinear RT bubbles (III C) and for nonlinear RT spikes (III D), and the properties of RT dynamics with variable acceleration (III E). We finalize the work with discussion and conclusion in Sec. IV.

Rayleigh**–**Taylor instability develops when fluids of different densities are accelerated against their density gradients and/or when the fluid interface moves with an acceleration directed toward the denser fluid and normal to the nearly planar interface.^{3,4} In the vastly different physical circumstances, at astrophysical and molecular scales, RT flows have a number of similar features of their evolution.^{11,12} RTI starts to develop when the flow fields and the interface are slightly perturbed near the equilibrium state.^{3,4} Small perturbations grow quickly (e.g., exponentially in time), and the flow transits to a nonlinear stage, where the growth rate slows (to a power-law in time). The interface is transformed to a composition of a large-scale structure of bubbles and spikes and small-scale vortical structures.^{11,12} The bubbles (spikes) are portions of the light (heavy) fluid penetrating the heavy (light) fluid. Their large-scale dynamics is coherent—as is set by the initial conditions. The small-scale vortical structures are produced by shear at the interface; their dynamics can be irregular. Due to the enhanced interactions and the coupling of scales, RT dynamics becomes complicated. Intense interfacial mixing ensues with time.^{11–21} The dynamics of Rayleigh–Taylor mixing is believed to be self-similar. Particularly, for a constant driving acceleration, the flow length scale in the acceleration direction grows quadratic in time.^{11,14–21}

RTI and RT mixing are extremely challenging to study in theory, experiments, and simulations.^{1,2} RT experiments use advanced technologies to meet tight requirements for the flow implementation, diagnostics, and control.^{13–18} RT simulations employ highly accurate numerical methods and massive computations to track unstable interfaces, capture small-scale processes, and enable a large span of scales.^{19–21} In theory, we have to reliably grasp non-equilibrium RT dynamics, identify universal properties of the asymptotic solutions, and capture the symmetries of the unstable dynamics.^{11,12} Significant success was achieved recently in the understanding of RTI and RT mixing with constant acceleration.^{1,2} In particular, the group theory approach^{11,12} found the multi-scale character of nonlinear RTI and the order in RT mixing and explained a broad set of experiments in plasmas and fluids, which observed that RT mixing may depart from the scenario of canonical turbulence and may keep order at high Reynolds numbers up to 3.2 × 10^{6}.^{13–18}

In realistic circumstances, RT relevant phenomena are often driven by variable acceleration:^{1,2,5–18,22} In supernova blast, the acceleration is induced by strong variable shocks, and RT mixing with variable acceleration enables the conditions for the synthesis of heavy and intermediate mass elements. In inertial confinement fusion, the shocks' induced accelerations and decelerations of the ablation front influence the formation of hot spots. In nano-fabrication, the transformation of materials at a high strain rate is governed by RTI with variable acceleration.

In this work, we focus on RTI subject to variable acceleration, and we consider accelerations with a power-law time-dependence.^{1,22,23} On the side of fundamentals, power-law functions are important to study since they may result in new invariant and scaling properties of the dynamics.^{10,22} As regards to applications, power-law functions can be used to adjust the acceleration's time-dependence in realistic environments and to ensure the practicality of theoretical results.^{5–9,22} To solve the long-standing problem of the scale-dependent early-time and late-time RT dynamics with variable acceleration, we apply the group theory approach.^{11,12,22} This approach worked remarkably well for constant and impulsive accelerations.^{11–14} In the case of variable acceleration, it revealed a new mechanism for energy accumulation and transport at small scales in supernovae.^{22} Furthermore, the group theory approach found the dependence of the dynamics on the exponent of the acceleration's power-law.^{22} Particularly, the scale-dependent interfacial dynamics is of Rayleigh–Taylor type and is driven by the acceleration for the acceleration exponents greater than (−2). It is of Richtmyer–Meshkov (RM) type and is driven by the initial growth rate for the acceleration exponents smaller than (−2).^{22,23}

Here, we consider RTI for a three-dimensional spatially extended periodic flow with hexagonal symmetry in the plane normal to the acceleration, and we solve the boundary value problem involving boundary conditions at the interface and at the outside boundaries of the domain.^{10–12,22} For the early-time linear dynamics, we provide the dependence of the RTI growth rate on the acceleration's parameters. For the late-time non-linear dynamics, we find a continuous family of asymptotic solutions, directly link the interface dynamics to the interface morphology and the interfacial shear, and derive the solutions for the regular bubbles and for the singular spikes. The family contains the special solutions for the bubbles and for the spikes, including the critical, convergence limit, Taylor, Layer-drag, Atwood, and flat solutions. The fastest Atwood bubble is regular and stable, whereas the fastest Atwood spike is singular and unstable. We also reveal the interfacial and multi-scale character of the scale-dependent RT dynamics and analyze the mechanism of transition from the scale-dependent dynamics to the self-similar mixing for RT bubbles and for RT spikes. Our theory agrees qualitatively with available observations, and elaborates extensive quantitative benchmarks for future research.^{1,2,13–21} The properties of the scale-dependent RT dynamics found by our theory—the fields of the velocity and pressure, the interplay of the interface morphology and shear with the flow acceleration, the interface growth, and the growth rate—serve for better understanding of RT driven plasma phenomena in nature and technology.^{1,2,5–10,17,18}

## II. METHOD

### A. Governing equations

RT dynamics of ideal fluids is governed by the conservation of mass, momentum, and energy,

with spatial coordinates $(x1,x2,x3)=(x,y,z)$, time $t$, fields of density, velocity, pressure, and density energy $(\rho ,v,P,E)$, $E=\rho (e+v2/2)$, and $e$ the specific internal energy.^{10,22}

We introduce a continuous local scalar function $\theta (x,y,z,t)$, whose derivatives $\theta \u0307$ and $\u2207\theta $ exist (the dot marks a partial time-derivative). The function $\theta $ is $\theta =0$ at the interface. The heavy (light) fluid is located in the region with $\theta >0$ ($\theta <0$), and its fields are $(\rho ,v,P,E)\u2192(\rho ,v,P,E)h(l)$ and are marked with the subscript $h(l)$.^{10–12,22,23} The ratio of the fluids' densities is parameterized by the Atwood number $A=(\rho h\u2212\rho l)/(\rho h+\rho l)$.^{3,4,10–23}

By using the Heaviside step-function $H(\theta )$, we represent the flow fields as $(\rho ,v,P,E)=(\rho ,v,P,E)hH(\theta )+\u2009(\rho ,v,P,E)lH(\u2212\theta )$.^{11,12,22} At the interface, the balance of fluxes of mass and normal and tangential components of momentum and energy obeys the boundary conditions,

where the jump of a quantity across the interface is denoted by $[\u2026]$; the normal and tangential unit vectors of the interface are $n$ and $\tau $ with $n=\u2207\theta /|\u2207\theta |$ and $(n\u22c5\tau )=0$; the mass flux across the moving interface is $j\u0303=\rho (n\theta \u0307/|\u2207\theta |+v)$; and the specific enthalpy is $W=e+P/\rho $.

The boundary conditions [Eq. (2)] are derived from the governing equations (1), assuming that the mass flux is conserved at the interface, $[j\u0303\u22c5n]=[j\u0303n]=0$, with $j\u0303n=j\u0303\u22c5n$. There is an important particular case when the conserved mass flux is zero at the interface, $j\u0303n|\theta =0\xb1=0$. This case corresponds to a so-called contact discontinuity.^{10–12,22–24} The condition $j\u0303n|\theta =0\xb1=0$ leads to the continuity of the normal component of velocity at the interface $[v\u22c5n]=0$ and transforms the condition for the conservation of the normal component of momentum at the interface to the continuity of pressure at the interface, $[P]=0$. For the zero mass flux at the interface, $j\u0303n|\theta =0\xb1=0$, the condition of continuity of the tangential component of momentum at the interface holds true for an arbitrary jump of the tangential velocity at the interface, $[v\u22c5\tau ]=arbitrary$, and the condition of continuity of the energy flux at the interface holds true for any jump of the enthalpy at the interface, $[W]=arbitrary$. In full consistency with the classic results,^{10} in the case of zero mass flux at the interface, $j\u0303n|\theta =0\xb1=0$, the boundary conditions at the interface are

Per Eq. (3), the normal components of velocity and pressure are continuous at the interface, and the tangential velocity component and the enthalpy can be discontinuous at the interface.^{10–12,22–24}

The outside boundaries have no influence on the dynamics, and the flow is free from mass sources. This leads to

The initial conditions are the initial perturbations of the interface and the flow's fields.

RTI is driven by the acceleration (the body force or the gravity) $g$, directed from the heavy to the light fluid, $g=(0,0,\u2212g)$, $g=|g|$. In the presence of the acceleration $g$, the pressure is modified as $P\u2192P\u2212\rho gz$.^{10–12,22–24} We consider accelerations that are power-law functions of time, $g=Gta$, $t>0$. Here, $a$ is the acceleration exponent, $a\u2208(\u2212\u221e,+\u221e)$, and $G$ is the pre-factor, $G>0$. Their dimensions are $[G]=m/s2+a$ and $[a]=1$.^{10–12,22,23}

We consider RT flow periodic in the plane $(x,y)$ normal to the $z$ direction of the acceleration $g$; this periodicity is set by the initial conditions.^{10–12,22,23} For ideal incompressible fluids, the initial conditions define the length-scale, time-scale, and symmetry of the dynamics. For the length-scale, we choose the spatial period (wavelength) $\lambda $ of the initial perturbation, $[\lambda ]=m$, or, for convenience, the wavevector $k$, with $[k]=m\u22121$ and $1/k\u223c\u2009\lambda $. For the given value of the wavevector $k$ and the acceleration parameters $a,G$, there are two natural time-scales: the time-scale $\tau G=(kG)\u22121/(a+2)$ set by the acceleration and the time-scale $\tau 0=(kv0)\u22121$ set by the initial growth-rate $v0$ at the initial time $t0$.^{22,23} For the exponent values $a>\u22122$, these time-scales relate as $\tau G\u2009\u226a\u2009\tau 0$, and $\tau G$ is the smallest time-scale corresponding to the fastest process. Hence, for Rayleigh–Taylor dynamics driven by the acceleration, the relevant time-scale is $\tau G$.^{22,23} For $a<\u22122$, the time-scales relate as $\tau 0\u2009\u226a\u2009\tau G$, and for Richtmyer–Meshkov dynamics driven by the initial growth-rate, the relevant time-scale is $\tau 0$.^{22,23} Since we consider Rayleigh–Taylor dynamics in this work, we set the time-scale to be $\tau =\tau G$, $[\tau ]=s$, and we study RT evolution for $t>t0$.^{22,23}

### B. Theoretical approaches

The mathematical problem of Rayleigh–Taylor instability requires one to solve the system of nonlinear partial differential equations in four-dimensional space-time, solve the boundary value problem at the nonlinear freely evolving discontinuity and outside boundaries, and also solve the ill-posed initial value problem, with account for singularities developing in a finite time [Eqs. (1)–(4)].^{11,12,22} This challenging problem can be handled in well-defined theoretical approximations; empirical models may be applied to interpolate experimental and numerical simulation data.^{11,12,25}

Having started from the classical works of Rayleigh^{3} and Taylor,^{26} significant progress was achieved in understanding the early-time RT dynamics in plasmas and in fluids in a broad range of conditions.^{27–30} Particularly, in the linear regime, the RTI growth rate was derived in a continuous approximation and through the solution of the boundary value problem [Eqs. (1)–(4)] with account for the effects of surface tension, viscosity, mass diffusion, mass ablation, magnetic fields, compressibility, thermal conductivity, stratification, finite size of the domain, and other factors.^{3,4,11,12,27–30} For RT flows with various symmetries and with variable acceleration, the group theory approach was employed.^{11–14,18,22,23} Weakly-nonlinear stages of RTI were studied theoretically and empirically with the focus on the effect of initial spectra on the RTI growth rate.^{30–32}

Nonlinear dynamics of RTI is a long-standing problem, investigated by Davies and Taylor,^{4} Layzer,^{33} Birkhoff,^{34} and Garabedian.^{35} The group theory approach analyzed a broad class of Rayleigh–Taylor problems, including dynamics of two- and three-dimensional RT flows with various symmetries.^{11,12,37} In the case of constant acceleration, it reconciled with one another the approaches^{4,33–37} and identified the properties of nonlinear RT dynamics, including the tendency to maintain isotropy in the plane normal to the acceleration, the discontinuity of the dimensional crossover, and the multi-scale character of the dynamics of nonlinear RT bubbles.^{11,12}

The interactions of multiple scales further trigger the transition of the scale-dependent RT dynamics to the self-similar mixing.^{13–18} The transition may occur via the growth of horizontal scales of the RT coherent structure, as found by the empirical models, including interpolation, buoyancy-drag, and bubble merge models.^{39–41,47} It may also occur due to the dominance of the vertical scale, as first found by the group theory approach through the momentum model.^{11,42}

By further assuming that the self-similar RT mixing is similar to canonical turbulence due to its large Reynolds number, the interpolation and turbulence models were developed; they substituted the growth rate of RT mixing in turbulent scaling laws and computed turbulence effects.^{25,43–47} By being harmonious with the theory of canonical turbulence^{10,48} and with models^{39–41,43–47} and by analyzing symmetries and invariants of RT dynamics, the group theory approach identified, through the momentum model, the properties of RT mixing not considered in other studies.^{11,13,25,42} Particularly, it found that in the case of constant acceleration, the self-similar RT mixing may exhibit order since it has stronger correlations, steeper spectra, weaker fluctuations, and stronger sensitivity to deterministic conditions when compared to canonical turbulence.^{11,12} These results explained the high Reynolds number experiments in plasmas and fluids.^{13–15,17,18}

### C. Group theory approach for Rayleigh–Taylor dynamics

In this work, we apply the group theory approach to yield the unified framework for the scale-dependent RT dynamics with variable acceleration in the early-time linear and late-time nonlinear regimes.^{11–14,22,23} We solve the governing equations (1)–(4) by employing the theory of space groups, by using canonical forms of the Fourier series and spatial expansions, and by accounting for non-locality and singularities of the interface dynamics.^{12,37,38,49} For the early-time regime, we provide the dependence of RT evolution on the acceleration parameters and initial conditions. For the late-time regime, we find nonlinear solutions for regular RT bubbles and for singular RT spikes and provide physics interpretations of mathematical attributes of the boundary value problem [Eqs. (1)–(4)]. By analyzing the properties of regular RT bubbles and singular RT spikes, we further outline the mechanism of the flow transition to the self-similar mixing. We address to the future the establishment of a direct link of the scale-dependent RT dynamics to the self-similar RT mixing within the group theory approach.

#### 1. Space groups

For spatially periodic flows, space groups can be applied since RT dynamics is invariant with respect to a group $G$ whose generators are translations in the plane, rotations, and reflections.^{11,12,49} These groups are also known as Fedorov and/or Schoenflies groups, and they are commonly used in physics.^{10,49} While there are 7 one-dimensional and 17 two-dimensional space groups, only some of these groups should be considered. Particularly, groups relevant to structurally stable RT dynamics must have the anisotropy in the acceleration direction and the inversion in the normal plane, such as groups of hexagon $p6mm$ or square $p4mm$ for a three-dimensional (3D) flow and group $pm11$ in a two-dimensional (2D) flow.^{10–12,49} We use international classification of the space groups, where $p$ stands for periodicity in one (two) direction(s), and for each of the spatial directions, $1$ stands for the unit element, $m$ stands for the mirror plane of reflection, and $6(4)$ stands for the $6(4)$-fold axis of rotation.^{10,49}

By using the techniques of group theory, we apply irreducible representations of a relevant group to expand the flow fields as Fourier series. We further make spatial expansions in a vicinity of a regular point at the interface (i.e., the tip of the bubble or the tip of the spike). The governing equations are then reduced to a dynamical system in terms of surface variables and moments. The system's solutions are sought, and their stability is studied.^{11,12,37,38} More details on the group theory approach can be found in works on RTI with constant acceleration.^{11,12,14,22,37,38}

#### 2. Fourier series and spatial expansions

We consider incompressible dynamics with negligible stratification and density variations. For incompressible ideal fluids, the spatial period (wavelength) $\lambda $ of the initial perturbation defines the length-scale of the early-time linear and late-time nonlinear dynamics. We focus on the large-scale coherent dynamics with length scale $\u223c\lambda $, presuming that shear-driven interfacial vortical structures are small with length scales $\u226a\lambda $.^{11,12}

For the large-scale coherent structure, the fluid motion is potential in the bulk, and the velocity of the heavy (light) fluid is $vh(l)=\u2207\Phi h(l)$. For incompressible ideal fluids, the equation for the conservation of mass in Eq. (1) leads to the Laplace equation $\Delta \Phi h(l)=0$ for $\theta >0\u2009(<0)$.^{10–12} At the outside boundaries of the domain, the velocity of the heavy (light) fluid vanishes, $\u2207\Phi h(l)|z\u2192+\u221e(\u2212\u221e)=(0,0,0)$ [Eq. (4)]. At the regular point of the interface—the tip of the bubble (spike)—with the coordinate $(0,0,z0(t))$, the velocity of the heavy (light) fluid is $\u2207\Phi h(l)|(0,0,z0)=(0,0,v(t))$, where $v(t)=\u2202z0/\u2202t$.

We consider a three-dimensional coherent structure with hexagonal symmetry in the plane normal to the acceleration (Fig. 1).^{10–12,22,23,49} In order to make Fourier series expansions of the flow fields, we recall that for group $p6mm$ with a hexagonal lattice, the spatial periods $a1(2)$ have equal lengths, $|a1(2)|=\lambda $, and are inclined relative to one another at the angle $2\pi /3$, with $a1\u22c5a2=\u22121/2$. The wave-vectors $k1(2)$ of the reciprocal lattice are defined by the relations $k1(2)\u22c5a2(1)=0$ and $k1(2)\u22c5a1(2)=2\pi $. The wave-vectors $k1(2)$ are inclined relative to one another at the angle $\pi /3$, with $k1\u22c5k2=1/2$, and have an equal length $k=4\pi /(\lambda 3)$. A linear combination of independent spatial periods $a1(2)$ defines the spatial period $a3=\u2212(a1+a2)$, and a linear combination of independent wave-vectors of the reciprocal lattice $k1(2)$ defines the wavevector $k3=k1\u2212k2$.^{11,37,49}

For convenience, we perform the derivations in the non-inertial frame of reference moving with the velocity $v(t)$ in the $z$-direction relative to the laboratory frame of Refs. 11, 12, 23, and 37. Then, the Fourier series expansions of the potentials $\Phi h(l)$ have the following form:

Here, $r=(x,y)$ and $\Phi m(\Phi \u0303m)$ are the Fourier amplitudes of the heavy (light) fluid, $fh(l)$ are time-dependent functions, and $m$ is the natural number.^{11,12,22,37} Cross terms appear in high orders. The Fourier series in Eq. (5) hold true upon group $p6mm$ transformations, including rotations on angles $\xb1\pi /3,\xb12\pi /3,\xb1\pi $ around the $z$-axis and reflections $x\u2192\u2212x$, $y\u2192\u2212y$, and $(x,y)\u2192\u2212(x,y)$ in the mirror planes along the $x$- and $y$-axes. The Fourier series in Eq. (5) ensure that in the moving non-inertial frame of reference, the tip of the bubble (spike) with coordinates $(0,0,0)$ is the stagnation point with $\u2207\Phi h(l)|(0,0,0)=(0,0,0)$ and at the outside boundaries $\u2207\Phi h(l)|z\u2192+\u221e(\u2212\u221e)=(0,0,\u2212v(t))$.

In the laboratory frame of reference near the tip of the bubble (spike), the interface has the form $Z=z0(t)+z*(r,t)$, and the interfacial function $\theta $ is $\theta =\u2212z+Z=\u2212z+z0(t)+z*(r,t)$. In the moving non-inertial frame of reference near the tip of the bubble (spike), the interface is described as $Z\xaf=z*(r,t)$ and the function $\theta $ is $\theta =\u2212z+Z\xaf=\u2212z+z*(r,t)$, where the function $z*(r,t)$ is defined locally as

Here, $N$ is the natural number, $\zeta N$ are the surface variables, and $\zeta 1=\zeta $ is the principal curvature at the bubble (spike) tip. Cross terms appear in high orders.^{11,12,22,23,37}

#### 3. Dynamical system

Upon substituting these expansions in the governing equations and further expanding the equations in the vicinity of the regular point of the interface—the tip of the bubble (spike)—we derive from the governing equations a dynamical system in terms of moments and surface variables.^{11,12,22,23,37} For group $p6mm$, at $N=1$, the interface is $z*=\zeta (x2+y2)$, and the dynamical system is

These equations represent the continuity of the normal component of velocity and the normal component of momentum at the interface, the discontinuity of the tangential component of velocity at the interface, and the absence of mass sources. The moments $Mn(M\u0303n)$ are

Since the tangential component of the velocity can have a jump at the interface, there can be the interfacial shear [Eq. (3)]. In order to quantify the interfacial shear, we define the shear function as the directional derivative of the jump of the tangential component of velocity at the interface in the direction of the tangential vector at the interface. In the vicinity of the regular point—the tip of the bubble (spike), at $N=1$, the shear function is $\Gamma =\Gamma x(y)$, with $\Gamma x(y)=\u2202\u230avx(y)\u230b/\u2202x(y)$ and $\Gamma =\u2212M1+M\u03031$ in Eq. (7).^{22} According to the boundary conditions in Eqs. (3) and (7), the shear value can be arbitrary.

In a dynamical system [Eq. (7)], the initial conditions at time $t0$ are the initial curvature $\zeta 0=\zeta (t0)$ and velocity $v(t0)$; the latter defines the initial growth-rate $v0=|v(t0)|$. As time evolves, bubbles move up and are concave down, and spikes move down and are concave up (Fig. 1). For Rayleigh–Taylor dynamics with $a>\u22122$, we set the length-scale in the system [Eq. (7)] to be $1/k$, where the wavevector $k$ is $k=4\pi /3\lambda $, and the time-scale to be $\tau =\tau G=(kG)\u22121/(a+2)$.^{22,23} We find solutions for RT bubbles and RT spikes in the early time linear regime, $(t\u2212t0)\u2009\u226a\u2009\tau $, and in the late time nonlinear regime, $(t\u2212t0)\u2009\u226b\u2009\tau $, with $t\u2009\u226b\u2009\tau $ and $t\u2009\u226b\u2009t0$.

## III. RESULTS

### A. Early-time linear dynamics

#### 1. Solutions for a linearized system

For the early-time dynamics with $(t\u2212t0)\u2009\u226a\u2009\tau $, only the first order harmonics are retained in the expressions for the momentum, $Mn=kn\Phi 1,\u2009M\u0303n=kn\Phi \u03031,\u2009n=0,1,2$ with $M0=\u2212M\u03030=\u2212v$.^{23} For Rayleigh–Taylor dynamics with $a>\u22122$, for small perturbations, the dynamical system is linearized, leading to

Its solution is

where $Ip$ is the modified Bessel function of the $p$ th order (Table I).^{22,23} For constant acceleration, $a=0$, this solution reproduces the classical result. Figure 2 illustrates the solutions in Eq. (10) for the early-time dynamics of RT bubbles and RT spikes at various values of the acceleration exponent $a$. In the linear regime, the growth and growth rates of the bubbles and the spikes are nearly symmetric.

Bubbles . | Spikes . |
---|---|

$v\u2009k\u2009\tau >0,\u2003\zeta k<0$ | $v\u2009k\u2009\tau <0,\u2003\zeta k>0$ |

$t\u2212t0\tau \u226a1,\u2003v=\u22124kddt\zeta k,\u2003\zeta k=C1t\tau \u2009I12s(As(t\tau )s)+C2t\tau \u2009I\u221212s(As(t\tau )s),\u2003s=a+22$ | |

$t\u2248t0,\u2003v\u2212v(t0)\u2248\u22124A(\zeta 0/k)(\tau k)\u22121(t0\tau )a(t\u2212t0\tau ),\u2003\zeta \u2212\zeta (t0)\u2248\u2212k4(t\u2212t0\tau )sgn(v(t0)v0)$ |

Bubbles . | Spikes . |
---|---|

$v\u2009k\u2009\tau >0,\u2003\zeta k<0$ | $v\u2009k\u2009\tau <0,\u2003\zeta k>0$ |

$t\u2212t0\tau \u226a1,\u2003v=\u22124kddt\zeta k,\u2003\zeta k=C1t\tau \u2009I12s(As(t\tau )s)+C2t\tau \u2009I\u221212s(As(t\tau )s),\u2003s=a+22$ | |

$t\u2248t0,\u2003v\u2212v(t0)\u2248\u22124A(\zeta 0/k)(\tau k)\u22121(t0\tau )a(t\u2212t0\tau ),\u2003\zeta \u2212\zeta (t0)\u2248\u2212k4(t\u2212t0\tau )sgn(v(t0)v0)$ |

Note that Richtmyer–Meshkov dynamics with $a<\u22122$ is driven by the initial growth-rate $v0$ and is independent of the acceleration; its early-time solution is

For $a>\u22122$, the contribution of these initial-growth-rate driven terms is negligible when compared to those induced by the acceleration in Eq. (10).

#### 2. Effect of initial conditions

Solutions for a linearized system [Eq. (9)] describe, as time evolves, the dynamics of bubbles with $\zeta \u22640,v\u22650$ and spikes with $\zeta \u22650,v\u22640$. To study how the bubbles and spikes are being formed, we consider the very early time dynamics, $t\u2009\u2248\u2009t0$.^{22,23} The solution for Eq. (9) linearized for $t\u2009\u2248\u2009t0$ is

This suggests the following changes of the morphology and the velocity of the interface near the tip.

For $v(t0)/v0>0$ and $\zeta 0k<0$, the interface becomes more curved and its velocity increases. For $v(t0)/v0>0$ and $\zeta 0k>0$, the interface flattens and the velocity decreases. For $v(t0)/v0<0$ and $\zeta 0k>0$, the interface becomes more curved and the velocity magnitude increases. For $\zeta 0k<0$ and $v(t0)/v0<0$, the interface flattens and the velocity magnitude decreases. For the acceleration-driven Rayleigh–Taylor dynamics with $a>\u22122$, the bubbles are formed for $\zeta 0k<0$ and the spikes are formed for $\zeta 0k>0$. The positions of RT bubbles and RT spikes are defined by the initial morphology of the interface, in excellent agreement with observations.^{11–23}

### B. Late-time nonlinear dynamics

#### 1. Asymptotic behavior of nonlinear solutions

For the late-time dynamics, we seek asymptotic solutions in the form of power-law functions,

where the exponents $\alpha ,\beta $ are defined from the dominant balance conditions in the system [Eq. (7)]. By substituting these expressions into the system [Eq. (7)], we find that for $a>\u22122$, the dominant balance conditions require the power-law exponents to be

This leads^{22} to asymptotic nonlinear solutions with

Note that for Richtmyer–Meshkov dynamics with $a<\u22122$, the dominant balance conditions require the power-law exponents to be $\beta =0,\alpha =\u22121$ and lead to asymptotic nonlinear solution with $\zeta \u223ck,\u2009v\u223c(k\tau )\u22121(t/\tau )\u22121,\u2009Mn,M\u0303n\u223ckn(k\tau )\u22121(t/\tau )\u22121$. For $a>\u22122$, the contribution of these terms is negligible, and Rayleigh–Taylor nonlinear dynamics is independent of the initial growth rate.

#### 2. Family of nonlinear asymptotic solutions

To find for $a>\u22122$ asymptotic nonlinear solutions with $\zeta \u223ck$ and $Mn,M\u0303n\u223ckn(t/\tau )a/2$, we retain higher order harmonics in the expressions for moments $Mn,M\u0303n$. Each of the moments is an infinite sum of the weighted Fourier amplitudes.^{11,12,22,36–38} To solve the closure problem, group theory is further applied. Particularly, the nonlinear asymptotic solutions are not unique and form a continuous family. Their multiplicity is due to the singular and non-local character of RT dynamics.^{35–38} The number of the family parameters is set by the flow symmetry. For group $p6mm$, the dynamics is highly isotropic, $z*\u223c\zeta (x2+y2)$, the interface morphology is captured by the principal curvature $\zeta $, and the family has one parameter.^{11,12,37,38} The multiplicity of the nonlinear asymptotic solutions is also due to the presence of shear $\Gamma $ at the interface [Eqs. (3) and (7)].^{22}

At $N=1$, for the family of the nonlinear solutions, the velocity $v$, shear function $\Gamma $, and Fourier amplitudes $\Phi m(\Phi \u0303m)$ depend on the curvature $\zeta $ as (see Table II for the explicit expressions)

Bubbles . | Spikes . |
---|---|

$v\u2009k\u2009\tau >0,\u2003\zeta k<0$ | $v\u2009k\u2009\tau <0,\u2003\zeta k>0$ |

$t\tau \u226b1,\u2003v=1k\tau (t\tau )a/2S,\u2003\zeta =\xi \u2009k,\u2003\Gamma =1\tau (t\tau )a/2Q,\u2003\Phi m\u2009(\Phi \u0303m)=1k\tau (t\tau )a/2\varphi m\u2009(\varphi \u0303m)$ | |

$S=\xb1(9\u221264\u2009\xi 2)2A\xi 48\xi \u2212A(9+64\xi 2)$ | |

$Q=12\u2009S9\u221264\u2009\xi 2,\u2003\varphi 1=\u22124S1+4\u2009\xi 3+8\u2009\xi ,\u2003\varphi 2=S1+8\u2009\xi 3+8\u2009\xi ,\u2003\varphi \u03031=4S1\u22124\u2009\xi 3\u22128\u2009\xi ,\u2003\varphi \u03032=\u2212S1\u22128\u2009\xi 3\u22128\u2009\xi $ |

Bubbles . | Spikes . |
---|---|

$v\u2009k\u2009\tau >0,\u2003\zeta k<0$ | $v\u2009k\u2009\tau <0,\u2003\zeta k>0$ |

$t\tau \u226b1,\u2003v=1k\tau (t\tau )a/2S,\u2003\zeta =\xi \u2009k,\u2003\Gamma =1\tau (t\tau )a/2Q,\u2003\Phi m\u2009(\Phi \u0303m)=1k\tau (t\tau )a/2\varphi m\u2009(\varphi \u0303m)$ | |

$S=\xb1(9\u221264\u2009\xi 2)2A\xi 48\xi \u2212A(9+64\xi 2)$ | |

$Q=12\u2009S9\u221264\u2009\xi 2,\u2003\varphi 1=\u22124S1+4\u2009\xi 3+8\u2009\xi ,\u2003\varphi 2=S1+8\u2009\xi 3+8\u2009\xi ,\u2003\varphi \u03031=4S1\u22124\u2009\xi 3\u22128\u2009\xi ,\u2003\varphi \u03032=\u2212S1\u22128\u2009\xi 3\u22128\u2009\xi $ |

Our dynamical system [Eq. (7)] can account for any number of harmonics $m$ in any order $N$. In higher orders, $N>1$, the solutions can be found similar to the $N=1$ case.^{11,12,37,38} Briefly, for $N>1$, in a broad range of curvature values $(\zeta /k)$, the asymptotic solutions exist and converge with the increase in $N$, the lowest-order amplitudes $|\Phi 1|,|\Phi \u03031|$ are dominant, and the values of amplitudes $|\Phi m|,|\Phi \u0303m|$ decay with the increase in $m$. Group theory analysis of RT dynamics with constant acceleration, $a=0$, finds that the $N=1$ nonlinear solutions properly capture the physical behavior.^{11,12,37,38} We address to the future the detailed investigation of mathematical properties of the nonlinear asymptotic solutions in RTI with variable acceleration at $N>1$.

The nonlinear asymptotic solutions in the family in Eq. (15) describe both the bubbles, which move up and are concave down, and the spikes, which move down and are concave up. Figure 3 illustrates the dependence of the velocity and shear function on the interface morphology (i.e., the curvature) for the nonlinear bubbles and the nonlinear spikes at various Atwood numbers in the family in Eq. (15). The curvature values are scaled with $\zeta cr=\u2212(3/8)k$ for the bubbles and with $\zeta \u0302cr=(3/8)k$ for the spikes. While the velocity and shear function of the asymptotic nonlinear solutions in Eq. (15) are time-dependent, when scaled with the values of the acceleration $g=Gta$ and the wavevector $k$, their magnitudes $|v/g/k|$ and $\Gamma /gk$ depend only on the interface morphology $|\zeta /k|$ and the Atwood number $\zeta /k$. Nonlinear solutions in Fig. 3 belong to the same family in Eq. (15), with $(\zeta /k)<0$ for bubbles and with $(\zeta /k)>0$ for spikes. Nonlinear bubbles are regular, with finite values $|v/g/k|$ and $|\Gamma /gk|$, whereas nonlinear spikes may be singular, with both the re-scaled values of the velocity $|v/g/k|$ and the interfacial shear $|\Gamma /gk|$ approaching infinity at some finite $(\zeta /k)$ (Fig. 3).

In Secs. III C–III E, we study in detail the properties of bubbles and spikes in the scale-dependent nonlinear dynamics of RTI with variable acceleration.

### C. Dynamics of nonlinear bubbles

#### 1. Nonlinear solutions

For nonlinear bubbles, we consider solutions with $v(\tau k)>0,(\zeta /k)<0$ in the family in Eq. (15) (Table II). In the family of asymptotic solutions with negative curvatures $\zeta /k<0$, the velocities are positive, $v(\tau k)>0$, and the shear functions are positive, $\Gamma \tau >0$. Nonlinear bubbles move down and are concave up, and vortical structures, which may appear in a vicinity of the bubble tip due to the interfacial shear, “rotate” from the heavy to the light fluid, in agreement with the experiment.^{14–18}

The solutions for RT bubbles with $(\zeta /k)<0,\u2009v(\tau k)>0,\Gamma \tau >0$ in Eq. (15) have the following properties. The domain of the function $v$ is $\zeta \u2208(\zeta cr,0),\u2009\zeta cr=\u2212(3/8)k$, where the subscript stands for critical. The range of the velocity function is $v\u2208(0,vmax)$, where the zero value $v=0$ is achieved at $\zeta =0$ and at $\zeta =\zeta cr$ and where the maximum value $v=vmax$ is achieved at $\zeta =\zeta max$ with $\zeta max\u2208(\zeta cr,0)$. The fastest solution in the family $(\zeta max,vmax)$ obeys the condition $\u2202v/\u2202\zeta |\zeta =\zeta max=0$ with $\u22022v/\u2202\zeta 2|\zeta =\zeta max<0$. The explicit dependence of the curvature $\zeta max$ and velocity $vmax$ on the Atwood number is complex. For $\zeta \u2208(\zeta cr,0)$, the shear function $\Gamma $ is the 1–1 function on the curvature $\zeta $ with $\Gamma \u2208(\Gamma min,\Gamma max)$. The minimum value $\Gamma min=0$ is achieved at $\zeta =0$, and the maximum value $\Gamma max$ is achieved at $\zeta =\zeta cr$. For $\zeta \u2208(\zeta cr,0)$ and $\Gamma \u2208(0,\Gamma max)$, the bubble velocity $v$ is a 1–1 function on the shear $\Gamma $. Explicit functions $v(\Gamma )$ and $\Gamma (v)$ are cumbersome.

For nonlinear RT bubbles, Figs. 4(a)–4(c) illustrate the dependence of the re-scaled velocity $v/vmax$ and the re-scaled shear function $\Gamma /\Gamma max$ on the re-scaled curvature $\zeta /\zeta cr$ and the dependence of the shear $\Gamma /\Gamma max$ on the velocity $v/vmax$ at some Atwood numbers $A$. The universality of nonlinear dynamics of RT bubbles and the link between the interface dynamics, interface morphology, and the interfacial shear are clearly seen from these plots. Figure 4(d) further illustrates the convergence of nonlinear RT bubbles with $\zeta =\zeta max$ by showing that the lowest-order amplitudes $|\Phi 1|,|\Phi \u03031|$ are dominant for a broad range of Atwood numbers $A$.

According to our results, the late-time dynamics of nonlinear RT bubbles is regular, in agreement with experiments (Figs. 3 and 4).^{14–18,22}

#### 2. Special solutions

The family of asymptotic solutions for nonlinear RT bubbles with $\zeta /k<0$ has some special solutions. We call these solutions the critical bubble, the convergence limit bubble, the Taylor bubble, the Layzer-drag bubble, the Atwood bubble, and the flat bubble. The properties of these solutions are illustrated in Fig. 5.

##### a. The critical bubble

For the critical bubble with curvature $\zeta =\zeta cr$, the solution is

This solution is special since the bubbles in the family of solutions cannot be more curved than the critical bubble. For the critical bubble, the velocity is zero and the shear function achieves its maximum value. The large shear can alone maintain the pressure at the interface (Fig. 5).

##### b. The convergence limit bubble

At $N=1$, the amplitudes of the heavy fluid relate as $|\Phi 1|>|\Phi 2|$ for $\zeta \u2208(\zeta cn,0)$ and as $|\Phi 1|=|\Phi 2|$ at $\zeta =\zeta cn$, where $\zeta cn=\u2212(5/24)k$ and sub-script stands for convergence. The amplitudes of the light fluid relate as $|\Phi \u03031|>|\Phi \u03032|$ for $\zeta \u2208(\zeta cr,0)$ and as $|\Phi \u03031|=|\Phi \u03032|$ at $\zeta =\zeta cr$. This defines the convergence limit bubble solution.

For the convergence limit bubble with curvature $\zeta =\zeta cn$, the solution is

At $N=1$, the convergence holds for $\zeta \u2208(\zeta cn,0)$ since the amplitudes of the heavy fluid are $|\Phi 1|>|\Phi 2|$ for $\zeta \u2208(\zeta cn,0)$ with $|\Phi 1|=|\Phi 2|$ at $\zeta =\zeta cn$ and the amplitudes of the light fluid are $|\Phi \u03031|>|\Phi \u03032|$ for $\zeta \u2208(\zeta cr,0)$ with $|\Phi \u03031|=|\Phi \u03032|$ at $\zeta =\zeta cr$. The ratio is $\zeta cn/\zeta cr=5/9$. For any Atwood number, the convergence limit bubble is less curved and has larger velocity and smaller shear when compared to the critical bubble (Fig. 5).

##### c. The Taylor bubble

The family has a solution with $\zeta =\u2212(1/8)k$. We call this solution the Taylor bubble since its curvature is the same as in the solution, which can be found in the study by Davies and Taylor (1950), upon modifying the wavevector $k=4\pi /3\lambda $ to $k=2\beta 1/\lambda $, where $\beta 1$ is the first zero of the Bessel function $J1$.^{4,11,12,37} For the Taylor bubble, the curvature, velocity, and shear function are

Note that for the Taylor bubble, the second Fourier amplitude of the heavy fluid is zero at $N=1$$(\Phi 2)T=0$. This property does not hold for $N>1$. The ratios are $\zeta T/\zeta cr=1/3$ and $\zeta T/\zeta cn=3/5$. For any Atwood number, the Taylor bubble is less curved and has larger velocity and smaller shear when compared to the convergent limit bubble and the critical bubble (Fig. 5).

##### d. The Layzer-drag bubble

The bubble family has a solution with velocity $v=(k\tau )\u22121(t/\tau )a/22A/(1+A)$. We call the solution the Layzer-drag bubble since the drag model applies this velocity re-scaling with the Atwood number to the single-mode Layzer's first-order approximation at $A=1$.^{33,39} Experiments and simulations tend to compare well with this re-scaling.^{33,39,47} For the Layzer-drag bubble, the curvature, velocity, and shear function are

The dependence of the Layzer-drag bubble curvature on the Atwood number is cumbersome (Fig. 5). For fluids with very different densities $A\u21921\u2212$, the solution for the Layzer-drag bubble is

For fluid with very similar densities $A\u21920+$, the solution for the Layzer-drag bubble is

The Layzer-drag bubble is more curved for fluids with similar densities $A\u21920+$ when compared to fluids with very different densities $A\u21921\u2212$.

For fluids with very different densities, $A\u21921\u2212$, the values of the curvature, velocity, and shear of the Layzer-drag bubble are the same as the values of the corresponding quantities of the Taylor bubble. For fluids with a finite density ratio, $A\u2208(0,1)$ the curvatures relate as $\zeta L/\zeta cr\u2208(9\u221243/3,1/3)$ and $\zeta L/\zeta T\u2208(9\u221243,1)$. For $0<A<1$, the Layzer-drag bubble is more curved and has smaller velocity and larger shear when compared to the Taylor bubble (Fig. 5).

##### e. The Atwood bubble

We call the fastest solution in the bubble family as the “Atwood bubble” to emphasize its complex dependence on the Atwood number, $\zeta A=\zeta max,\u2009vA=vmax$. The Atwood bubble solution has a remarkable invariant property,^{11,12,22,37,38}

For fluids with very different densities $A\u21921\u2212$, the curvature, velocity, and shear function of the Atwood bubble are

For fluid with very similar densities $A\u21920+$, the curvature, velocity, and shear function of the Atwood bubble are

For fluids with very different densities, $A\u21921\u2212$, the values of the curvature, velocity, and shear of the Atwood bubble are the same as the values of the corresponding quantities of the Taylor and the Layzer-drag bubbles. For fluids with a finite density ratio, the curvatures relate as $\zeta A/\zeta cr\u2208(0,1/3)$, $\zeta A/\zeta L\u2208(0,1)$, and $\zeta A/\zeta T\u2208(0,1)$. For $0<A<1$, the Atwood bubble is less curved and has larger velocity and smaller shear when compared to the Layzer-drag and the Taylor bubble (Fig. 5).

##### f. The flat bubble

For the flat bubble with curvature $\zeta =\zeta f$, the solution is

The curvature, velocity, and shear function of this bubble equal zero. The ratio is $\zeta f/\zeta cr=0$. Note that while the velocities of both the flat and critical bubbles are zero, the shear of the flat bubble is zero, whereas the shear of the critical bubble has a maximum value.

#### 3. Stability of solutions

To analyze the asymptotic stability of the solutions in the family, we slightly perturb $\zeta \u2192\zeta +\Delta \zeta $, $M\u2192M+\Delta M$, and $M\u0303\u2192M\u0303+\Delta M\u0303$ with $|\Delta \zeta /\zeta |\u226a1$, $|\Delta M/M|\u226a1$, and $|\Delta M\u0303/M\u0303|\u226a1$, such that $\Delta \zeta /\zeta \u223c\psi (t/\tau )$, $\Delta M/M\u223c\psi (t/\tau )$, and $\Delta M\u0303/M\u0303\u223c\psi (t/\tau )$, where $\psi $ is some function of time $(t/\tau )$, $|\psi |\u226a1$, to be defined. Upon substituting these expressions to the governing equations (7) and by linearizing the equations for small perturbations, we find that at $N=1$,

where $\chi $ is the function of the Atwood number and the bubble curvature. The real part of the function is negative $Re[\chi ]<0$ for stable solutions and is positive $Re[\chi ]>0$ for unstable solutions.^{12,37,38}

For $a>\u22122$, the function $\chi $ does not depend on the acceleration parameters. The explicit expression for this function $\chi $ is cumbersome. The detailed investigation of the dependence of $\chi $ on the Atwood number, and the bubble curvature suggests that for a finite Atwood number, the flattened bubbles are unstable with $Re[\chi ]>0$, whereas the curved bubbles are stable with $Re[\chi ]<0$, including the Atwood bubble, the Layzer-drag bubble, the Taylor bubble, and the convergence limit bubble.^{22} The $N>1$ analysis is similar to Refs. 37 and 38, and we address it to the future.

### D. Dynamics of nonlinear spikes

#### 1. Nonlinear solutions

For nonlinear spikes, we consider solutions with $v(\tau k)<0,(\zeta /k)>0$ in the family in Eq. (15) (Table II). In the family of asymptotic solutions with positive curvatures $\zeta /k>0$ in Eq. (15), velocity values are negative, $v(\tau k)<0$, and the shear function values are negative, $\Gamma \tau <0$. Nonlinear spikes move down and are concave up; vortical structures, which may appear in a vicinity of a spike due to the interfacial shear, rotate from the light to the heavy fluid, in agreement with experiments.^{14–18}

The solutions with $(\zeta /k)>0$, $v(\tau k)<0$, and $\Gamma \tau <0$ have the following properties. The domain of the velocity function $v$ is $\zeta \u2208(\zeta \u0302max,\zeta \u0302cr)$ with $\zeta \u0302max\u2208\u230a0,\zeta \u0302cr\u230b$. For $\zeta \u2208(\zeta \u0302max,\zeta \u0302cr)$, the range of the function $v$ is $v\u2208(v\u0302max,0)$ with $v\u21920$ for $\zeta \u2192\zeta \u0302cr$ and with $v\u2192v\u0302max$ for $\zeta \u2192\zeta \u0302max$ (Fig. 3). The critical curvature value is $\zeta \u0302cr=(3/8)k$. The curvature value $\zeta \u0302max$ bounds the interval of values $\zeta $, in which the function $v$ is real and non-positive. The value $\zeta \u0302max$ is the pole of the velocity function, as well as the amplitudes and shear function, and it is $\zeta \u0302max/k=3(1\u22121\u2212A2)/(8A)$. It corresponds to the spike with the largest value of the velocity magnitude $v\u0302max$, with $v\u0302max\u2192\u2212\u221e$. The shear function is $\Gamma \u2208(\u2212|\Gamma \u0302max|,\u2212|\Gamma \u0302min|)$ for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cr)$, achieving the value with the minimum magnitude $\Gamma \u2192\Gamma \u0302min$ and $\Gamma \u0302min<0$, for $\zeta \u2192\zeta \u0302cr$, and achieving the value with the maximum magnitude of the shear $\Gamma \u2192\Gamma \u0302max$ and $\Gamma \u0302max<0$, approaching negative infinity, $\Gamma \u0302max\u2192\u2212\u221e$, for $\zeta \u2192\zeta \u0302max$ (Fig. 3).

For nonlinear RT spikes, Figs. 6(a)–6(c) illustrate the dependence of the re-scaled velocity $v/v\xafmax$ and the re-scaled shear function $\Gamma /\Gamma \xafmax$ on the re-scaled curvature $\zeta /\zeta \u0302cr$ and the dependence of the shear $\Gamma /\Gamma \xafmax$ on the velocity $v/v\xafmax$ for some Atwood numbers $A$, respectively. Here, the scaling values are regularized as $v\xafmax=v\u0302max((\zeta \u2212\zeta \u0302max)/k)1/2$ and $\Gamma \xafmax=\Gamma \u0302max((\zeta \u2212\zeta \u0302max)/k)1/2$. The universality of nonlinear dynamics of RT spikes and the link between the interface dynamics, interface morphology, and the interfacial shear are clearly seen from these plots. Figure 6(d) further illustrates the convergence of solutions for nonlinear of RT spikes with $\zeta =\zeta \u0302max$ by showing that for a broad range of values of the Atwood numbers $A$, the lowest-order amplitudes $|\Phi 1|,|\Phi \u03031|$ are dominant even at the (singular) point $\zeta =\zeta \u0302max$.

According to our results, late-time dynamics of nonlinear RT spikes may exhibit singular behavior, in agreement with experiments (Figs. 3 and 6).^{11–15,17,18}

#### 2. Special solutions

The family of asymptotic solutions for nonlinear RT spikes with $\zeta /k>0$ and $v(\tau k)<0$ has some special solutions. We call these solutions the critical spike, the convergence limit spike, the Taylor spike, the Layzer-drag spike, the Atwood spike, and the flat spike. The properties of these solutions are illustrated in Fig. 7.

##### a. The critical spike

For the critical spike with curvature $\zeta =\zeta \u0302cr$, the solution is

For the critical spike, the velocity is zero for any density ratio, and the shear function has a non-positive value, which is zero for $A\u21920$ and is finite for $0<A<1$ and which is singular approaching negative infinity for $A\u21921$ (Fig. 7).

The magnitudes of the curvature and velocity of the critical spike are the same as those for the critical bubble, with $\zeta \u0302cr=\u2212\zeta cr$ and $v\u0302cr=vcr=0$. The magnitude of the shear of the critical spike differs from that of the critical bubble: The shear's magnitude of the critical spike approaches infinity as $A\u21921\u2212$, whereas the shear's magnitude of the critical bubble is finite for any $A$ (Figs. 5 and 7).

##### b. The convergence limit spike

At $N=1$, the amplitudes of the heavy fluid relate as $|\Phi 1|>|\Phi 2|$ for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cr)$, with $|\Phi 1|=|\Phi 2|=0$ at $\zeta \u0302cr=(3/8)k$. The light fluid amplitudes relate as $|\Phi \u03031|>|\Phi \u03032|$ for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cn)$, with $|\Phi \u03031|=|\Phi \u03032|$ at $\zeta \u0302cn=(5/24)k$ (sub-script stands for convergence). This defines the solution for the convergence limit spike.

For the convergence limit spike with curvature $\zeta =\zeta \u0302cn$, the solution is

At $N=1$, the convergence holds for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cn)$ since the amplitudes are $|\Phi \u03031|>|\Phi \u03032|$ and $|\Phi 1|>|\Phi 2|$ for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cn)$ with $|\Phi \u03031|=|\Phi \u03032|$ at $\zeta =\zeta \u0302cn$, whereas the amplitudes $|\Phi 1|>|\Phi 2|$ for $\zeta \u2208(\zeta \u0302max,\zeta \u0302cr)$ (Figs. 6 and 7).

For the convergence limit spike, the velocity and shear function have non-positive values, which are zero for $A\u21920+$ and finite for $0<A<45/53$ and which become singular approaching negative infinity for $A\u219245/53\u22480.85$. The ratio is $\zeta \u0302cn/\zeta \u0302cr=5/9$, and the convergence limit spike is less curved when compared to the critical spike (Fig. 7).

The magnitude of the curvature of the convergence limit spike is the same as that for the convergence limit bubble, with $\zeta \u0302cn=\u2212\zeta cn$. The magnitudes of the velocity and shear of the convergence limit spike differ from their respective values of the convergence limit bubble: The magnitudes of the velocity and shear of the spike approach infinity for $A\u219245/53$, whereas the magnitudes of the velocity and shear of the bubble are finite for any $A$. Note that the convergence limit spike is defined by the condition $|\Phi \u03031|=|\Phi \u03032|$ for the amplitudes of the light fluid, whereas the convergence limit bubble is defined by the condition $|\Phi 1|=|\Phi 2|$ for the amplitudes of the heavy fluid (Figs. 3, 5, and 7).

##### c. The Taylor spike

The family has a solution with $\zeta =(1/8)k$. We call this solution the Taylor spike, by analogy with the Taylor bubble.^{11,12,37}

For the Taylor spike with curvature $\zeta =\zeta \u0302T$, the solution is

For this special solution, the second Fourier amplitude of the light fluid is zero $(\Phi \u03032)T=0$ at $N=1$. This property does not hold for $N>1$ (Fig. 7).

For the Taylor spike, the velocity and shear function have non-positive values, which are zero for $A\u21920$ and finite for $0<A<3/5$ and which approach negative infinity and become singular for $A\u21923/5=0.6$. Since the ratios are $\zeta \u0302T/\zeta \u0302cr=1/3$ and $\zeta \u0302T/\zeta \u0302cn=3/5$, the Taylor spike is less curved when compared to the critical spike and the convergence limit spike (Fig. 7).

The magnitude of the curvature of the Taylor spike is the same as that for the Taylor bubble, with $\zeta \u0302T=\u2212\zeta T$. The magnitudes of the velocity and shear of the Taylor spike differ from their respective values of the Taylor bubble: The velocity and shear magnitudes for the spike approach infinity for $A\u21923/5$; the velocity and shear magnitudes for the bubble are finite for any $A$ (Figs. 5 and 7).

##### d. The Layzer-drag spike

The RT family of spikes has a solution with velocity $v=\u2212(k\tau )\u22121(t/\tau )a/22A/(1\u2212A)$. We call this solution the Layzer-drag spike since the drag model applies such re-scaling with the Atwood number for the velocity of the nonlinear spike to the Layzer first order approximation of the velocity of the nonlinear bubble at $A=1$.^{33,39} Experiments and simulations tend to compare well with this re-scaling.^{39,47} For the Layzer-drag spike, the curvature, velocity, and shear function are

The dependence of the Layzer-drag spike curvature on the Atwood number is cumbersome. The velocity and shear function of the Layzer-drag spike have non-positive values, which are zero for $A\u21920$ and finite for $0<A<1$ and which become singular approaching negative infinity for $A\u21921$.

For fluids with very different densities $A\u21921\u2212$, the solution for the Layzer-drag spike is

The velocity and shear of the Layzer-drag spike are singular for $A\u21921\u2212$, with the shear being even more singular than the velocity.

For fluid with very similar densities $A\u21920+$, the solution for the Layzer-drag spike is

The Layzer-drag spike is more curved for fluids with very different densities $A\u21921\u2212$ than for fluids with very similar densities $A\u21920+$. For a two fluid system with $A\u2208(0,1)$, the ratios are $\zeta \u0302L/\zeta \u0302cr\u2208(9\u221243/3,1)$ and $\zeta \u0302L/\zeta \u0302T\u2208(9\u221243,3)$, and the Layzer-drag spike is less curved than the critical spike and is more curved than the Taylor spike (Fig. 7). The Layzer-drag spike has the smaller magnitude of the velocity and smaller magnitude of the shear when compared to those of the Taylor spike. The curvatures of the Layzer-drag spike and convergence limit spike equal one another $\zeta \u0302L/\zeta \u0302cn=1$ at $A=235/451\u2009\u2248\u20090.52$, and for larger values of the Atwood numbers, the Layzer-drag spike solutions are not convergent (Fig. 7).

For fluids with very similar densities, $A\u21920$, the magnitude of the curvature of the Layzer-drag spike is the same as that for the Layzer-drag bubble, with $\zeta \u0302L=\u2212\zeta L$. The magnitudes of the velocity and shear of the Layzer-drag spike differ from those of the Layzer-drag bubble: The velocity and shear magnitudes of the spike approach infinity for $A\u21921$, whereas the velocity and shear magnitudes are finite for any $A$ (Figs. 5 and 7).

##### e. The Atwood spike

The fastest spike in the family of solutions is $(\zeta \u0302max,v\u0302max)$. We call this solution the “Atwood spike” due to its complex dependence on the Atwood number, $\zeta \u0302A=\zeta \u0302max,\u2009v\u0302A=v\u0302max$. This solution is

For the Atwood spike, the curvature is regular, $\zeta \u2192\zeta \u0302A$, whereas the velocity and shear function are singular, with $v\u0302A\u2192\u2212\u221e,\Gamma \u0302A\u2192\u2212\u221e$ for any Atwood number. The functions $\kappa ,\phi ,\gamma $ for the curvature, velocity, and shear are, respectively,

For fluids with very different densities $A\u21921\u2212$, these functions behave as

Note the singular character of the function $\gamma $ and hence $\Gamma \u0302A$ for $A\u21921\u2212$. For fluids with very similar densities $A\u21920+$, these functions behave as

The Atwood spike is less curved than the Layzer-drag spike for $A\u2208[0,1)$ and has the same curvature as the Layzer-drag spike $\zeta \u0302A/\zeta \u0302L=1$ at $A=1$. The Atwood spike attains the same curvature as the critical spike $\zeta \u0302A/\zeta \u0302cr=1$ at $A=1$, the same curvature as the Taylor spike $\zeta \u0302A/\zeta \u0302T=1$ at $A=3/5\u2009\u2248\u20090.6$, and the same curvature as the convergence limit spike $\zeta \u0302A/\zeta \u0302cn=1$ at $A=45/53\u2009\u2248\u20090.85$ (Fig. 7). For $A>45/53$, the Atwood spike solutions may not converge.

The properties of the curvature, velocity, and shear Atwood spike differ dramatically from those of the Atwood bubble: For any density ratio $A$, the magnitudes of the velocity and shear of the Atwood spike are singular, whereas the magnitudes of the velocity and shear of the Atwood bubble are regular and finite. While the magnitudes of the curvature of both the Atwood spike and the Atwood bubble are finite, the Atwood spike is more curved (less curved) when compared to the Atwood bubble for fluids with very different (very similar) densities $A\u21921\u2212$ ($A\u21920+$) (Figs. 5 and 7).

##### f. The flat spike

For the flat spikes with curvature $\zeta =\zeta \u0302f$, the solution is

The curvature, velocity, and shear function of this spike equal zero. The ratio is $\zeta \u0302f/\zeta \u0302cr=0$. Note that while the value $\zeta \u0302A\u2192\zeta \u0302f$ for $A\u21920+$, the flat spike remains an isolated solution, in contrast to the flat bubble.

#### 3. Stability of solutions

The stability of asymptotic solutions in the family of spikes can be studied similarly to the case of the family of bubbles [Eq. (22)]. The stability properties of the spikes differ significantly from those of the bubbles. Particularly, in the family of spikes, for a given Atwood number, there is a narrow interval of stable solutions with finite curvatures, whereas fast solutions with $\zeta \u2192\zeta \u0302A$ are unstable. This suggests that nonlinear spikes can be in transient and move faster, with $|v\tau k\u2009|\u223c(t/\tau )\delta ,\u2009\delta >a/2$, than it is prescribed by the nonlinear asymptotic dynamics, with $|v\tau k\u2009|\u223c(t/\tau )\u2009a/2$.

### E. Properties of Rayleigh–Taylor dynamics

#### 1. Interplay of acceleration, morphology, and shear

By analyzing the properties of the linear and nonlinear dynamics [Eqs. (7)–(29)] for $a>\u22122$, we find that in the linear regime, RT dynamics is acceleration driven: the growth rate of the interface perturbations is defined by the acceleration, whereas the initial morphology of the interface defines the positions of bubbles and spikes [Eqs. (9) and (10) and Table I].^{22,23} The nonlinear dynamics is more complex. From the properties of the solutions' family [Eq. (15)], including the special solutions for RT bubbles [Eqs. (16)–(21)] and for RT spikes [Eqs. (23)–(29)], we find that the asymptotic nonlinear RT dynamics is influenced by the interplay of the acceleration, the interface morphology, and the interfacial shear (Table II and Figs. 3–7).

##### a. RT bubbles

From the properties of the solutions describing RT bubbles, it follows that the curved bubbles move faster than the flattened bubbles due to the larger effective acceleration [Eqs. (15) and (16)–(21) and Fig. 4]. The larger curvature leads to larger interfacial shear. The larger shear, in turn, stronger decelerates the bubble. For a well curved bubble, the shear is large, the deceleration is strong, and the bubble velocity decreases (Fig. 4). The shear is the 1–1 function on the curvature, and for the critical bubble, the very large shear can alone maintain the pressure at the interface (Figs. 4 and 5).

For any Atwood number, the bubble velocity as the function of the interfacial shear is very nearly linear for small shear values, then it achieves its maximum value corresponding to the Atwood bubble, and then, it sharply drops to zero for large shear values [Eq. (20) and Fig. 4]. In a vicinity of the fastest moving Atwood bubble solution, the velocity has a very steep dependence on the interfacial shear. This suggests the need in highly accurate methods of numerical modeling and experimental diagnostics for accurate capturing of RT dynamics.^{13–21}

##### b. RT spikes

Similar to RT bubbles, RT spikes move slower for larger magnitudes of their curvatures [Eqs. (15) and (23)–(29) and Fig. 6]. For fluids with a finite density ratio, the shear of the strongly curved critical spike can alone maintain the pressure at the interface, reducing the spike velocity to zero value (Fig. 6). In contrast to RT bubbles, RT spikes are singular (Figs. 3 and 6). The velocity and shear values of the convergence limit and the Taylor spikes become singular approaching negative infinity at some Atwood numbers (Fig. 7). Those of the Layzer-drag spike approach negative infinity for the Atwood number $A\u21921\u2212$, with the shear being more singular than the velocity (Fig. 7).

The fastest moving Atwood spike is singular for any Atwood number [Eqs. (27) and (28) and Fig. 7]. For the Atwood spike, both the velocity and shear approach infinite values in the case of fluids with a finite density ratio $0<A<1$, and the shear is even more singular than the velocity in the case of fluids with highly contrasting densities $A\u21921\u2212$. The singular character and instability of the Atwood spike solution suggest that RT spike, upon achieving a finite curvature value, can become unstable and move faster than it is prescribed by the nonlinear asymptotic dynamics [Eq. (14)].

#### 2. Single-scale and multi-scale dynamics

By analyzing the properties of the early-time linear dynamics [Eq. (10)], we find that for ideal incompressible fluids and for the interface perturbations with very small initial amplitude RT, dynamics is singe-scale and is defined by the horizontal length scale, which is the spatial period (the wavelength) of the interface perturbation (Fig. 2).

The nonlinear dynamics is more complex. From the properties of the solutions' family in Eq. (15), we find that, the asymptotic nonlinear RT dynamics is essentially multi-scale and is defined by the two macroscopic scales—the horizontal length scale (the wavelength) and the vertical length scale (the amplitude). This can be understood by viewing the RT coherent structure as a standing wave, whose nonlinear dynamics is influenced by its wavelength and by its growing amplitude.

##### a. RT bubbles

In the family of the asymptotic solutions describing nonlinear RT bubbles, we define the fastest stable solution as the physically significant solution, $\zeta \u223c\zeta A,v\u223cvA,\Gamma \u223c\Gamma A$. The Atwood bubble has the invariant property, $vA\tau \u2009k\u2009(t/\tau )\u2212a/2(8|\zeta A|/k)\u22123/2=1$ [Eq. (20)]. The left-hand side of the expression is a function of time and of the vertical and the horizontal length scales and their derivatives, whereas the right-hand side is unity. This invariance implies that the two macroscopic scales—the wavelength $\lambda $ and amplitude $z0$—contribute to the multi-scale dynamics of nonlinear RT bubbles.^{11,12,22,37,38}

##### b. RT spikes

The fastest Atwood spike has no invariant property similar to the fastest Atwood bubble. The singular character of the dynamics of the (unstable) fastest Atwood spike with $|v\u0302A|(Gta/k)\u22121/2\u2192\u221e$ suggests that the nonlinear RT spike can move faster than the nonlinear asymptotic dynamics prescribes [Eqs. (27) and (28)]. As such, it is influenced even more by the vertical scale (the amplitude) when compared to the nonlinear RT bubbles. The dynamics of the nonlinear RT spike is thus multi-scale, with the two macroscopic length scale contributing—the vertical length scale (amplitude) and the horizontal length scale (wavelength). To our knowledge, this is the first report of the multi-scale character of nonlinear dynamics of RT spikes based on the properties of nonlinear asymptotic solutions for the boundary value problems [Eqs. (1)–(4)].

#### 3. Transition to the self-similar mixing

The linear and nonlinear RT dynamics are scale-dependent. As time evolves, the flow transits to the scale-invariant regime, in which the vertical scale grows self-similarly with time.^{3,4,11–22,39–47} The group theory approach studies the properties of self-similar RT dynamics within the frame of the momentum model.^{11,14,22,42,55} This model has the same symmetries and scaling transformations as the governing equations and yields asymptotic solutions for the scale-dependent and self-similar dynamics, up to constant.

The important outcome of the momentum model is that RT dynamics can transit from the scale-dependent nonlinear regime to the scale-invariant self-similar mixing regime when the vertical scale (the amplitude) is the dominant scale, in excellent agreement with high Reynolds number experiments.^{11,14,42} This mechanism complements the traditional merger mechanism, in which the transition to the self-similar mixing is believed to occur via the growth of horizontal scales.^{39–41} Moreover, for RT coherent structures with hexagonal symmetry, the growth of horizontal scales may not be feasible since it may lead to anisotropy of the highly isotropic hexagonal pattern.^{56} Reference 56 provides the detailed classification of the formation of patterns in RT flows and possible discrete transitions with the increase in the wavelength(s) of RT coherent structures under the influence of modulations. These results excellently agree with the existing experiments.^{13–18}

We address to the future the identification of a direct link between the analytical results in the present paper and momentum model. Here, we analyze how the multi-scale character of the nonlinear dynamics is associated with and may lead to the acceleration of RT bubbles and RT spikes.

##### a. RT bubbles

Recall that in the laboratory frame of reference, the interface is described near the tip of the bubble as $Z=z0+\zeta (x2+y2)$ with the velocity $v=z\u03070$. For the fastest Atwood bubble, the curvature is $\zeta \u223c\zeta A$ and the velocity is $v\u223cvA$, and their values are $|\zeta A|\u223ck\u223c\lambda \u22121$ and $|vA|\u223c(k\tau )\u22121(t/\tau )a/2$, as found by the asymptotic balance $v2\u223cg\u2009|\zeta |\u223cg/\lambda $, leading to $|v|\u223cG\lambda ta$ and $|z0|\u223ctG\lambda ta$ [Eqs. (14) and (15)]. The invariant property of the Atwood bubble [Eq. (20)] suggests that the amplitude $z0$ and spatial period $\lambda $ contribute independently to the multi-scale dynamics of the nonlinear bubble. Hence, the amplitude $z0$ can be a dominant scale defining the rate of loss of specific momentum (i.e., the drag force per unit mass) as $\u223cv2/|z0|$.^{42,55} The effective reduction of the drag force from $\u223cv2/\lambda $ to $\u223cv2/|z0|$ can lead (for $|z0|>\lambda $) to the bubble acceleration and transit the bubble from the scale-dependent nonlinear regime with $|v|\u223cG\lambda ta$ and $|z0|\u223ctG\lambda ta$ to self-similar mixing with $|v|\u223cGta+1$ and $|z0|\u223cGta+2$.^{11,14,22,42}

##### b. RT spikes

Similar to RT bubbles, in the laboratory frame of reference, the interface is described near the tip of the spike as $Z=z0+\zeta (x2+y2)$ with the velocity $v=z\u03070$. For the Atwood spike, the curvature is $\zeta \u0302\u223c\zeta \u0302A$ and the velocity is $v\u223cv\u0302A$ with the finite value of the curvature $|\zeta \u0302A|\u223ck\u223c\lambda \u22121$ and with the singular value of the velocity approaching infinity $|v\u0302A|\u223c(k\tau )\u22121(t/\tau )a/2(k/(\zeta \u2212\zeta \u0302A))1/2\u2192\u221e$ [Eqs. (27) and (28)]. This suggests that for the Atwood spike, the amplitude $|z0|$ can increase faster than it is prescribed by the asymptotic balance $v2\u223cg\u2009|\zeta |/k2\u223cg\lambda $ with $|v|\u223cG\lambda ta$ and $|z0|\u223ctG\lambda ta$ in Eqs. (14) and (15). Thus, the amplitude $z0$ can be the dominant scale defining the rate of loss of specific momentum (i.e., the drag force per unit mass) as $\u223cv2/|z0|$. The reduction of the rate of loss of specific momentum (the drag force per unit mass) from $\u223cv2/\lambda $ to $\u223cv2/|z0|$ can lead (for $|z0|\u226b\lambda $) to the acceleration of the spike and transits the spike from the scale-dependent nonlinear dynamics to the self-similar mixing with $|z0|\u223cGta+2$ and $|v|\u223cGta+1$.^{11,14,22,42,55}

For fluids with very different densities, $A\u21921\u2212$, the Atwood spike is highly singular, suggesting that its dynamics can only be self-similar, in agreement with experiments.^{4,13–19} To our knowledge, this is the first analysis of the mechanism of transition to self-similar mixing on the basis of properties of rigorous solutions for nonlinear dynamics of RT bubbles and RT spikes in Eqs. (1)–(4).

#### 4. Interfacial dynamics

By accurately accounting for the interplay of harmonics and by systematically connecting the interface velocity to the interfacial shear in a broad range of the acceleration parameters, our analysis finds that RT dynamics is essentially interfacial: In both the linear and nonlinear regimes, for both RT bubbles and RT spikes, the flow has intense motion of the fluids in a vicinity of the interface and has effectively no motion away from the interface. The velocity field is potential in the fluid bulk. Due to the presence of shear at the interface, shear-driven vortical structures develop at the interface. Figure 8 illustrates the velocity fields in the vicinity of the tips of RT bubbles and RT spikes.

This velocity pattern is in excellent agreement with experiments and simulations.^{13–21} For realistic fluids, in order to accurately grasp the formation and the dynamics of vortical structures, one should also account for the influence of viscosity, compressibility, surface tension, mass flux across the interface, magnetic field, and other factors. We address these tasks to the future.

#### 5. Effect of interfacial shear

Our theory finds that the interfacial shear plays the important role in RT dynamics. The presence of interfacial shear is fully consistent with the conversation laws, including the continuity of the tangential component of momentum at the interface and, in the case of zero interfacial mass flux, the discontinuity of the tangential component of velocity at the interface [Eqs. (1)–(3)]. The presence of shear at the interface leads to the appearance of the shear-driven vortical structures at the interface (Fig. 8).^{13–18} It is a common wisdom that RT evolution is accompanied by shear-driven Kelvin–Helmoltz instability (KHI) developing due to the jump of the tangential component of velocity at the interface.^{13–21} Our results can be used to explain and predict the properties of KHI-driven vortical structures in RT flow.

First, according to our results, the jump of tangential velocity at the interface is zero at the tips of the bubble and spike, and it increases away from these tips. Hence, the interface is free from vortical structures exactly at the tips of the bubble and spike. KHI may develop away from the tips of the bubble or spike, leading to the appearance of shear-driven vortical structures on the sides of evolving bubbles and spikes. This theoretical explanation agrees excellently with experiments.^{13–18}

Second, the shear function, which we establish in our work, can be used to estimate the growth rate of KHI, as $\omega KHI\u223c|\Gamma |$ in the region of bubbles and $\omega \u0302KHI\u223c|\Gamma \u0302|$ in the region of spikes, and to find the dependence of $\omega KHI,\u2009\omega \u0302KHI$ on the density ratio, interface morphology, flow symmetry, and acceleration parameters. According to our results, the shear function is regular in the region of bubbles and is singular in the region of spikes. Hence, one expects a more intense production of vortical structures for RT spikes when compared to RT bubbles, in excellent agreement with the existing observations.^{13–21} To our knowledge, this is the first theoretical explanation of the bubble/spike asymmetry of vortical structures in RT flow.

#### 6. Formal properties of nonlinear solutions

By directly linking the interface velocity, the flow fields, the interfacial shear function, and the interface morphology in RT nonlinear dynamics, our theory finds the physics interpretation and resolves the long-standing problem of multiplicity of solutions for nonlinear RTI.^{33–37} According to our results, for RT dynamics with variable acceleration, there is a continuous family of asymptotic nonlinear solutions, with each solution having its own curvature, velocity, the Fourier amplitudes, the flow fields, and the interfacial shear.^{22} The multiplicity of nonlinear asymptotic solutions for the boundary value problem is due to the presence of shear at the interface and is the essential property of RT dynamics associated with the governing equations and the boundary conditions at the interface [Eqs. (1)–(4)].^{11,22,35} In the continuous family of nonlinear solutions, asymptotic solutions with negative curvature and positive velocity correspond to regular RT bubbles (which move up and concave down), whereas asymptotic solutions with positive curvature and negative velocity correspond to singular RT spikes (which move down and concave up) (Figs. 1 and 8). To our knowledge, our work is the first to report that asymptotic solutions describing dynamics of the nonlinear RT spikes and the nonlinear RT bubbles with variable acceleration form a continuous family.

#### 7. Effect of the acceleration parameters

According to our results, for a broad range of acceleration parameters, $a>\u22122,G>0$, the interface dynamics is described by tabular special functions, such as modified Bessel functions and power-law functions. The early-time linear dynamics is super-exponential (i.e., quicker than exponential) for $a>0$, is exponential at $a=0$, and is sub-exponential (i.e., slower than exponential) for $\u22122<a<0$ [Eq. (10), Table I, and Fig. 2]. For the late-time nonlinear dynamics, the velocity near the regular point of the interface increases with time for $a>0$, is steady at $a=0$, and decreases with time for $a<0$ [Eq. (15), Table II, and Fig. 3]. Furthermore, except for the difference in the time-dependence, nonlinear dynamics depends only on the interface morphology and flow symmetry.

These properties enable a comparative study of RT dynamics for various acceleration parameters, $a>\u22122,G>0$. Particularly, by analyzing the properties of fast RT dynamics for large exponents $a>0$, one can deduce the properties of slow RT dynamics for small exponents $\u22122<a<0$. Such analysis is especially important in studies of RTI in high energy density plasmas (for instance, in supernovae and fusion), where RTI is driven by an explosion or an implosion with acceleration exponents $\u22122<a<\u22121$ defined by blast waves.^{22,50–54} Our results thus provide new extensive theory benchmarks for the future.

#### 8. Comparisons with the existing observations

Significant efforts are undertaken in the RT community to better understand the scale-dependent RT dynamics and self-similar RT mixing. See edited research books^{1,2} and experimental and numerical papers for details.^{11–25,45–47,55–60} As discussed earlier in Secs. III A 1–III E 7, our theory is in excellent qualitative agreement with the existing observations^{11–25,45–47,55–60} and provides new extensive quantitative benchmarks for experiments and simulations. To date, to the best of the authors' knowledge, accurate experiments and simulations are free from systematic investigations of RT dynamics driven by variable acceleration for three-dimensional flow with hexagonal symmetry.^{13–21}

We consider RT dynamics with variable acceleration for three-dimensional flow with hexagonal symmetry since hexagonal symmetry is important to study from the physics perspectives and since it has a number of advantages when compared to other space groups. For instance, the hexagonal group is the most isotropic when compared to other space groups; for structures with hexagonal symmetry, transitions to super-structures with larger periods lead to a loss of isotropy and may thus be infeasible.^{10–12} At the same time, structures with hexagonal symmetry may be challenging to accurately implement in experiments and simulations. To the best of the authors' knowledge, the existing experiments and numerical simulations are absent for such a setup.^{13–21} Our theory results provide new benchmarks for the future research.

For instance, the agreements and departures of our scale-dependent solutions for the linear and nonlinear dynamics of RT bubbles and RT spikes from those in experiments and numerical simulations can serve to better understand RT dynamics in realistic fluids. Note that our theory focuses on the large-scale dynamics and presumes that interfacial vortical structures are small. This assumption is applicable for fluids with very different densities $A\u21921\u2212$ and with finite density ratios $0\u226aA<1$. For fluids with very similar densities $A\u21920+$, other approaches can be employed.^{11,12,22}

Our theory finds the effect of the interplay of the acceleration, the interface morphology, and the interfacial shear on the scale-dependent dynamics of RTI with variable acceleration. Our results suggest that this interplay can be more complex than it is traditionally seen by the empirical Layzer-type buoyancy-drag models.^{25,39,41,47} A direct link of the scale-dependent RT dynamics and the self-similar RT is thus the important task. The group theory approach can be applied to solve this task, to be done in the future.^{11,14,22}

Our theory focuses on RT dynamics in ideal incompressible fluids. In realistic circumstances, other factors should be accounted for, e.g., viscosity, compressibility, surface tension, and magnetic field. Qualitative and quantitative influences of these factors should be systematically studied.

For instance, viscosity may influence the instability growth rate in the linear regime and the realization of asymptotic solutions in the nonlinear regime. Particularly, while the Layzer-drag bubble has a smaller velocity when compared to the Atwood bubble, it also has a larger shear and may thus “survive” under the viscous effects (Fig. 5). Likewise, viscosity may “stabilize” for a finite time and realize the nonlinear spike with asymptotic velocity $|v|\u223cG\lambda ta$ before the spike transits to the self-similar mixing (Fig. 7).

#### 9. Theory benchmarks

Our analysis elaborates qualitative and quantitative benchmarks of RT dynamics not discussed before. These are the fields of velocity and pressure, the interface morphology and the curvatures of the bubble and spike, the interfacial shear and its link to the velocity and the curvature of the bubble and the spike, the spectral properties of the velocity and pressure, and the interface growth and growth rate. By diagnosing the dependence of these quantities on the density ratio, the flow symmetry, and the acceleration's exponent and strength, by identifying their universal properties, and by accurately measuring departures of data in real fluids from theoretical solutions in ideal fluids, one can further advance the knowledge of RT dynamics in realistic environments, achieve better understanding of RT relevant processes, from supernovae to fusion, and improve methods of numerical modeling and experimental diagnostics of interfacial dynamics in fluids, plasmas, and materials under conditions of low and in high energy densities.

## IV. DISCUSSION AND CONCLUSION

In this work, we focused on the scale-dependent dynamics of Rayleigh–Taylor instability driven by variable acceleration with a power-law time-dependence for the acceleration exponents greater than (−2) and studied the spatially periodic three-dimensional flow with hexagonal symmetry in the plane normal to the acceleration direction. We applied the group theory approach to solve the boundary value problem, by expanding the flow fields with the use of irreducible representations of the group, by deriving the dynamical system from the governing equations, and by finding solutions for the system to describe the dynamics of the RT coherent structure of bubbles and spikes in the linear and nonlinear regimes [Eqs. (1)–(29), Tables I and II, and Figs. 1–8].

For the early-time linear dynamics, we provide the dependence of the RTI growth rate on the acceleration parameters and show that the positions of bubbles and spikes in RT flow are defined by the initial morphology of the interface [Eqs. (10) and (11), Table I, and Fig. 2]. For the late-time nonlinear dynamics, we find the continuous family of asymptotic solutions, directly link the interface dynamics to the interfacial shear, and scrupulously study the properties of these solutions [Eqs. (15)–(29), Table II, and Figs. 3–8]. The essentially multi-scale and interfacial character of the coherent dynamics is demonstrated. The former can be understood by viewing the RT coherent structure as a standing wave with the growing amplitude. The latter implies that RT flow has effectively no motion of the fluids away from the interface and intense motion of the fluids near the interface, with shear-driven vortical structures appearing at the interface (Fig. 8). Our theory provides the unified framework for the scale-dependent RT dynamics [Eqs. (1)–(29)]. Our solutions describe in the linear and nonlinear regime the RT bubbles, which move up and concave down, and the RT spikes, which move down and concave up. We find that in the linear regime, RT bubbles and RT spikes grow nearly symmetrically (Fig. 2). In the nonlinear regime, the dynamics of RT bubbles is regular and the dynamics of RT spikes is singular (Fig. 3). We also reveal that RT bubbles and RT spikes may further accelerate due to the dominance of the vertical length scale and thus transit from the scale-dependent to the self-similar regime. Our results are congruent with previous theoretical studies^{11,12,22–48} and with available observations^{13–21} and elaborate extensive theory benchmarks for future experiments and simulations.

The group theory approach finds the multiplicity of nonlinear solutions for three-dimensional coherent dynamics of RT bubbles and RT spikes with variable acceleration and associates this multiplicity with the presence of the interfacial shear. Our approach can be further applied to find solutions for three-dimensional (3D) flows with other symmetries and for two-dimensional (2D) flows.^{12,37,49,56} We address the study of these problems to the future. As a brief summary, 3D flows tend to conserve isotropy in the plane, 3D highly symmetric dynamics is universal, and the dimensional 3D–2D crossover is discontinuous.^{12,37,49,56} Within the frame of the group theory approach, we can also directly link the scale-dependent dynamics and self-similar mixing, to be done in the future.^{11,22}

Our work studies the dynamics of ideal incompressible fluids. Our approach can be further applied to systematically analyze the influence of surface tension, viscosity, compressibility, and other factors on RT dynamics with variable acceleration. We address these studies to the future. Our work focuses on the acceleration-driven RT dynamics with the acceleration exponents $a>\u22122$. For $a<\u22122$, the dynamics is driven by the initial growth rate and is of Richtmyer–Meshkov (RM) type, whereas at $a=\u22122$, the transition occurs from RT type to RM type dynamics with varying acceleration strengths.^{22,23} We provide in the future the detailed analysis of 3D and 2D flows with various symmetries for acceleration exponents $a\u2264\u22122$.^{22}

Possible applications of our theory in high energy density plasmas include supernovae remnants, inertial fusion, and nano-fabrication, among others.^{5–9,57–60} For instance, the explosion of a core-collapse supernova is driven by a blast wave, which can be viewed as a variable strong shock.^{5,22} The blast-wave-induced acceleration is a power-law function of time, with the acceleration exponents $\u22122<a<\u22121$.^{50–54} According to our results, the scale-dependent interfacial dynamics driven by such accelerations is of RT type, and it is sub-exponential in the linear regime and decelerates in the nonlinear regime (Figs. 2 and 3). In the self-similar mixing, the dynamics can be of Richtmyer–Meshkov type and be sub-diffusive.^{22} The latter implies that the energy transport at small scales may occur via localization and trapping, thus influencing the synthesis of heavy mass elements in the supernova blast.^{22} Our results can thus be applied to explain the abundance of chemical elements and the richness of structures observed in supernovae remnants.^{5}

Our group theory approach finds that Rayleigh–Taylor dynamics can, in principle, be controlled, in excellent agreement with experiments in fluids and plasmas at Reynolds number up to $\u223c3.2\xd7106$.^{14–18,58} For inertial fusion, this suggests new opportunities for the control of plasma flows by means of variable acceleration and deterministic conditions.^{6,14,57–60} Particularly, in the inertial confinement fusion, it may be worth it to scratch the target in order to pre-impose appropriate deterministic conditions and to gain better control of fluid instabilities and the interfacial mixing.^{6,57–60} This may help the existing methods, which are focused on a fine polishing of ICF targets in order to fully eliminate RTI.^{6,57}

The results of the group theory approach can also be applied in nano-fabrication, by describing the dynamics and morphology of RT unstable interface in the scale-dependent regime and the scaling and properties, as well as the sensitivity to deterministic conditions in the self-similar regime.^{9,57–60}

To conclude, we studied the problem of RTI with variable acceleration by applying the group theory approach. We linked the interface dynamics to the interfacial shear, revealed the interfacial and multi-scale character of RT dynamics, and elaborated new theory benchmarks for future research.

## AUTHORS' CONTRIBUTIONS

S.I.A. contributed to conceptualization, formal analysis, methodology, investigation, project administration, resources, supervision, and writing, and K.C.W. contributed to investigation and preparation of figures.

## ACKNOWLEDGMENTS

The authors thank for the 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 Australian Government Research Training Program Scholarship at The University of Western Australia [K.C.W].

## DATA AVAILABILITY

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