On the Rayleigh-Taylor unstable dynamics of 3D interfacial coherent structures with time-dependent acceleration

Rayleigh-Taylor instability (RTI) occurs in a range of industrial and natural processes. Whereas the vast majority of existing studies have considered constant acceleration, RTI is in most instances driven by variable acceleration. Here we focus on RTI driven by acceleration with a power-law time-dependence, and by applying a group theoretic method find solutions to this classical nonlinear boundary value problem. We deduce that the dynamics is dominated by the acceleration term and that the solutions depend critically on the time dependence for values of the acceleration exponent greater than $-2$. We find that in the early-time dynamics, the RTI growth-rate depends on the acceleration parameters and initial conditions. For the later-time dynamics, we link the interface dynamics with an interfacial shear function, and find a continuous family of regular asymptotic solutions and invariant properties of nonlinear RTI. The essentially interfacial and multi-scale character of the dynamics is also demonstrated. The velocity field is potential in the bulk, and vortical structures appear at the interface due to interfacial shear. The multi-scale character becomes clear from the invariance properties of the dynamics. We also achieve excellent agreement with existing observations and elaborate new benchmarks for future experimental work.


I. INTRODUCTION
Rayleigh-Taylor instability (RTI) occurs whenever fluids of different densities are accelerated against their density gradient, and leads to intense interfacial Rayleigh-Taylor (RT) movement of the fluids [1][2][3] . The phenomenon plays a key roles in many processes, for instance, supernova explosions and inertial confinement fusion, and has received significant attention over the past seven decades. A reliable theory of RTI is required in order to expand our knowledge of non-equilibrium dynamics, to better understand RT-relevant phenomena, and to aid the development and improvement of processes in the areas of energy production and the environment, among many others 4,5 . Here we study the long-standing problem of RTI subject to a variable acceleration and use discrete group theory to solve the boundary value problem for early-and latetime RT evolution 4,6,7 . We directly link the interface dynamics and interfacial shear, discern its invariance properties for a broad range of acceleration parameters, and reveal the interfacial and multi-scale character of RT dynamics. We also elaborate extensive theory benchmarks that can be used in future experiments and simulations.
RT flows, whilst occurring in distinct physical circumstances, have similar evolution features 2,3,5-10 . RTI develops when the flow fields and/or the interface are slightly perturbed from their equilibrium state 1 . The interface is transformed into a composition of small-scale shear driven vortical structures and a large-scale coherent structure of bubbles and spikes, with a bubble (spike) being a portion of the light (heavy) fluid penetrating into the heavy (light) fluid 2,3,5-13 . Intense interfacial fluid mixing ensues with time 5,6,9-13 . a) Electronic mail: des.hill@uwa.edu.au b) Electronic mail: akbhowmi@andrew.cmu.edu c) Electronic mail: snezhana.abarzhi@uwa.edu.au RTI and RT mixing are challenging to study in experiments and simulations, and to investigate theoretically 4,5 . RT experiments in fluids and plasmas use advanced technologies to meet tight requirements for flow implementation, diagnostics and control [8][9][10] . RT simulations employ highly accurate numerical methods and massive computations to track unstable interfaces, capture small-scale processes and permit a large span of scales [11][12][13] .
Here we study RTI subject to a variable acceleration. Only limited information is currently available on RT dynamics under these conditions, suggesting the need for a systematic approach 4,22 . We consider accelerations with powerlaw time-dependence. These are important to study because they may result in new invariant and scaling properties of the dynamics 23 . They can be tuned to better model realistic environments and thus ensure practicality of our results 4,22,24-28 . We consider RTI in a 3D spatially extended periodic flow and apply group theory to solve the relevant boundary value problem, which involves boundary conditions at the interface and at the outside boundaries 7,20-22 . For early-time dynamics we identify the dependence of the RTI growth-rate on the acceleration parameters and initial conditions. For late-time dynamics, we directly link the interface dynamics to interfacial shear, find a continuous family of regular asymptotic solutions, and discover invariance properties of nonlinear RTI. The parameters of the critical, Atwood, Taylor and flat bubbles are identified, including their velocities, curvatures, Fourier amplitudes, and interfacial shear functions. We also reveal the essentially interfacial and multi-scale character of RT dynam-ics. The former is exhibited by the velocity field having intense fluid motion near the interface and effectively no motion in the bulk. The latter follows from the invariance properties of the dynamics set by the interplay of the two macroscopic length-scales -the wavelength and the amplitude of the interface.

