Confinement-induced stabilization of the Rayleigh-Taylor instability and transition to the unconfined limit

Sufficient confinement can completely suppress the Rayleigh-Taylor instability between two density-inverted miscible fluids.


Linear-Stability Analysis with Boundaries Single boundary
To probe how the most unstable wavelength is affected by confinement we incorporate two walls into the linear-stability analysis of Chandrasekhar (9), as illustrated in Fig. S1A. Two fluids with densities ρ 1 and ρ 2 , and dynamic viscosities η 1 and η 2 form an interface located between two walls a distance b apart at a location cb, with 0 < c < 1. To begin we set c = 0.5 to match the majority of the simulations reported in the manuscript. Note that we do not include the effect of mass diffusion in this analysis. We here show the dispersion relation and the wavelengths that result from this linear-stability analysis; details of the calculation of the dispersion relation are provided in the next section.
For now we restrict ourselves to the case of ν 1 = ν 2 , where ν = η/ρ. Fig. S1B shows the dispersion relation for varying b, where n is the instability growth rate and k is the wavenumber of a disturbance at the fluid-fluid interface. The change in color from black to red denotes a decrease in b. The location of the maximum growth rate, where ∂ k n = 0, indicates the most unstable wavenumber that will be selected by the system. The magnitude of the growth rate decreases and the most unstable wavenumber increases as we make b smaller. The most unstable wavelength, λ = 2π/k, increases linearly with the plate spacing at small b and plateaus at larger b, as shown in Fig. S1C. In the linear regime, λ ≈ 1.28b.
We compare λ to the most unstable wavelength λ ∞ selected when there is no boundary (9).
In this unconfined case, for ν 1 = ν 2 , the dispersion relation has the approximate form (33,34) where A = (ρ 2 −ρ 1 )/(ρ 1 +ρ 2 ) is the Atwood number. This yields the most unstable wavelength We use λ ∞ to rescale both λ and b in Fig. S1C. The plateau wavelength at large b is precisely λ ∞ ; the transition from the linear regime to open space limit is reached when 1.28b ≈ λ ∞ .
By varying the value of c we probe how the proximity of the interface to either boundary affects the wavelength. In exploring this we restrict ourselves to the confined regime, b < λ ∞ .
The wavelength exhibits a maximum at c = 0.5 and decreases linearly as c approaches either 0 or 1, as shown in Fig. S1D. These two limits correspond to the limit where only one boundary matters and the interface does not feel the effect of the second wall. In this limit λ ≈ 3.7cb, where cb is the distance to the nearest boundary.
As b decreases, the growth rate for the most unstable mode also decreases, as seen in Fig. S1B. The timescale for the growth of the instability thus increases and eventually becomes comparable to the diffusive timescale at small enough b. Comparing these two timescales allows us to estimate the critical plate spacing b c below which the instability no longer occurs. We note that the growth rate for the full theory is always bounded by the approximate form from Hide (33), seen in Eqn. (1). At large k beyond the most unstable wavenumber, the growth rate is bounded by Since the most unstable wavelength is λ ≈ 1.28b, we can estimate the growth rate n c for this most unstable mode to be The characteristic timescale for diffusion is b 2 /4D. By equating these two timescales we can estimate at what value of b they become comparable; this occurs when If we interpret this length scale as the critical gap, This value is of the same order as the limiting case of b c ≈ 15b * shown in the main text. To reach a quantitative match, mass diffusion must be included into the linear-stability analysis.

Derivation of the dispersion relation
Here we show how to incorporate boundaries into the linear-stability analysis for the Rayleigh-Taylor instability. To reduce the amount of algebra and clearly identify the effect of incorporating a boundary we first show, in detail, the derivation for a single wall. This is followed by stating the result for two walls, adapting the same procedure as for one wall.

Dispersion matrix derivation: single boundary
To get the dispersion matrix we closely follow the work of Chandrasekhar (9) and Zetina et al. (31). We have the linearized Navier-Stokes equation: where w is the z-component of the velocity field. The solutions of this equation are of the form: where q 2 = k 2 + n/ν. The problem we address here is different from that in reference (9), in that we have a boundary located close to the unstable interface. The wall is at a distance H from the free boundary located at z = 0, see Fig. S2A. There are a few boundary conditions we need to satisfy. The first is that on any wall we require no-slip for the velocity field. This gives us two boundary conditions, w = 0 and ∂ z w = 0. We require this at z = H and z = −∞. From these boundary conditions we have that with index 1 and 2 referring to the lower and upper fluids. Using the boundary conditions on z = H we get By multiplying equation (12) by k or q 2 and subtracting equation (13) we solve for C 2 and D 2 in terms of A 2 and B 2 . This gives where the β's are introduced for brevity.
To solve for the dispersion matrix we consider the boundary conditions on the unstable fluid-fluid interface. From Chandrasekhar (9) we have the following boundary conditions at z = 0: where [w] 0 and [∂ z w] 0 are the common value at the interface (for w this could be the average value of w at the interface; but any choice of w 1 , w 2 , or the average are consistent). These four boundary conditions correspond to continuity of velocity, normal stress, tangential stress, and the kinematic boundary condition respectively.
By using the ansatz solutions for w 1 and w 2 and the substitutions for C 2 and D 2 we get the following equations from each boundary condition: We define the variables: Using Eqns. (20) -(23) we look for non-trivial A 1 , B 1 , A 2 , B 2 that satisfy them. We express these equations in matrix form such that where M is To ensure that the solution for C is non-trivial we need to find an M satisfying det(M) = 0.
Calculating this determinant yields the dispersion relation that relates k and n. The calculations of the dispersion relation are done numerically.

Dispersion matrix derivation: two boundaries
We follow the same process for calculating the matrix M for the case of two boundaries, where the velocity for w 1 now has additional terms. We define the height above the interface as H 2 and the distance below the interface as H 1 , as shown in Fig. S2B, and get the matrix with H → H 2 in the β terms and the following definitions for the γ terms (note that they are the same as the β terms but with indices switched from 2 to 1 for the fluid identification)