II. THE METHOD OF SOLUTION
A. The governing equations The dynamics of ideal fluids is governed by conservation of mass, momentum and energy: where (x 1 , x 2 , x 3 ) = (x, y, z) are the spatial coordinates, t is time, (ρ, v, P, E) are the fields of density ρ , velocity v , pressure P and energy E = ρ(e + 1 2 v 2 ), where e is the specific internal energy 20 .
We consider immiscible, inviscid fluids of differing densities, separated by a sharp interface. It is required that momentum must be conserved at the interface and that there can be no mass flow across it. Hence the boundary conditions at the interface are where [· · · ] denotes the jump of functions across the interface; n and τ are the normal and tangential unit vectors of the interface with n = ∇θ |∇θ | and n · τ = 0; w = e + P ρ is the specific enthalpy; θ = θ (x, y, z,t), is a local scalar function, with θ = 0 at the interface and θ > 0 (θ < 0) in the bulk of the heavy (light) fluid, indicated hereafter by subscript h(l).
The heavier fluid sits above the lighter fluid and the entire system is subject to a time-dependent downwards acceleration field, directed from the heavy to the light fluid and is the power-law function of time g = (0, 0, −g) where g = Gt a . Here a is the acceleration exponent, and G > 0 is the acceleration pre-factor [17][18][19] . Their dimensions are [G] = ms −(a+2) and [a] = 1. This modifies the pressure field: P → P + ρgz. We assume that there are no mass sources and hence the boundary conditions There are two natural time scales in the problem, these are τ g = (kG) −1/(a+2) and τ 0 = 1 kv 0 , where v 0 is some initial growth rate and 1/k is the length scale. We consider here acceleration-driven RT dynamics a > −2. In this case, the former time scale is fastest, τ g τ 0 . We set the time scale to be τ = τ g . Time is t t 0 with t 0 τ, and the Atwood number

B. Large-scale coherent structures
These are arrays of bubbles and spikes periodic in the plane normal to the acceleration direction. At such scales the flow can be assumed to be irrotational at these large scales. We also assume that the fluids are incompressible and hence that the velocities are expressible in terms of scalar potentials Φ h (x, y, z,t) and Φ l (x, y, z,t). Because the fluids are ideal these will be harmonic. That is, We focus on bubbles propagating in the z-direction. For convenience our calculations are performed in the frame of reference moving with velocity v(t) in the z-direction, where v(t) = ∂ z 0 /∂t and z 0 are the velocity and position of the bubble in laboratory reference frame. The interface shape is θ (x, y, z,t) = z−z * (x, y,t) = 0, and the interface conditions are then The vertical far-field boundary conditions are C. The dynamical system The periodic nature of the large-scale coherent structure can be accommodated by appealing to the theory of discrete groups [5][6][7] . We first identify groups enabling structurally stable dynamics. These are, for example, the group p6mm for hexagonal symmetry, p4mm for square symmetry, p2mm for rectangular symmetry. The relevant symmetry group (in our case, p4mm) dictates a specific Fourier series (an irreducible representation of the group) which can be used to solve the nonlinear boundary value problem Eqs. 4,5. We then make spatial expansions in the vicinity of the tip of a bubble. This approach reduces the governing equations to a dynamical system of ordinary differential equations in terms of interface variables and Fourier moments [5][6][7][19][20][21][22] .
For three-dimensional flow with square symmetry, the potentials are where α mn = √ m 2 + n 2 , m and n are integers, k = 2π λ is the wavenumber, Φ mn andΦ mn are the Fourier amplitudes for the heavy and light fluids respectively, and Φ 00 =Φ 00 = 0. Symmetry requires that Φ mn = Φ nm andΦ mn =Φ nm .
In order to examine the local behavior of the interfacial dynamics in the vicinity of the bubble tip, we expand the interface function in a power series about (x, y) = (0, 0). In the moving frame of reference, this is where ζ i j (t) = ζ ji (t) due to symmetry, ζ (t) = ζ 10 (t) is the the principal curvature at the bubble tip, and N = i + j is the order of the approximation. To lowest order (that is, The Fourier series and interface function are substituted into the interface conditions and the resulting expressions expanded as Taylor series. This yields a system of ordinary differential equations for Φ m (t),Φ m (t) and ζ i j (t). We may express the potentials in terms of moments M a,b,c (t) = ∑ mn Φ mn (t)(mk x ) a (nk y ) b α c mn and their tilde equivalents. We note that by symmetry, M a,b,c = M b,a,c and M a+2,b,c + M a,b+2,c = M a,b,c+2 and similarly forM. The vertical far-field conditions give M 0 = −M 0 = −v(t). For N = 1, we abbreviate the series to second order in x and y, and first order in z since z * (x, y,t) is quadratic in x and y. The interface conditions becomė where In this regime, the system can be linearised and only one harmonic is needed, that is, the moments retain only one Fourier amplitude. The initial conditions at time t 0 are the initial curvature ζ 0 = ζ 1 (t 0 ) and velocity v 0 = v(t 0 ).
For a broad class of initial conditions, integration of the governing equations is a challenge. The solution can be found 2,3,22 when the amplitude of the initial perturbation is small τk|v 0 | 1, and the interface is nearly flat |ζ 1 /k| 1. The system reduces tȯ 2. Later time, t − t 0 τ In the later-time situation, the behaviour is nonlinear and multiple harmonics must be retained. We find asymptotic solutions for the relevant equations and determine their stability. To leading order in time, regular asymptotic solutions will have the following time-dependence: We investigate the stabiity of the asymptotic solutions by including the perturbations where ξ (t) is to be determined. Since there are in total four equations (three conservation equations and the auxilliary condition M 0 = −M 0 ), we consider perturbations of four of the variables (three Fourier amplitudes and the bubble curvature) to analyze the solution stability. The results are insensitive to which harmonic we choose to leave unperturbed. The interface conditions require thatξ = βt a 2 ξ , where β is to be determined. Perturbations are stable if Re[β ] < 0, otherwise unstable. The solution of this differential equation is For a = 0 we have exponential solutions and the family of solutions are exponentially stable/unstable. For a = −2 we have the power-law solution ξ (t) = t β and the family of solutions are asymptotically (but not exponential) stable/ unstable. Eq. 9 now yields a quadratic equation for β which does not depend on the value of a, but is algebraically unwieldy.

B. The later-time regime
At later times, spikes are singular (the singularity is finitetime), whereas bubbles are regular [5][6][7] . For t τ, higher order harmonics are retained in the expressions for the moments, and regular asymptotic solutions can be derived. For N = 1, the first two harmonics are retained and we arrive at a one-parameter family of solutions (as there are four equation in five unknowns). We choose the bubble curvature ζ to parametrize the family. Substitution of the asymptotic forms Eq. 11 into Eqs. 8 and 9, employing a dominant balance agrument and solving the resulting set of equations leads to the solution (15) which is valid for σ ∈ (0, σ cr ) where σ cr = 3 8 with corresponding ζ cr = − 3 8 k. Fig. 1 shows the bubble tip velocity as a function of the bubble curvature. We observe that the bubble tip velocity is larger for larger values of the Atwood number A, and occurs at a steeper curvature. The Fourier amplitudes are Solutions for N > 1 can likewise be calculated. These solutions converge for increasing N and in each case the lowest order harmonics are dominant. Figs. 2 and 3 demonstrate that the second Fourier amplitude is much smaller than the first for σ < σ cr . At N = 1, the solution with maximum growth rate is a curved bubble, having a curvature which depends on the Atwood number. The curvature σ max which maximizes the velocity satisfies This is a universal relation between the curvature and the maximum velocity.

C. The effect of shear
The multiplicity of these solutions is also due to the presence of shear at the interface. We define shear function Γ to be the spatial derivative of the jump in the tangential velocity across the interface. We find that in the vicinity of the bubble tip it is Γ =M 1 − M 1 . Specifically, and is a strictly monotone function of σ , rising from Γ = 0 at σ = 0 to Γ = Γ max at σ = σ cr . Fig. 4 shows the interface shear as a function of the bubble curvature. We note that the shear function is larger for larger values of the Atwood number A, and in each case tends towards a constant value as the curvature increases. Fig. 5 shows the variation of the bubble tip velocity with the interface shear function. We note that the velocity achieves a maximum at some point and that this peak occurs at higher values of the shear function for larger Atwood numbers. The peak velocity itself is also larger for larger Atwood numbers. We particularly note that very soon after the velocity curve reaches its peak, it drops sharply. This rapid change presents significant problems for numerical simulations of RTI.

The flat bubble
The solution corresponding to a flat bubble is simply

The Atwood bubble
The fastest member of the family we refer to as the 'Atwood bubble' to emphasise its complex dependence on the Atwood number. The solution is where and the Fourier amplitudes are given by Eq. 16. We note that theΓ A (V A ) function is very nearly linear for V A ∈ (0, 1). In the limit A → 1, the solution is and the Fourier amplitudes Φ 10A1 , Φ 20A1 , Φ 10A1 ,Φ 20A1 are, respectively, In the limit A → 0, the solution is (21) and the Fourier amplitudes are

The 'Taylor' bubble
We refer to this bubble as a 'Taylor bubble' since its curvature is as in 2 except for a difference in the wavenumber value. The solution is

The critical bubble
The solution is (23) and Fourier amplitudes are E. Stability AkG. In the limits A → 1 − and A → 0 stability is achieved when the curvature reaches respectively. All bubbles are stable for A ≥ 1 24 and hence the Taylor and critcal bubbles are stable. The Atwood bubble is stable when A = 1, having σ = − 4 3 √ kG, −2 √ kG. The Atwood bubble for A = 0.1 and A = 0.01 are also stable and it would appear that the Atwood bubble is stable for any value of Atwood number A at N = 1. The N > 1 analysis is to be presented elsewhere. As is the case for a = 0, we expect that the stability interval will narrow sharply for N > 1.  [5][6][7][19][20][21]30,31 . Thus, our results excellently agree with existing observations [8][9][10][11][12][13]22 . Our theory is focused on large-scale dynamics and presumes that interfacial vortical structures are small-scale. This assumption is applicable for fluids very different densities and with finite density ratios. For fluids with very similar densities 0 ≈ A 1 other approaches can be employed [5][6][7][20][21][22] .

F. The velocity field
By accurately accounting for the interplay of harmonics and systematically connecting the interfacial velocity and shear for a broad range of acceleration parameters, we have found that RT dynamics is essentially interfacial: It has intense fluid motion in the vicinity of the interface, effectively no motion away from the interface and shear-driven vortical structures at the interface. The velocity is potentisl in the bulk of each fluid. This velocity pattern is observed in experiments and simulations, demonstrating excellent agreement with our results [8][9][10][11][12][13]22 .

IV. DISCUSSION
We have found solutions for acceleration driven RTI in a 3D spatially extended periodic flow in both the early-time regime and later-timer regime. In the early-time regime the dynamics is faster than exponential for a > 0 and slower than exponential for −2 < a < 0. In the later-time regime, bubbles are accelerated for a > 0, steady at a = 0, and decelerated for −2 < a < 0. Since the dependencies are given by standard formulae, we can compare various acceleration exponents. The dynamics for small exponents are usually a diagnostic challenge becaue it is slow. Specific examples are fusion, supernovae and nano-fabrication. The dynamics for fast exponents can be easily diagnosed, and the deduced properties can be applied to the slow-dynamics case. Hence we have obtained an important practical result.
The nonlinear bubble velocity and shear, when rescaled as kτ(t/τ) −a/2 v and τ(t/τ) −a/2 Γ, depend only on the interface morphology and flow symmetry. Hence, by analyzing properties of RT bubbles for fast dynamics and large exponents a > 0, we can obtain properties of those for slow dynamics and small exponents −2 < a < 0 . These are especially convenient for studies of RTI in high energy density plasmas in astrophysics and fusion, where RTI is driven by an explosion or an implosion with blast wave acceleration exponents −2 < a < −1 in these cases 22,23,29 .
Our analysis also elaborates upon diagnostic quantities which have not been discussed before. These are the velocity and pressure fields, the interface morphology and bubble curvature, the interfacial shear and its link to the bubble velocity and curvature, the spectral properties of the velocity and pressure, along with the interface growth and growth-rate. By determining the dependence of these quantities on the density ratio, flow symmetry, and acceleration exponent and strength, and by identifying their universal properties, and comparing all of these to data obtained for real fluids, we can further advance our knowledge of RT dynamics in realistic environments, thereby achieving a better understanding of RT relevant processes and also improve methods of numerical modeling and experimental diagnostics of interfacial dynamics in fluids, plasmas and other materials.
Our techniques can also be used to study Richtmyer-Meshkov instability (that is, a < −2). The results obtained are very different from those of the Rayleigh-Taylor instability (that is, a > −2).
We note also that the expressions obtained in the constant acceleration a = 0 case are unexpectedly very similar in form to those of the −2 < a < 0 case.