\section{Introduction}

Rayleigh-Taylor instability (RTI) occurs when the interface between two fluids of different density is subjected to a normal pressure gradient.
Geometric complexity is an essential feature of (RTI), especially in three space dimensions. When a fluid undergoes RTI, vorticity is deposited on the  material interface  separating the two fluids. This interface vorticity then initiates the Kelvin-Helmholtz instability (KHI), which causes the interface to roll up into extremely complex shapes. Accurately simulating such behavior requires maintaining a sharp discontinuity in density across the interface, while at the same time faithfully representing the interface geometry. The rectangular cells of Eulerian differencing schemes are ill-suited to the highly curved geometry of RTI problems, and often introduce spurious KHI spirals which obscure the form of the primary instability. Indeed, a look at the bewildering variety of behaviors observed in the excellent comparative study of Liska and Wendroff \cite{LiWe2003} demonstrates that high-order methods such as WENO and PPM, set on structured Cartesian meshes, can be susceptible to spurious small-scale KHI.

To improve accuracy, methods such as marker-and-cell \cite{Da1967} and volume-of-fluid \cite{HiNi1981} use either marker particles or a marker function in addition 
to an Eulerian grid to track the interface. Similarly, front-tracking methods, which use a triangular mesh that advects with the fluid velocity, have shown considerable 
success \cite{TrUn1991}. When the fluids are inviscid, incompressible, and irrotational in their respective bulks, the 3D Euler equations can be written as boundary 
integral equations involving only the interface position and the vortex sheet density on the interface. In three dimensions, these nonlocal integral equations were first derived and simulated by Baker, Meiron, and Orszag in \cite{BaMeOr1984}, following the earlier derivation by Birkhoff \cite{Bi1962} and Rott \cite{Ro1956} for  vortex sheets in 2D flows.

Discretization of Birkhoff's 2D boundary integral equation for vortex sheets using a trapezoidal rule yields an even earlier method, known as the 
\emph{point vortex method}, where a vortex sheet is approximated by finitely many point vortices. The earliest simulations of vortex sheet motion, by 
Rosenhead \cite{Ro1931}, computed the motion of twenty-five point vortices to approximate a vortex sheet; however, doubts were raised by 
Birkhoff and Fisher about 
the convergence of the point vortex method to solutions of the 2D Euler equations, and they used the Hamiltonian nature of systems of points vortices to predict eventual chaotic motion \cite{BiFi1959}. Subsequent attempts to reproduce Rosenhead's vortex sheet rollup at higher resolutions verified Birkhoff and Fisher's suspicions \cite{Ta1964}. A number of numerical approximations
of this problem have been attempted, but the most successful numerical algorithm has been Krasny's desingularization method \cite{Kr1986}\cite{Kr1986b}, based on the \emph{vortex blob method} of Chorin \cite{Ch1973}. The convergence properties of Chorin's approximation were studied in a series of papers by Hald \cite{Ha1978}\cite{Ha1979}\cite{Ha1987}, and the convergence properties of Krasny's approximation were studied by Beale and Majda \cite{BeMa1982}\cite{BeMa1985}.


The boundary integral formulation of the Euler equations, relying on a parametrization of the material interface, has several attractive features for simulating RTI. It is 
naturally adaptive, concentrating grid points near regions of physical interest, it avoids the numerical diffusion of grid-based methods, and it involves fewer independent 
variables. However, these benefits come at a significant computational cost, due to the nonlinear and nonlocal form of the time derivatives in this formulation. The use of standard time-stepping
algorithms for evolution equations in which the time-derivatives are linear and isolated cannot be used.  Rather, 
simulating such equations requires
 either the solution of a nonlinear Fredholm integral equation of the second kind for the time-derivatives, or the use of an iterative predictor-corrector method to evolve forward in time. Additionally, the nonlocal expression requires $O(N^2)$ operation to evaluate (where $N$ is the number of interface points) making such predictor-corrector methods extremely expensive. 
 

 Herein, we introduce three interface models for RTI problems which we shall refer to as $z$-models, with $z=(z^1,z^2,z^3)$ denoting a parameterization of the 
material interface.  These three $z$-models
are nonlocal evolution equations for the interface parameterization $z$ and the vorticity restricted to the interface $\mu$ (whose precise definition will be given below), 
and will be presented in what we term low-order, medium-order, and high-order variants.  The order of these three models coincides with the amount of 
nonlocality each model retains from the original Euler system.  Specifically, 
 with the aim of preserving the geometric complexity of RTI, we do not obtain such model equations using typical asymptotics in either amplitude or slope, but rather, we  use an asymptotic expansion  in the nonlocality of the  interaction potential between interfacial particle positions. 
 The
 resulting $z$-models are then subjected to single-mode and (random) multi-mode RTI numerical tests and those results are compared to classical RTI experiments.
It is verified that the $z$-models reproduce experimental results on single-mode RTI and match theoretical growth rates for Rayleigh-Taylor mixing layers, all while achieving the same level of accuracy as Eulerian finite differencing with roughly 600 times less computational effort \cite{Reisner}.


\subsection{Outline}
In Section 2, we write the irrotational and incompressible Euler equations as a system of singular integral equations for the fluid interface parameterization
and the vorticity measure on the interface.   By expanding the leading order terms in an asymptotic expansion of the interaction kernel, we derive the
 $z$-models for the interface motion and describe the numerical discretization that is used for simulations.


In Section 3, the $z$-models are used for numerical simulations of a variety of classical Rayleigh-Taylor and Kelvin-Helmholtz test problems.  These
problems are initiated using both
 single-mode initial data (Section \ref{sec::single-mode-data}) as well as random multi-mode initial data (Sections \ref{sec::rocketrig}--\ref{sec::tiltedrig}). 
 Because of the speed of our simulations, we are able to consider ensemble averages of numerous runs with with random initial data and show that such
 ensembles are able to model viscous mixing layers with optimal mixing rates (Section \ref{sec::mixing}).   The $z$-models are not restricted to horizontally-periodic 
 interface motion, and can simulate the motion of fluid interfaces that are represented by closed surfaces (such as spheres).   We simulate the motion of
a rising bubble, and obtain  qualitative agreement with high-speed photography experiments  (Section \ref{sec::bubble}).

\section{Interface models}\label{sec::models}

We consider two (open) volumes of inviscid,  incompressible, irrotational fluid $\Omega^+(t)$ and $\Omega^-(t)$, separated by a time-dependent evolving
 material interface $\Gamma(t)$ (the free surface),  where
$t$ denotes time.  The  evolution of the fluid  is restricted to a time interval $0 \le t \le T$, and is modeled by the incompressible and irrotational Euler equations,
which can be written as
\begin{align} 
\left.
\begin{array}{r}
\rho ^\pm( \partial_t u^\pm + u^\pm \cdot \nabla u^\pm) + \nabla p^\pm = - g \rho^\pm e_3\\
\partial_t \rho^\pm + u^\pm \cdot \nabla \rho^\pm =0
\end{array}
\right\}
 \ \text{ in } \ \Omega^\pm(t) \,, \ \ \ 0< t\le T \,, 
 \label{euler00}
\end{align} 
with the constraints that
\begin{align*} 
\operatorname{div} u^\pm =0 \ \text{ and } \ \operatorname{curl} u^\pm =0 \ \text{ in } \ \Omega^\pm(t) \,, \ \ \ 0\le t\le T \,.
\end{align*} 
The initial interface $\Gamma(0)$ is specified and hence so too are the initial domains $\Omega^\pm(0)$, and initial conditions for
velocity and density are given:
\begin{align*} 
u^\pm(x,0) = \mathring u^\pm(x) \,, \ \ \  \rho^\pm (x,0)= \mathring \rho^\pm (x)   \ \text{ in } \ \Omega^\pm(0) \,.
\end{align*} 
Here, $u^\pm$ denotes the velocity  and $p^\pm$ denotes the pressure in $\Omega^\pm(t)$, while $g$ represents the (assumed) constant gravitational
acceleration, and the unit vector $e_3$ is $(0,0,1)$.   We assume constant  initial density functions 
$\mathring\rho^\pm$ in
$\Omega^\pm(0)$, which implies that $\rho^\pm$ remain constant functions in $\Omega^\pm(t)$.
Our focus is on flows in
which $\Gamma(t)$ denotes a vortex sheet, and hence we supplement the Euler equations with the following jump conditions on $\Gamma(t)$:
\begin{subequations} 
\label{jump}
\begin{align} 
\left.
\begin{array}{r}
\jump{p} =0 \\
\jump{u \cdot n} =0 
\end{array} 
\right\} \ \text{ for } \ 0\le t\le T \,,
\label{jump-unp}
\end{align} 
where  $n$ denotes the  unit normal to $\Gamma(t)$, and $\jump{f} = f^+ - f^-$ on $\Gamma(t)$.     We assume that $\mathring\rho^+ \neq \mathring\rho^-$ so that
\begin{align} 
\jump{\rho} \neq 0 \ \text{ on } \ \Gamma(t) \ \text{ for } \ 0\le t\le T \,. \label{jump-rho}
\end{align} 
Because of \eqref{jump-unp} and \eqref{jump-rho}, even if $\mathring u^+ = \mathring u^-$ on $\Gamma(0)$, the 
tangential component of velocity becomes discontinuous across $\Gamma(t)$ and we have that 
\begin{align} 
\jump{u \cdot \tau_\alpha} \neq 0 \ \text{ on } \ \Gamma(t)  \ \text{ for } \ 0< t\le T \,, \ \ \alpha=1,2\,. \label{jump-utau}
\end{align} 
\end{subequations} 
Here,  $\tau_1$ and $\tau_2$ are  the unit tangent vectors to $\Gamma(t)$, chosen such that $(\tau_1,\tau_2, n)$ form a right-handed orthonormal basis.   To complete the description of the dynamics,
the motion of the interface $\Gamma(t)$ is governed by the normal component of the fluid velocity.  Letting
$\mathcal{V}(\Gamma(t))$ denote the normal speed of the interface,
\begin{align} 
\mathcal{V} (\Gamma(t)) = u \cdot n \,.  \label{Gamma-speed}
\end{align} 





Let $(x^1,x^2,x^3)$ denote the standard Euclidean coordinates on $\mathbb{R}^3$, and let $(s^1,s^2)$ denote coordinates on $ \mathbb{R}^2  $, used to
parameterize $\Gamma(t)$.  Specifically, the time-dependent interface
$\Gamma(t)$ is parametrized by a smooth function $z:\mathbb{R}^2 \times [0,T] \to\mathbb{R}^3$ and
\begin{align} 
\Gamma(t) = \{ z(s^1,s^2,t) \ : (s^1,s^2) \in \mathbb{R}^2  \,, \ t \in [0,T]\} \,, \ \ \ z=(z^1,z^2,z^3) \,. \label{z-param}
\end{align} 
We use Latin indices for coordinates in Euclidean space, and Greek indices for coordinates on $\Gamma(t)$, and we will apply Einstein's summation convention without further comment. Euclidean space is endowed with the standard diagonal metric $\delta_{ij}$, and the induced metric on $\Gamma(t)$ is given by
$$h=h_{\alpha\beta}ds^\alpha\otimes ds^\beta\,,\qquad\qquad h_{\alpha\beta}=\partial_\alpha z\cdot\partial_\beta z\,.$$
We set
$$|h|=\det h\,,\qquad\qquad \sqrt h=\sqrt{\det h}\,. $$
With respect to the parameterization \eqref{z-param}, the  time-dependent unit normal is given by
 $$n= \frac{\partial_1z\times\partial_2z}{\sqrt h} \,. $$

As noted above, we assume that  the vorticity vanishes in the open sets $\Omega^+(t)$ and $\Omega^-(t)$,
 but on the material interface $\Gamma(t)$,  the tangential component of velocity experiences a jump discontinuity, resulting in a vorticity measure 
 $\omega(s,t)$ concentrated 
 on the interface. The velocity in $\Omega^+(t)\cup \Omega^-(t)$ is computed from this vorticity measure using the Birkhoff-Rott integral
\begin{equation}\label{BR1}
	u(x,t)=\frac{1}{4\pi}\iint_{\mathbb{R}^2}\omega(s,t)\times\frac{x-z(s,t)}{|x-z(s,t)|^3}ds \,,
\end{equation}
where 
$$\omega=\omega^\alpha\partial_\alpha z$$ 
and $\omega^\alpha$, $\alpha=1,2$ denote the components of the vorticity measure with respect to the basis $\partial_\alpha z$.   To be more precise, 
$\omega(s,t)=\omega(s^1,s^2,t)$
denotes the the vortex sheet {\it density} and is often referred to as the {\it amplitude of vorticity}  or simply the vorticity measure on $\Gamma(t)$.
When evaluated on $\Gamma(t)$ (or equivalently, along the parameterization $z(s,t)$), the Birkhoff-Rott integral exists in the principal value sense, 
and gives the average velocity 
\begin{align} 
\bar u=  \tfrac{1}{2} (u^++u^-) \circ z \,. \label{baru-def}
\end{align} 
In particular, by choosing the Lagrangian parameterization $z(s,t)$ such that $\partial_t z(s,t) = \tfrac{1}{2} (u^++u^-) ( z(s,t),t)$ with $z(s,0)=s$, we have that
\begin{equation}\label{BR2}
	\partial_tz(s,t)=\bar u(s,t):=\frac{1}{4\pi}\iint_{\mathbb{R}^2}  \omega(s',t)\times\frac{z(s,t)-z(s',t)}{|z(s,t)-z(s',t)|^3}ds' \,.
\end{equation}
Let us remark that according to \eqref{Gamma-speed}, it is only the normal component of $u$ that determines the shape of $\Gamma(t)$, and there exists an 
infinite-dimensional tangential reparameterization symmetry of the interface which does not change its shape.   This means that a different parameterization 
$\mathsf{z}(s,t)$ 
such that $\partial_t \mathsf{z} \cdot n ( \mathsf{z},t)  = \bar u(\mathsf{z},t) \cdot  n(\mathsf{z},t)$ 
but 
 $\partial_t \mathsf{z} \cdot \tau_\alpha(\mathsf{z},t) \neq \bar u(\mathsf{z},t) \cdot  \tau_\alpha(\mathsf{z},t)$ would provide the same shape for $\Gamma(t)$,
but would alter the distribution of particles along the interface.  For our purposes, the Lagrangian parameterization \eqref{BR2} is a convenient choice, as the resulting evolution equation \eqref{fe2} takes a particularly simple form. 

We may write the vorticity measure $\omega$ on $\Gamma(t)$ in terms of the velocity jump
\begin{align} 
w=(u^+-u^-)\circ z \,, \label{u-jump}
\end{align} 
in the form
\begin{equation}\label{vort}
	\omega=(w\cdot\partial_2z)\partial_1z-(w\cdot\partial_1z)\partial_2z \,,
\end{equation}
so that $w\cdot\partial_2z= \omega^1$ and $w\cdot\partial_1z= - \omega^2$.
Some of the difficulty of vortex methods stem from the computation of the Birkhoff-Rott velocity; specifically, 
 a straightforward numerical quadrature results in chaotic motion of the interface \cite{Kr1986}. Let us describe our method of evaluation. Equation \eqref{BR2} may be rewritten as
\begin{align} 
\bar u(s,t)=\iint_{\mathbb{R}^2  }\omega(s',t)\times\nabla G(z(s,t)-z(s',t))ds' \,, \label{u-bar}
\end{align} 
where $G(x)=\tfrac{1}{4\pi|x|}$ is the fundamental solution of the Laplacian in $\mathbb{R}^3$. Note that \eqref{BR1} is defined everywhere in the bulk 
$\Omega^+(t)\cup\Omega^-(t)$ and has  arguments $(x^1,x^2,x^3,t)\in\mathbb{R}^3 \times [0,T]$, while \eqref{BR2} is well-defined
on the interface $\Gamma(t)$ and has arguments $(s^1,s^2,t)\in\mathbb{R}^2\times [0,t]$. To avoid the singularity in the integral \eqref{u-bar}, we replace 
$G$ by a desingularized function
$$G^\epsilon(x)=\frac{1}{4\pi\sqrt{\epsilon^2+|x|^2}},$$
where $0<\epsilon<1$ and $G^\epsilon$ converges weakly to $G$ as $\epsilon\to 0$. This results in the desingularized Birkhoff-Rott velocity
\begin{equation}\label{ubarepsilon}
	\bar u^\epsilon(s,t)=\frac{1}{4\pi}\iint_{\mathbb{R}^2}\omega(s',t)\times\frac{z(s,t)-z(s',t)}{\big(\epsilon^2+|z(s,t)-z(s',t)|^2)^{3/2}}ds'.
\end{equation}
This particular desingularization was chosen for algebraic simplicity and falls under the well-known umbrella of \emph{vortex blob methods} (see, for example,
\cite{Kr1986b}).

 These are methods where the standard Green's function, satisfying $\Delta G=\delta$, is replaced with a desingularized Green's function  satisfying 
$$\Delta G^\epsilon(x)=\psi^\epsilon(x):=\frac{1}{\epsilon^3}\psi\left(\frac{x}{\epsilon}\right),$$
where $\psi$ is a smooth, nonnegative function with most of its mass near zero and $\int_{\mathbb{R}^3}\psi=1$. This kind of desingularization has the effect of replacing a singular vorticity distribution $\omega(x)=\omega_0\delta(x-x_0)$ with a smooth `vortex blob' of finite size $\omega(x)=\omega_0\psi^\epsilon(x-x_0)$. For our choice of desingularization, 
\begin{align} 
\psi(x)=\tfrac{3}{4\pi}\tfrac{1}{(1+|x|^2)^\frac{5}{2}}.
\label{kernel-desing}
\end{align} 
With this desingularization employed, we may evaluate $\bar u^\epsilon$ with any standard 2D quadrature method. This  allows us to obtain arbitrarily high-order spatial discretizations.

\subsection{The irrotational and incompressible  Euler equations using potential functions}
Thanks to the irrotationality of the two fluids, we may write the velocities $u^\pm$ in $\Omega^\pm(t)$  in terms of velocity potentials,  $u^\pm = \nabla \varphi^\pm$,
where $\varphi^\pm$ are governed by Bernoulli's law,
\begin{equation}\label{bern}
	\begin{cases}
	\partial_t\varphi^+ +\tfrac{1}{2}|u^+|^2+gx^3=-p^+/\rho^+&\text{in }\Omega^+(t)\\
	\partial_t\varphi^- +\tfrac{1}{2}|u^-|^2+gx^3=-p^-/\rho^-&\text{in }\Omega^-(t)
\end{cases} \,,
\end{equation}
and where (we recall that) $p^\pm$ denotes the pressure functions, $\rho^\pm$ denotes the density functions,  and  $g$ denotes the gravitational acceleration. 
We assume that the two fluids have infinite extent, so that $\Omega^+(t)\cup \Gamma(t)\cup \Omega^-(t)=\mathbb{R}^3$. The jump conditions on $\Gamma(t)$
 associated to  \eqref{bern}  are given by \eqref{jump}.
 
 To determine how the potential jump $\varphi^+-\varphi^-$  varies along $\Gamma(t)$, we compose \eqref{bern} with the parameterization $z$ and apply the chain rule:
\begin{align*}
	\partial_t((\varphi^+-\varphi^-)\circ z)&=(\nabla\varphi^+\circ z-\nabla\varphi^-\circ z)\cdot\partial_tz+\partial_t(\varphi^+ -\varphi^-)\circ z=w\cdot\bar u
	+\partial_t(\varphi^+-\varphi^-)\circ z \,, \\
	\partial_t((\varphi^++\varphi^-)\circ z)&=(\nabla\varphi^+\circ z+\nabla\varphi^-\circ z)\cdot\partial_tz+\partial_t(\varphi^+ +\varphi^-)\circ z=2|\bar u|^2
	+\partial_t(\varphi^++\varphi^-)\circ z \,. 
\end{align*}
Here,  $w=u^+-u^-$ is the velocity jump on $\Gamma(t)$. Note that $\bar u$, defined by \eqref{u-bar}, and $w$  are well-defined on $\Gamma(t)$, and 
are thus functions of $(s^1,s^2)$ and $t$.
We introduce the {\it Atwood number}
$$A=\tfrac{\rho^+-\rho^-}{\rho^++\rho^-}\,.$$
By using \eqref{bern}, we write
\begin{equation}\label{presdiff}
	\left(\frac{p^+}{\rho^+}-\frac{p^-}{\rho^-}\right)+A\left(\frac{p^+}{\rho^+}+\frac{p^-}{\rho^-}\right)=(A+1)\frac{p^+}{\rho^+}+(A-1)\frac{p^-}{\rho^-}=2\frac{p^+-p^-}{\rho^++\rho^-}=0.
\end{equation}
Using \eqref{bern}--\eqref{presdiff} and the identity
\begin{align*}
	|u^\pm|^2=\left|\bar u\pm\tfrac{1}{2}w\right|=|\bar u|^2\pm\bar u\cdot w+\tfrac{1}{4}|w|^2 \,,
\end{align*}
we find that along $\Gamma(t)$, 
\begin{subequations} 
\label{diffmean}
\begin{align}
\left(\frac{p^+}{\rho^+}-\frac{p^-}{\rho^-}\right) \circ z &= \left(\partial_t(\varphi^+-\varphi^-) +\tfrac{1}{2}\left(|u^+|^2-|u^-|^2\right)\right) \circ z
	       =\partial_t((\varphi^+-\varphi^-)\circ z) \,, \\
\left(\frac{p^+}{\rho^+}+\frac{p^-}{\rho^-} \right) \circ z &=\left( \partial_t(\varphi^++\varphi^-)+\tfrac{1}{2}\left(|u^+|^2+|u^-|^2\right)\right) \circ z
      +2gz_3 \notag \\
      &=\partial_t((\varphi^+-\varphi^-)\circ z)-|\bar u|^2+\tfrac{1}{4}|w|^2+2gz^3 \,.
\end{align}
\end{subequations} 
Substitution of \eqref{diffmean} into \eqref{presdiff} results in 
\begin{equation}\label{potEuler}
	\partial_t((\varphi^+-\varphi^-)\circ z)+A\partial_t((\varphi^++\varphi^-)\circ z)+A\left(-|\bar u|^2+\tfrac{1}{4}|w|^2+2gz^3\right)=0.
\end{equation}
It is convenient to introduce a ``rotated'' version of the vorticity components,
$$\mu_\alpha=w\cdot\partial_\alpha z\,. \ \ \  (\alpha=1,2)\,, $$
Differentiating \eqref{potEuler} with respect to $s_\alpha$ then yields
\begin{equation}\label{Euler}
	\partial_t(\mu_\alpha+2A\bar u\cdot\partial_\alpha z)=A\partial_\alpha\big(|\bar u|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha\mu_\beta-2gz^3\big) \,,
\end{equation}
where $h^{\alpha\beta}$ is the $(\alpha,\beta)$ component of the inverse of the induced metric on $\Gamma(t)$:
$$h^{\alpha\beta}h_{\beta\gamma}={\delta^\alpha}_\gamma.$$
The coupled system of  equations \eqref{BR2} and \eqref{Euler} are written as
\begin{subequations} 
\label{full-euler}
\begin{align}
	\partial_tz& =\bar u \,,  \label{fe1} \\
	\partial_t(\mu_\alpha+2A\bar u\cdot\partial_\alpha z)& =A\partial_\alpha\big(|\bar u|^2-\tfrac{1}{4}h^{\beta\gamma}\mu_\beta  \mu_\gamma-2gz^3\big) \,, \label{fe2}
\end{align}
\end{subequations} 
and are an equivalent form 
of the 3D irrotational and incompressible  Euler equations \eqref{euler00}, giving the interface position $z=(z^1,z^2,z^3)$ and the components of velocity 
jump $\mu=(\mu_1,\mu_2)$ 
(or equivalently, the vorticity measure).  The system \eqref{full-euler} is supplemented with the initial conditions $z(s,0) =\mathring z(s)$ and $\mu_\alpha(s,0)
= \mathring \mu_ \alpha (s)$.

\subsection{The Need for Approximation}

The numerical solution of the singular integral form of the Euler equations \eqref{full-euler} is difficult because of  the time-derivative term $\partial_t(2A\bar u\cdot\partial_\alpha z)$ in \eqref{fe2}. Expanding this term  using \eqref{fe1} shows that
\begin{align*}
	\partial_t(2A\bar u\cdot\partial_\alpha z)&=\frac{A}{2\pi}\iint_R(\mu_2(s)\partial_1\bar u(s)-\mu_1(s)\partial_2\bar u(s)+\partial_t\mu_2(s)\partial_1z(s)-\partial_t\mu_1(s)\partial_2z(s))\times\frac{z-z(s)}{|z-z(s)|^3}\\
	&-(\mu_2(s)\partial_1z(s)-\mu_1(s)\partial_2z(s))\times(z-z(s))\frac{3(\bar u-\bar u(s))\cdot(z-z(s))}{|z-z(s)|^5}ds \,.
\end{align*}
where we have not written the time-dependence in the integrand.
Hence we see that \eqref{full-euler} takes the form of a system of  nonlinear and nonlocal integral equations for the time-derivatives 
$\partial_t(z^1,z^2,z^2,\mu_1,\mu_2)$. 
In addition to the difficulty of solving \eqref{full-euler} caused by the term $\partial_t(2A\bar u\cdot\partial_\alpha z)$, yet another challenge stems from
fact that the Birkhoff-Rott integral can be time-consuming to evaluate: a naive implementation takes $O(N^2)$ operations for an $N$-point discretization
of $\Gamma(t)$, although this can be  improved greatly by the use of  fast multipole methods.  We are not aware of any successful attempts to directly simulate the singular integral form of the Euler equations \eqref{full-euler}, despite numerous simulations  of these equations in two space dimensions and with axisymmetry.

In two dimensions, Granero-Belinch\'on and Shkoller \cite{GrSh2017} remedied the computational difficulty of the Birkhoff-Rott velocity by introducing an approximate velocity which can be computed as a Hilbert transform, requiring only $O(N\log N)$ evaluations:
\begin{align}
	\widetilde u=\frac{H\omega}{2|\partial_{1}z|}n,&& \widehat{Hf}(k)=-i\,\text{sgn}(k)\hat f(k).  \label{approx-u2d}
\end{align}
The system \eqref{approx-u2d} is the approximation which results from taking a limit of ``small nonlocality'' (explained in detail in the sequel), which is 
accurate when the interface is not too curved; in particular, it was shown in \cite{GrSh2017} that
\begin{equation}
	\max_{s \in \mathbb{R}  }|\bar u(s,t)-\widetilde u(s,t)|\leq\frac{2}{\pi}\sqrt{3K} \int_{\Gamma(t)} \omega(s,t)^2  \sqrt{h} ds\,,
\end{equation}
where $K$ is the maximum curvature of the interface. Inserting the approximate velocity \eqref{approx-u2d} into the 2D version of \eqref{full-euler} provides
a simple set of evolution equations for $z$ and $\omega$, called the 2D \emph{lower-order $z$-model}. It was subsequently verified in \cite{CaDeFrGoRaReSh2020} that the lower-order $z$-model agrees well with experimental data of Rayleigh-Taylor problems, forming a sort of ``envelope'' for the interface rollup.

\subsection{The 3D Lower-Order $z$-Model}

The aim of this section is to introduce a 3D generalization of the  2D lower-order $z$-model of \cite{GrSh2017}, as well as two more accurate models (with greater
nonlocality), which we shall refer to as the ``medium-order'' and ``higher-order'' $z$-models. While already useful for 2D simulations (in respect to the speed-up over
traditional Euler solvers), the need for such a simplification is even greater for 3D simulations.

For the remainder of this section, we omit writing the explicit dependence on time $t$.
The derivation begins with the observation that the dominant contribution to BR fluid velocity $\bar u(s)$ in \eqref{BR1} arises from the
singular integrand in a small neighborhood of $s$.    
Expanding $z(s')$ in a small neighborhood of $s$ yields
\begin{align} 
z(s') = z(s) + \partial_ \alpha z(s) ({s'}^ \alpha - s^ \alpha ) + \tfrac{1}{2} \partial_ {\alpha \beta} z(s) ({s'}^ \alpha - s^ \alpha ) ({s'}^ \beta - s^ \beta )  + O(|s-s'|^3) \,,
\label{temp0}
\end{align} 
and the expansion of $\omega(s')$ about $\omega(s)$ is given by
\begin{align} 
\omega(s') & = \mu_2(s') \partial_1z(s') - \mu_1(s') \partial_2 z(s')  \notag \\
& = \mu_2(s') \partial_1z(s) - \mu_1(s') \partial_2 z(s)+  \mu_2(s')  \partial_{1 \alpha }z(s) ({s'}^\alpha -s^ \alpha )  - \mu_1(s')  \partial_{2 \alpha } z(s)  ({s'}^ \alpha - s ^ \alpha )
+O(|s-s'|^2) \,. \label{temp1}
\end{align}
In order to expand $ \tfrac{1}{|z(s)-z(s')|}$ about the singular point $s=s'$,  we expand $z(s)-z(s')$ about $s=s'$, factor the linear term of this expansion, and use
the series expansion for $ \tfrac{1}{1- \xi} $ about $ \xi$; in particular, we have that
\begin{align}
	\frac{1}{|z(s)-z(s')|^3}&=\frac{1}{|\partial_\alpha z(s)(s'^\alpha-s^\alpha)|^3}\left(\frac{|z(s)-z(s')|}{|\partial_\alpha z(s)(s'-s)|}\right)^{-3} \notag \\
	&=\frac{1}{|\partial_\alpha z(s)(s'^\alpha-s^\alpha)|^3}
	\left(1+\frac{1}{2}\frac{\partial_{\alpha\beta}z(s)(s'^\alpha-s^\alpha)(s'^\beta-s^\beta)}{|\partial_\alpha z(s)(s'^\alpha-s^\alpha)|}+\cdots\right)^{-3}\notag \\
	&=\frac{1}{|\partial_\alpha z(s)(s'^\alpha-s^\alpha)|^3}\left(1-\frac{3}{2}\frac{\partial_{\alpha\beta}z(s)(s'^\alpha-s^\alpha)(s'^\beta-s^\beta)}{|\partial_\alpha z(s)(s'^\alpha-s^\alpha)|}+\cdots\right)
	\label{temp3}
\end{align}
The expansions \eqref{temp0}--\eqref{temp3}, together with a somewhat tedious computation, show that the Birkhoff-Rott integrand 
$\omega(s')\times\tfrac{z(s)-z(s')}{|z(s)-z(s')|^3}$ has the following expansion about $s=s'$:
\begin{align}
&\frac{(s^\alpha-s'^\alpha)}{\big(h_{\alpha\beta}(s)(s^\alpha-s'^\alpha)(s^\beta-s'^\beta)\big)^\frac{3}{2}}
\Big( \mu_\alpha(s') 
+(s^\beta-s'^\beta) \big( \mu_\alpha(s') \partial_\beta\sqrt{h}-  \tfrac{1}{2}\partial_{\alpha\beta}z\cdot\partial_\gamma z)h^{\gamma\delta}(s)\mu_\delta(s') \big)
\Big)  n \notag \\
&
\qquad 
+\frac{(s^\alpha-s'^\alpha)(s^\beta-s'^\beta)}{\big(h_{\alpha\beta}(s)(s^\alpha-s'^\alpha)(s^\beta-s'^\beta)\big)^\frac{3}{2}}\Big(\mu_\alpha(s')
\sqrt h\partial_\beta n
-\tfrac{1}{2}(\mu_2(s')\partial_1z-\mu_1(s')\partial_2z)\times n(\partial_{\alpha\beta}z\cdot n)\Big) \notag \\
&
\qquad
-\tfrac{3}{2}\frac{\mu_\alpha(s')(s^\alpha-s'^\alpha)(s^\beta-s'^\beta)(s^\gamma-s'^\gamma)(s^\delta-s'^\delta)(\partial_{\alpha\beta}z\cdot\partial_\gamma z)\partial_\delta z}{\big(h_{\alpha\beta}(s)(s^\alpha-s'^\alpha)(s^\beta-s'^\beta)\big)^\frac{5}{2}}n+O(1) \,.
	\label{series0}
\end{align}
Isolating the dominant contribution, we thus have that
\begin{align}
\omega(s')\times\frac{z(s)-z(s')}{|z(s)-z(s')|^3}
&=\frac{\mu_1(s')(s^1-{s'}^1)+\mu_2(s')(s^2-{s'}^2)}{\big(h_{\alpha\beta}(s)({s'}^\alpha-s^\alpha)({s'}^\beta-s^\beta)\big)^{3/2}}
\sqrt{h(s)}n(s)+O(|s-s'|^{-1}) \,. \label{series1}
\end{align}
We can then integrate the leading order term  in \eqref{series1}
and obtain an approximate (average) velocity along the interface;  however,  the resulting velocity  field  is difficult to write in terms of simple Fourier multipliers
(which is our objective), and therefore maintains the computational expense of the original Birkhoff-Rott velocity.   In order to derive a velocity field which
can be written in terms of Fourier multipliers and is hence extremely fast to compute, 
we make the further assumption that the metric $h$ is isotropic, by which we mean that the components $h_{ \alpha \beta}$ of the
metric satisfy
 $$ h_{11} = h_{22}\,, \ \ \  h_{12} =h_{21} =0 \,, $$
 in which case
\begin{align} 
\frac{\sqrt{h(s)}\mu_\alpha(s')(s^\alpha-{s'}^\alpha)}{\big(h_{\alpha\beta}(s')({s'}^\alpha-s^\alpha)({s'}^\beta-s^\beta)\big)^{{\frac{3}{2}} }} n(s)
=\frac{\mu_\alpha(s')(s^\alpha-{s'}^\alpha)}{|h(s)|}n(s)\,. \label{temp2}
\end{align} 
Replacing the original integrand in the Birkhoff-Rott integral \eqref{BR2} with the first-order term \eqref{temp2} in its expansion
 results in an approximate velocity $\widetilde u$,  computable in terms of \emph{Riesz transforms} $R^ \alpha $:
\begin{align}\label{utilde}
	\widetilde u=\frac{R^\alpha\mu_\alpha}{2|h|}n,&&R^\alpha f(s)=\frac{1}{2\pi}\iint_{\mathbb{R}^2}\frac{(s^\alpha-s'^\alpha)f(s')}{|s-s'|^3}ds' \,, \ \alpha =1,2 \,.
\end{align}
Riesz transforms are Fourier multipliers: $\widehat{R^\alpha f}(k)=-ik^\alpha\hat f(k)/|k|$, where $\hat f(k)$ denotes the
Fourier transform of $f(x)$.   These operators are the multi-dimensional generalization of the Hilbert transform used in the definition of the lower-order 2D
$z$-model in \eqref{approx-u2d}.



The first term of the Laurent series for 
$\omega(s')\times\tfrac{z(s)-z(s')}{|z(s)-z(s')|^3}$ is given by \eqref{series1} with an error which is $O(|s-s'|^{-1})$ as $s' \to s$.   
The $O(|s-s'|^{-1})$ term in this expansion is what remains in  \eqref{series0}, and upon integration in $s'$,  produces a leading-order estimate  for the error or 
{\it remainder} $\mathcal{R} $ between the lower order velocity $\tilde u$ and the full 
Birkhoff-Rott velocity $\bar u$. Utilizing our isotropy assumption shows that
\begin{align}
	\mathcal R=\frac{\partial_\beta(\sqrt h n)R^\beta R^\alpha\mu_\alpha}{|h|^{3/2}}
	&-\frac{1}{2}\frac{h^{\gamma\delta}(\partial_{\alpha\beta}z\cdot\partial_\gamma z)R^\beta R^\alpha\mu_\delta}{|h|^{3/2}}n \notag\\
	&-\frac{1}{2}\frac{(n\cdot \partial_{\alpha\beta}z)\big((\partial_1z\times n)R^\alpha R^\beta\mu_2-(\partial_2z\times n)R^\alpha R^\beta\mu_1\big)}{|h|^{3/2}}\nonumber\\
	&-\frac{3}{2}\frac{(R^\alpha R^\beta R^\gamma R^\delta\mu_\alpha)(\partial_\beta z\cdot\partial_{\gamma\delta} z)}{|h|^{5/2}}n \,.
	\label{eqn::error}
\end{align}
Using the fact that the Riesz transforms have unit norm as bounded operators on $L^2$, H\"{o}lder's inequality shows that
\begin{equation}
	\label{eqn::errorestimate}
	\left(\iint_{\mathbb{R}^2} \mathcal{R} (s)^2ds\right)^{{\frac{1}{2}} }\leq C\left(\iint_{\mathbb{R}^2}|\mu(s)|_h^2 ds\right)^{\frac{1}{2}}\sum_{\alpha,\beta=1}^2\sup_{s\in \mathbb{T}^2}\frac{|\partial_{\alpha\beta}z|}{1+|h|},
\end{equation}
where $|\mu|_h^2=h_{\alpha\beta}\mu^\alpha\mu^\beta$ is the magnitude of $\mu$ with respect to the metric $h$ on $\Gamma(t)$.
 In more concise notation, 
\begin{align} 
 \|\mathcal{R} \|_{L^2}\leq\|\mu\|_{L^2} \left\| \tfrac{D^2z}{|h|^{\frac{3}{2}}} \right\|_{L^\infty} \,.  \label{crude-bound}
\end{align} 
 Since the mean curvature vector on $\Gamma(t)$ is defined as
 $$
 \kappa(s,t)  n= \tfrac{h^{ \alpha \beta }}{ \sqrt h} \partial_{ \alpha \beta } z \cdot (\partial_1 z \times \partial_2 z) \,,
 $$
 The bound \eqref{crude-bound} becomes small when the mean curvature is very small, as $\|\mu\|_{L^2}$ remains bounded under this scenario.
 


Inserting $\widetilde u$ into the Euler equations \eqref{Euler} in place of the Birkhoff-Rott velocity gets rid of the problematic time-derivative term 
$(\bar u\cdot\partial_\alpha z)_t$, as $\widetilde u$ is normal to $\Gamma(t)$, and so  $\widetilde u \cdot \partial_\alpha z=0$. This yields the following system of equations:
\begin{subequations} 
\label{LO}
\begin{alignat}{2}
	\partial_tz&=\widetilde u \qquad && 0< t\le T \,,\label{LO_Z} \\ 
	\partial_t\mu_\alpha&=A\partial_\alpha\big(|\widetilde u|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha \mu_\beta-2gz^3\big)
	\qquad && 0< t\le T  \,,\label{LO_W} \\
	z(s,0) & =\mathring z(s) \,, \ \ \ \mu_\alpha(s,0) = \mathring \mu_\alpha(s)   \,,
\end{alignat}
\end{subequations} 
where $\mathring z(s)$  and  $ \mathring \mu_\alpha(s)$ denote the initial conditions for the interface parameterization and vorticity, respectively.
We refer to \eqref{LO} as the \emph{lower-order $z$-model}.
\begin{figure}[htb]
	\centering
	\includegraphics[width=5in]{bubble_lo.png}
	\caption{\footnotesize{Evolution of the lower-order $z$-model, starting from a Gaussian perturbation of a flat interface, run on a $192\times192$ mesh in 18s. }}
	\label{bubble_lo}

\end{figure}
 It is an extremely efficient model, which tracks the interface accurately until about the time of turnover, after which it forms a sort of ``envelope'' for the rollup in the near-symmetric regime (see Figure \ref{bubble_lo}). This envelope behavior accurately captures the large-scale behavior of RTI for interfaces with
 small curvature and small anisotropy, 
 although our medium- and higher-order $z$-models (which are introduced below) are better suited for capturing the small-scale structure of rollup.  
 The computational speed of the  lower-order $z$-model creates an efficient and  effective  tool for determining large-scale features of RTI structures, such
 as  interface growth rates and bubble and spike positions.  On average, 
the lower-order $z$-model runs  about 15,000 times  faster than standard numerical schemes for 3D hydrodynamics \cite{Reisner}.

We note that a very nice feature of the lower-order $z$-model \eqref{LO} is that the velocity $\widetilde u$ is given by a Fourier multiplier.  As such, the numerical
implementation of \eqref{LO} can take advantage of the Fast Fourier Transform which is both fast and easy to implement, and is one of advantages of this
model in comparison to the higher-order models to be introduced next.






\subsection{The 3D Medium-Order and Higher-Order $z$-Models}

The isotropy assumption  made in the derivation of the lower-order  $z$-model is too restrictive for a large class of initial data.   As such,
we turn our attention to the medium-order $z$-model \eqref{MO} and higher-order $z$-model \eqref{HO}.  
The medium-order $z$-model  is obtained by replacing the localized velocity \eqref{utilde} with the desingularized Birkhoff-Rott velocity \eqref{ubarepsilon} {\it only}
in  the $z$-evolution equation in \eqref{LO_Z} while keeping the vorticity equation \eqref{LO_W} unchanged:
\begin{subequations} 
\label{MO}
\begin{alignat}{2}
	\partial_tz&=\bar u^\epsilon \qquad && 0< t\le T \,,\label{MO_Z} \\ 
	\partial_t\mu_\alpha&=A\partial_\alpha\big(|\widetilde u|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha \mu_\beta-2gz^3\big)
	\qquad && 0< t\le T  \,,\label{MO_W} \\
	z(s,0) & =\mathring z(s) \,, \ \ \ \mu_\alpha(s,0) = \mathring \mu_\alpha(s)   \,,
\end{alignat}
\end{subequations} 
For the higher-order $z$-model, we replace the localized velocity with the desingularized Birkhoff-Rott velocity in both equations of \eqref{LO}
to obtain
\begin{subequations} 
\label{HO}
\begin{alignat}{2}
	\partial_tz&=\bar u^\epsilon \qquad && 0< t\le T \,,\label{HO_Z} \\ 
	\partial_t\mu_\alpha&=A\partial_\alpha\big(|\bar u|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha \mu_\beta-2gz^3\big)
	\qquad && 0< t\le T  \,,\label{HO_W} \\
	z(s,0) & =\mathring z(s) \,, \ \ \ \mu_\alpha(s,0) = \mathring \mu_\alpha(s)   \,,
\end{alignat}
\end{subequations} 
The medium- and higher-order $z$-models are able to capture the fine-scale structures of RTI, including the Kelvin-Helmholtz roll-up regions (see 
Figure \ref{bubble_hi}), but are more costly to simulate numerically than the lower-order $z$-model. Nevertheless, the medium- and higher-order $z$-models run $600$ times faster than standard numerical methods for 3D gas dynamics codes \cite{Reisner}.

\begin{figure}[h]
	\centering
	\includegraphics[width=5in]{bubble_hi.png}
	\vspace{-.1 in}
	\caption{{\footnotesize Evolution of the higher-order $z$-model, starting from a Gaussian perturbation of a flat interface, run on a $192\times192$ mesh in 420s. Note that the lower-order $z$-model (Figure \ref{bubble_lo}) traces a sort of `envelope' for the higher-order $z$-model.  }}
	\label{bubble_hi}

\end{figure}

The evolution equations for the higher-order $z$-model differ from the integral form of the Euler equations by only one term.  In particular, a comparison of
the right-hand sides of
\eqref{fe2} with  \eqref{HO_W} shows that the higher-order $z$-model does not have the problematic nonlinear and nonlocal
time-derivative term $2A\partial_t(\bar u\cdot\partial_\alpha z)$.   As we will explain in Section \ref{sec::vorticity}, the reason that the higher-order $z$-model
maintains accuracy is that this time-derivative term tends to be small.   We note that for the numerical simulations we consider in this work, results from
the  medium-order $z$-model and the higher-order $z$-model  are somewhat indistinguishable.  For this reason,  in the sequel,  we have chosen to
report results exclusively from the higher-order $z$-model. 

\begin{figure}[h]
	\centering
	\includegraphics[width=2.9in]{f2.pdf}
	\includegraphics[width=2.9in]{f3.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Cross-section of the 3D lower-order $z$-model (\textit{left}) and the 3D higher-order $z$-model (\textit{right}) showing that the lower-order $z$-model fails to correctly model single-mode RTI.}}
	\label{singlemode_comp}
\end{figure}

In contrast to the lower-order $z$-model, the higher-order $z$-model exhibits significant roll-up. More importantly, as shown in Figure \ref{singlemode_comp},  the higher-order $z$-model is able to correctly simulate the ``double rollup'' which distinguishes 3D RTI from 2D RTI.  This important feature of the 3D dynamics will be further discussed in Section \ref{sec::single-mode-data}.




An important feature of Rayleigh-Taylor instability is the differing behavior of `bubbles' of light fluid rising into the heavy fluid and `spikes' of heavy fluid falling into the light fluid. Our model accurately captures this differing behavior with the spikes being thinner and slightly longer than the bubbles. Figure \ref{bubblespike} shows the evolution of the initial conditions
$$z^1=s^1,\qquad z^2=s^2,\qquad z^3=\pm0.05e^{-9|s|^2},\qquad \mu_1=\mu_2=0$$
at two, three, and four characteristic times. 
\begin{figure}[h]
	\centering
	\includegraphics[width=5.5in]{bubble07.pdf}
	\includegraphics[width=5.5in]{spike07.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Cross-sections of the 3D higher-order $z$-model, run on a $64\times 64$ grid, with Atwood number $0.7$ and a runtime of 50s.}}
	\label{bubblespike}
	\vspace{-0.5in}
\end{figure}
\vfill
\clearpage

\subsection{Dynamics of Vorticity}
\label{sec::vorticity}

The relationship between the components of vorticity $\omega^\alpha$ and the interface rollup is quite straightforward. Each `roll' of the interface corresponds to 
a sequence of successively stronger concentric `ridges' in the magnitude of vorticity, and the vector field $(\omega^1,\omega^2)$ circulates around these ridges. 
Figure \ref{vorticity-components} shows a simulation of the higher-order $z$-model \eqref{HO} in which a single ring of vorticity (top row, pre-turnover) splits 
into two rings of vorticity, the inner ring stronger than the outer (bottom row, post-turnover). 
Figure \ref{vorticity-magnitude} shows the corresponding magnitude of vorticity for the interfaces in   Figure \ref{vorticity-components}. These spikes in the vorticity are essential to the rollup of the interface $\Gamma(t)$, but they are also a source of numerical instability. To mitigate this instability, we introduce a  smooth, vorticity-scaled,  nonlinear artificial viscosity, outlined in Section \ref{sec::numerics}.


\begin{figure}[htb]
	\centering
	\includegraphics[width=6in]{vort1.pdf}
	\includegraphics[width=6in]{vort2.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Nondimensionalized components of vorticity for the higher-order $z$-model, compared with the full interface. Color on the interface indicates magnitude of vorticity (red is higher, blue is lower).}}
	\label{vorticity-components}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=6in]{modvort.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Heatmap of the nondimensionalized magnitude of vorticity for the two interfaces shown in Figure \ref{vorticity-components}).}}
	\label{vorticity-magnitude}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=6in]{vort_low.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Nondimensionalized components of vorticity for the lower-order $z$-model, compared with the full interface. Color on the interface indicates magnitude of vorticity (red is higher, blue is lower).}}
	\label{vorticity-lower}
\end{figure}


The lower-order $z$-model fails to roll up because only one {\it sufficiently strong} ring of vorticity can form during the entire evolution. 
A secondary (weaker) ring of vorticity may form inside the first, but its amplitude is too small to initiate secondary turnover. This weak secondary vorticity produced by the lower-order $z$-model is shown in Figure \ref{vorticity-lower}.


We can now provide a heuristic argument as to why the problematic time-derivative term $2A\partial_t(\bar u\cdot\partial_\alpha z)$ on the right side
of \eqref{fe2} does not significantly alter the true dynamics, and thus allows for the higher-order $z$-model evolution to accurately simulate solutions of the Euler equations.
Integrating \eqref{fe2}, we have that
\begin{align} 
\mu_\alpha(s, t)  & =  \mathring\mu_\alpha(s) 
+ A \int_0^t \partial_\alpha\big(|\bar u|^2-\tfrac{1}{4}h^{\beta\gamma}\mu_\beta  \mu_\gamma-2gz^3\big) dt'  \notag \\
&  \qquad \qquad \qquad + 2A \big( \bar u^\epsilon(s,0) \cdot \partial_\alpha \mathring z (s) - \bar u(s,t)\cdot\partial_\alpha z(s,t) \big) \,.  \label{h-arguement1}
\end{align} 
 From  \eqref{series1} and
the definition of $\bar u$ given in \eqref{u-bar}, we see that for $\alpha =1,2$, 
\begin{align*} 
\bar u \cdot \partial_\alpha z
&= \iint_{\mathbb{R}^2  }
\frac{(s^\nu-s'^\nu)(s^\beta-s'^\beta)}{|s-s'|^3}\Big(\mu_\nu(s')
\sqrt h\partial_\beta n
-\tfrac{1}{2}(\mu_2(s')\partial_1z-\mu_1(s')\partial_2z)\times n(\partial_{\nu\beta}z\cdot n)\Big) \cdot  \partial_\alpha z  ds' \\
&
=-\iint_{\mathbb{R}^2  }\frac{(s^\nu-s'^\nu)(s^\beta-s'^\beta)}{|s-s'|^3}\Big(\mu_\nu(s')
h^{\delta \nu}\partial_\delta z
+\tfrac{1}{2}(\mu_2(s')\partial_1z-\mu_1(s')\partial_2z)\times n\Big) \cdot  \partial_\alpha z \ (\partial_{\nu\beta}z\cdot n)ds',
\end{align*} 
where we have used the identity
$$
\partial_\beta n = - h^{\delta \nu} ( \partial_{\nu \beta}z \cdot n) \partial_\delta z \,,
$$
and the fact that
$$
|s-s'|^3= \big(h_{\nu\beta}(s)(s^\nu-s'^\nu)(s^\beta-s'^\beta)\big)^\frac{3}{2} \,.
$$
Since the inner-product of the first-order term in \eqref{series1} with $\partial_\alpha z$ vanishes, 
$$
\frac{(s^\nu-{s'}^\nu)}{ |s-s'|^3} \mu_\nu(s')\sqrt{h(s)}n(s)\cdot \partial_\alpha z =0\,,
$$
the second-order correction is thus given by
\begin{align} 
\frac{(s^\nu-s'^\nu)}{|s-s'|^3} (s^\beta-s'^\beta)\Big(\mu_\nu(s')
h^{\delta \nu}\partial_\delta z
+\tfrac{1}{2}(\mu_2(s')\partial_1z-\mu_1(s')\partial_2z)\times n\Big) \cdot  \partial_\alpha z \ (\partial_{\nu\beta}z\cdot n) \,.   \label{h-arguement2}
\end{align} 
Now, either the curvature of $\Gamma(t)$ is small in which case $\partial_{\nu\beta}z\cdot n$ is small and hence $\bar u \cdot \partial_\alpha z$ is small, or the
curvature is large, in which case  the vorticity measure forms highly concentrated peaks (as seen in Figure \ref{vorticity-components}), in which the support
of the peaks, which is a good approximation for $|s^\beta-s'^\beta|$, is very small.   Hence, we see from \eqref{h-arguement2} that either small curvature
or strongly concentrated vorticity peaks produces small tangential velocity components.




\subsection{Linear Analysis}
The kernel desingularization \eqref{kernel-desing} is responsible for frequency damping.  In order to explain this, 
we shall linearize our $z$-models  about the flat interface $z=(s^1,s^2,0)$ with zero vortex strength $\mu_1=\mu_2=0$.  In the absence of desingularization, all three
$z$-models  (as well as the Euler equations) have the same linearization:
\begin{align*}
	\partial_t z & =-\tfrac{1}{2}R^\alpha \mu_\alpha\hat e_3\,,  \\
	 \partial_t\mu&=2Ag\nabla_sz^3 \,, 
\end{align*}
where $\nabla_s=(\partial_1,\partial_2)$. Differentiating the equation for $z^3$ once again yields
$$\partial_t^2z_{3}=-AgR^\alpha\partial_\alpha z^3=Ag(-\Delta)^\frac{1}{2}z^3.$$
or equivalently $\hat z^3_{tt}(k)=Ag|k|\hat z^3(k)$ in Fourier space. Provided we are in the Rayleigh-Taylor regime (i.e. $A>0$), this implies that the mode with wavenumber $k$ has amplitude
\begin{align*}
	a_k(t)\sim e^{\sigma(k)t},&&\sigma(k)=(Ag|k|)^\frac{1}{2} \,.
\end{align*}
This growth rate is characteristic of the Rayleigh-Taylor instability, and is what makes the RTI initial value problem linearly ill-posed. When the desingularization is added, our linearization has the form
\begin{align*}
	\partial_t z & =\frac{\hat e_3}{4\pi}\iint_{\mathbb{R}^2} w(s')\cdot\nabla_s\frac{1}{\sqrt{\epsilon^2+|s-s'|^2}}ds' \,, \\
	 \partial_t \mu&=2Ag\nabla_sz^3 \,.
\end{align*}
Making use of the Fourier transform identity
$\mathscr F\left\{\tfrac{1}{\sqrt{\epsilon^2+|s|^2}}\right\}=\tfrac{e^{-\epsilon|k|}}{|k|}$,
we find that $\partial_t^2\hat z^3(k)=-Ag|k|e^{-\epsilon|k|}\hat z_3(k)$. This yields a growth rate which is exponentially damped as $|k|\to\infty$:
$$\sigma^\epsilon(k)=\big[Ag|k|e^{-\epsilon|k|}\big]^\frac{1}{2}=e^{-\epsilon|k|/2}\sigma(k).$$
Hence the desingularization damps the growth of the highest-frequency modes, so that the fastest-growing mode has frequency $k=\tfrac{2}{\epsilon}$. 

\subsection{Numerical Approximation of the $z$-Model}
\label{sec::numerics}

All computations are set on a uniform spatial grid $\{(i\delta,j\delta):i,j=-n,\dots,n-1,n\}$, and the dynamical variables are stored as matrices of values on this grid, e.g.,
\begin{align*}
	z_{ij}(t)=z(i\delta,j\delta,t),&&i,j=-n,\dots,n-1,n.
\end{align*}
Spatial derivatives are computed using fourth-order centered-difference operators $(D_1,D_2)$. These are used to compute the discretized form of the normal vector $N=(D_1z)\times(D_2z)$ and the metric determinant $|h|=|N|^2$.

For the lower-order $z$-model,  we compute $\tilde u$ via the fast Fourier transform:
$$\tilde u_{jk}=\mathscr F^{-1}\left(\frac{il(\mathscr F\mu_1)_{lm}+im(\mathscr F\mu_2)_{lm}}{(l^2+m^2)^{1/2}}\right)_{jk}\frac{N_{jk}}{|h|_{jk}^{3/2}} \,. $$
This still leaves $(\mathscr F\tilde u)_{00}$ undefined, so we set it to zero. To mitigate the large spikes in vorticity discussed in Section \ref{sec::vorticity}, 
we introduce a nonlinear artificial viscosity operator $\nu \nabla\cdot\big(\frac{c\nabla}{\sup_s c}\big)$  on the right side of  the vorticity equations, where
$$c=(1-\nu\delta^2\Delta)^{-1}|\omega| \,, $$
where $\delta$ is the grid size and $\nu >0$ denote the artificial viscosity parameter.
The function $c$ is a  smoothed version of the magnitude of vorticity $|\omega|$, which provides a smooth localization for the addition of nonlinear artificial viscosity,  
 motivated by the nonlinear (space-time smooth) artificial viscosity method introduced in \cite{ReSeSh2013, RaReSh2019a, RaReSh2019b} but 
employing an elliptic solver at each time-step rather than the solution of a  parabolic reaction-diffusion equation. 
This yields the following set of ODEs for the discretized lower-order $z$-model:
\begin{align*}
	\frac{dz_{jk}}{dt}&=\widetilde u_{jk} \,, \\
	\frac{d\mu_{jk}}{dt}&=A\left[\binom{D_1}{D_2}\left(|\widetilde u|^2-\tfrac{1}{4} h^{\alpha\beta}\mu_\alpha\mu_\beta-2gz^3\right)\right]_{jk}+\nu\left[\sum_{\alpha=1}^2D_\alpha\left(\frac{cD_\alpha\mu}{\sup c}\right)\right]_{jk} \,,
\end{align*}
which are solved using a total variation decreasing (TVD) third-order Runge-Kutta method in time. 

For  both the medium and higher-order $z$-models, we compute the Birkhoff-Rott velocity $\bar u$ using a Simpson  $\tfrac{3}{8}$-rule quadrature, which has 
the form
\begin{equation*}
	\bar u^\epsilon_{ij}=\tfrac{1}{4\pi}\sum_k\sum_lM_{kl}\omega_{kl}\times\frac{z_{ij}-z_{kl}}{(\epsilon^2+|z_{ij}-z_{kl}|^2)^{3/2}},\qquad\omega_{kl}=\big(\mu_{2,kl}(D_1z)_{kl}-\mu_{1,kl}(D_2z)_{kl}\big) \,,
\end{equation*}
for some matrix of weights $(M_{ij})$. The discretized medium-order $z$-model is given by
\begin{align*}
	\frac{dz_{jk}}{dt}&=\bar u^\epsilon_{jk} \,, \\
	\frac{dw_{jk}}{dt}&=A\left[\binom{D_1}{D_2}\left(|\widetilde u|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha\mu_\beta-2gz^3\right)\right]_{jk}+\nu\left[\sum_{\alpha=1}^2D_\alpha\left(\frac{cD_\alpha\mu}{\sup c}\right)\right]_{jk} \,,
\end{align*}
and the discretized higher-order $z$-model is
\begin{align*}
	\frac{dz_{jk}}{dt}&=\bar u^\epsilon_{jk} \,, \\
	\frac{d\mu_{jk}}{dt}&=A\left[\binom{D_1}{D_2}\left(|\bar u^\epsilon|^2-\tfrac{1}{4}h^{\alpha\beta}\mu_\alpha\mu_\beta-2gz^3\right)\right]_{jk}+\nu\left[\sum_{\alpha=1}^2D_\alpha\left(\frac{cD_\alpha\mu}{\sup c}\right)\right]_{jk} \,.
\end{align*}
\begin{figure}
	\centering
	\includegraphics[width=3.5in]{benchmark.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Measured computation times for the Birkhoff-Rott velocity.}}
	\label{benchmark}
\end{figure}

The  time required for the computation of the Birkhoff-Rott velocity $\bar u$ can be drastically shortened with the use of the fast multipole method. 
Evaluating $\bar u^\epsilon_{jk}$ for all $j$ and $k$ requires $O(N^2)$ operations, where $N=(2n+1)^2$ is the total number of grid points in the simulation. 
In contrast, the fast multipole method evaluates in $O(N)$ operations, given a fixed relative accuracy goal.  As shown in Figure \ref{benchmark}, we have 
verified this performance for interface grid resolutions up to size $N=1024\times1024$ using the open-source package 
\texttt{fmm3d}.\footnote{https://github.com/flatironinstitute/FMM3D}



\section{Simulation of classical RTI Test Problems}

In this section we use our higher-order $z$-models to study a variety of classical test problems for the Rayleigh-Taylor instability.




\subsection{Single-Mode Initial Data}
\label{sec::single-mode-data}

Sinusoidal perturbations of a flat interface were the first to be studied in the context of RTI, dating back to the work of Lord Rayleigh in 1883, although his work 
was only in the linear regime. Since then, much of the analytical and experimental study of RTI has focused on understanding the evolution of single-mode initial 
data beyond the linear regime (see, for instance, \cite{ChEmWa1960}\cite{CoTa1973}\cite{BaMeOr1984}\cite{Go2002}). 
Good numerical simulations of RTI required the development of special-purpose numerical methods, since the rectangular cells of Eulerian differencing schemes 
are ill-suited to the highly curved geometry of RTI problems, and often introduce spurious KHI spirals which obscure the form of the primary instability. 
To improve accuracy, methods such as marker-and-cell \cite{Da1967} and volume-of-fluid \cite{HiNi1981} use either marker particles or a marker function in addition to an Eulerian grid to track the interface. Similarly, front-tracking methods, which use a triangular mesh that advects with the fluid velocity, have shown considerable success \cite{TrUn1991}. Of course, in most application scenarios, the interface is perturbed with many different frequencies, but to understand the single-mode problem is the first step towards understanding the more general problem.


\begin{figure}[htb]
	\centering
	\includegraphics[width=5in]{k=2.pdf}
	\vspace{-.15 in}
	\caption{{\footnotesize Single-mode initial data with two wavelengths along each side of the tank, run using the higher-order $z$-model on a $192\times192$ grid, with Atwood number 0.15 and a computation time of 380s. }}
	\label{onemode}

\end{figure}

For our first test, we use our higher-order $z$-model to simulate Wilkinson/Jacobs' isopropyl alcohol/water experiment ($A=0.15$). We initialize our $z$-model with initial conditions of the form
$$z^1=s^1,\quad z^2=s^2,\quad z^3=a_0\cos(\pi k s^1/\lambda)\cos(\pi k s^2/\lambda),\quad \mu_1=\mu_2=0,$$
where $-\lambda\leq s^1,s^2\leq\lambda$ and $k$ is a positive integer.   We may regard $\lambda$ as a characteristic length for this problem, and the corresponding characteristic time scale is given by $\tau=(\lambda/Ag)^{1/2}$.  This type of initial data was studied experimentally by Wilkinson and Jacobs,
 who used planar laser-induced fluorescence (PLIF) to visualize a diagonal cross-section of the fluid as it evolved, showing the 
{\it double rollup} which distinguishes 3D Rayleigh-Taylor instability from 2D \cite{JaWi2007}.
Results for the case $k=2$ are shown in Figure \ref{onemode}.


Starting from a single-mode perturbation of a flat interface, 3D RTI forms a periodic pattern of bubbles and spikes in such a way that \emph{two} vortices form on the side of each bubble or spike. In 2D, single-mode RTI produces a stationary row of counter-rotating KH vortices, which manifests as a single row of rollup regions. However, in 3D, single-mode RTI produces two counter-rotating grids of vortex \emph{rings}, which move apart from each other at roughly constant speed. 
Upon taking a cross-section along the diagonal of the tank, we obtain the double-spiral pattern seen here. 

\begin{figure}[h]
	\centering
	\includegraphics[width=4.5in]{comp160_high.pdf}
	\includegraphics[width=4.5in]{comp160_med.pdf}
	\vspace{-.2 in}
	\caption{ {\footnotesize \textit{Top:} Higher-order $z$-model. \textit{Bottom:} Medium-order $z$-model. Both were run on $160\times160$ meshes, with Atwood number 0.15 and a runtime of 270s. Experimental photographs from Wilkinson and Jacobs \cite{JaWi2007}. Copyright 2007, AIP Publishing. Reproduced with permission. }}
	\label{onemode_comparison}
\end{figure}
\begin{figure}[h]
	\centering
	\vspace{-5mm}
	\includegraphics[width=4.5in]{comp_low.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Lower-order $z$-model, run on a $130\times130$ grid, with Atwood number $0.15$ and a runtime of $9$s. Experimental photographs from Wilkinson and Jacobs \cite{JaWi2007}. Copyright 2007, AIP Publishing. Reproduced with permission.}}
	\label{onemode_lo}

\end{figure}


Using the parameters from their experiments and matching the wavelength of the initial data as best we could, we were able to achieve remarkably good agreement between our higher- and medium-order $z$-models and their experimental data -- see Figure \ref{onemode_comparison}. As shown in Figure \ref{onemode_lo}, the lower-order $z$-model matches the experiments quite well at early times.  At later times,  however, the lower-order
$z$-model cannot capture the full extent of the rollup, as discussed in Section \ref{sec::vorticity}.






\subsection{The 3D Rocket Rig Experiment}
\label{sec::rocketrig}

The growth rate of the mixing layer between two fluids subjected to random multi-mode perturbations has been the subject of significant interest and study, and it forms a useful benchmark for our $z$-models. Starting with a small perturbation of amplitude $h_0$ of a flat interface, theory predicts that the amplitude $h(t)$ of the resulting mixing layer has the form
$$
\frac{h-h_0}{\lambda}=\alpha\left(\frac{t}{\tau}\right)^2 \,,
$$
where $\alpha$ is a universal growth rate constant. The growth rate of the interface was tested experimentally by Read \cite{Re1984} and 
Youngs \cite{Yo1989}, and numerically by Youngs \cite{Yo1984}\cite{Yo1989}.  Their experimental setup consists of a $150\times150$mm tank on a vertical sled, accelerated downward using rocket motors, which gives this experimental setup the nickname ``rocket rig''.   
Read's results suggest that $\alpha$ is between $0.6$ and $0.7$, although individual experiments have given values as low as $0.58$ and as high as $0.77$.  
To test the higher-order $z$-model in this scenario, we simulate Read's experiment using sodium iodide (NaI) solution $(\rho^+=1.9\text{ g/cm}^3)$ and pentane $(\rho^-=0.63\text{ g/cm}^3)$ , which results in an Atwood number $A\approx 0.5$. We initialize our model a random perturbation of a flat interface, of the form
$$z_1=s_1,\qquad z_2=s_2,\qquad z_3=A_0\text{Re}\sum_{j_1=-n}^n\sum_{j_2=-n}^na_{j_1j_2}e^{i(j_1s_1+j_2s_2)}$$
where the complex coefficients $a_{j_1j_2}$ are drawn from a standard normal distribution, and $A_0$ is chosen so that $(\iint|z_3(s,0)|^2ds)^{1/2}=0.05$. 
\vspace{-0.2in}
\begin{figure}[htb]
	\centering
	\includegraphics[width=3.5in]{rocketrig_growth.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize Comparison of higher-order $z$-model with theoretical growth rate. The orange dots show the average of 100 different runs, each with random initial data, and the gray bands show one and two standard deviations. }}
	\label{growth}
	\vspace{-0.2in}
\end{figure}

Figure \ref{growth} shows results from 100 simulations of the higher-order $z$-model replicating Read's NaI/Pentane experiment, with an acceleration of around 25 g's. Both the experiment and the simulation were run for 67.4ms, or $2.7\tau$.  
These results show excellent agreement with the experimental growth rate $\alpha=0.06$. The computation time was 100s/run, using $100\times100$ meshes. 
To our knowledge, this is the first interface model which can accurately reproduce the growth rate of disordered RT mixing layers.



\begin{figure}[h]
	\centering
	\includegraphics[width=3in]{growth_lo1.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize  Comparison of lower-order $z$-model with theoretical growth rate. }}
	\label{growth_lo}
\end{figure}

The lower-order model does quite well also,  but the rate of interface growth slows after about two characteristic times, because the lower order $z$-model's is not able to fully capture interface roll up. This is shown in Figure \ref{growth_lo}. 


\subsection{Fluid mixing in the Rocket Rig experiment}
\label{sec::mixing}
In addition to measuring the growth rate of the mixing layer, the speed of our high-order $z$-model simulations makes it feasible to compute ensemble averages of  density for randomized initial data as shown in Figure \ref{rocketrig_mixing} for the case of an ensemble of 100 runs.   When computing ensembles using the invariant (Gibbs) measure associated to the Euler equations, it is known (at least in two space dimensions \cite{AlCr1990}) that ensemble averaging of runs with data using (Gaussian) random 
frequency fluctuations of the interface produce viscous effects (whose size is inversely proportional to the frequency of perturbation).    Such viscous effects can be computed
from a uniform ensemble\footnote{For the purpose of numerical simulations is much easier to compute the uniform ensemble than the ensemble averaging using the invariant
measure, but the viscous effects are extremely similar.} of runs, and our results for such ensembles is in good agreement with traditional simulations of RTI using
the Navier-Stokes equations.



\begin{figure}[h]
	\centering
	\includegraphics[width=6.5in]{rr_mixing.pdf}
	\vspace{-.2 in}
	\caption{{\footnotesize  Ensemble averaged density for the 100 runs of the rocket rig analyzed above, showing a cross-section at $y=0$.  }}
	\label{rocketrig_mixing}
\end{figure}

As explained in \cite{Th2012}, the mixing of a passive scalar (in our case the density $\rho$),  advected by a velocity field, can be measured quantitatively using the homogeneous Sobolev norms
$$\|\rho\|_{\dot H^s}=\left(\int_{\mathbb{R}^3}|k|^{2s}|\hat\rho(k)|^2dk\right)^{1/2}$$
with $s<0$. Of particular interest, are the so-called {\it mixing norms} corresponding to the exponents $s=-1$ and $s=-\tfrac{1}{2}$.  Letting
$\bar \rho$ denote the spatial average of $\rho$, a fluid  becomes mixed when varying density field $\rho(x,t)$ converges to its average $\bar \rho$.   Optimal mixing occurs when $ \rho(x,t)$ converges  exponentially fast to its average in one of these mixing norms; that is, optimal mixing occurs if
$$\|\rho( \cdot , t) - \bar \rho \|_{\dot H^s}\sim e^{-rt}\,, \ \ \ s=-1 \text{ or } -\tfrac{1}{2} $$
for some constant $r>0$ as time  $t\to\infty$ (see, for example,  \cite{Th2012}). 
By computing the $\dot H^{-\frac{1}{2}}$ and $\dot H^{-1}$ norms of the ensemble-averaged rocket rig density in a small strip around the plane $z=0$, we find that the Sobolev norms obey the optimal decay rate (Figure \ref{mixingnorms}).   In other words, the ensemble average of 100 $z$-model simulations of interfaces that are randomly perturbed produces optimal mixing of the two
fluids.
\begin{figure}
	\centering
	\includegraphics[width=3in]{mixingnorms4.pdf}
	\caption{{\footnotesize  Mixing norms of the ensemble-averaged density shown in Figure \ref{rocketrig_mixing}, shown on a log scale. }}
	\label{mixingnorms}
\end{figure}

\subsection{The 3D Tilted Rig}
\label{sec::tiltedrig}

Similar to the standard rocket rig, the tilted rocket rig consists of the small tank mounted on the rocket rig, but this time the tank is tilted at around $5.7^\circ$ from vertical. This causes the instability to grow more rapidly at the edges of the tank and increase the slope of the fluid interface. The purpose of the this experiment 
was to excite the instability in such a strong fashion as to make the fluids interchange position, so that the heavy fluid ends up on the bottom of the tank
and the light fluid 
on top. For our simulation, the higher-order $z$-model was initialized with a random perturbation off the tilted interface,
$$z^1=s^1,\qquad z^2=s^2,\qquad z^3=\tfrac{s^1}{10}+\text{Re}\sum_{k_1,k_2=-n}^nc_{k_1,k_2}e^{2\pi i(k_1s^1+k_2s^2)}.$$ 
Our parameters are matched to Youngs' NaI/Hexane experiment; that is, $A=0.5$, $g=15\times 9.81\text{ m/s}^2$. We use a grid resolution of $160\times160$.

No analytical predictions that we could test have been made about this scenario, but  as shown in Figure \ref{tiltedrig_comparison}, our models show excellent qualitative agreement with the experimental and numerical results from \cite{Yo1989}.
\begin{figure}[htb]
	\centering
	\includegraphics[width=6in]{tr.png}
	\vspace{-.2 in}
	\caption{{\footnotesize Side view of the material interface computed using the higher-order $z$-model, compared with  high-speed photography  from Youngs \cite{Yo1989}. Copyright 1989, Elsevier. Reproduced with permission.} }
	\label{tiltedrig_comparison}
\end{figure}




\subsection{Rising Bubble}
\label{sec::bubble}
The examples treated thus far have consisted of horizontally-periodic interface motion, however, the $z$-model can be employed for simulations in which
the interface is a closed two-dimensional surface.   An important test problem requires modeling the motion of an initially-spherical bubble, rising in an
inviscid fluid. 

The motion of an initially spherical gas bubble in water was one of the earliest to be studied in the literature, as it can be treated as a perturbation problem 
 of potential flow around a sphere. However,  outside of such a perturbative regime, the problem becomes very difficult, both analytically and numerically. Grid-based simulations of rising bubbles suffer from the same spurious KH instability as Rayleigh-Taylor simulations, so even early simulations used adaptive, boundary-fitted grids \cite{LeRy1984}. Front-tracking methods where the rising bubble's surface is modeled by a triangular mesh have met with considerable success (see, for example, \cite{TrBuEsJuAl2001}), but are complicated to implement and have poorly-understood convergence properties. 
\begin{figure}[h]
	\centering
	\includegraphics[height=1.45in]{sphere_exp1.png}
	\includegraphics[height=1.45in]{sphere_exp2.png}
	\includegraphics[height=1.45in]{sphere_exp3.png}
	\includegraphics[height=1.5in]{sphere1.png}
	\includegraphics[height=1.5in]{sphere2.png}
	\includegraphics[height=1.5in]{sphere3.png}
	\vspace{-.2in}
	\caption{{\footnotesize \emph{Top:} High-speed photography from the experiments of Bhaga and Weber \cite{BhWe1981}, showing the oblate spherical cap, oblate ellipsoidal cap, and skirted regimes respectively. Copyright 1981, Cambridge University Press. Reproduced with permission. \emph{Bottom:} Snapshots of the higher-order $z$-model, showing the same three regimes, run on a $160\times80$ grid, with a computation time of 150s.}}
	\label{bubble_comparison}
\end{figure}

We use the higher-order $z$-model to simulate rising-bubble motion, where the initial interface parameterization of the bubble given by
$$z^1=0.5\cos(\pi s^1)\cos( \tfrac{\pi s^2}{2}),\qquad z^2=0.5\sin(\pi s^1)\cos( \tfrac{\pi s^2}{2}),\qquad z^3=0.5\sin( \tfrac{\pi s^2}{2}),$$
with density inside the bubble given by $\rho_{\operatorname{in} } =0.0012\text{ g/cm}^3$ and density outside the bubble given by $\rho_{\operatorname{out} }=1.3\text{ g/cm}^3$ and
hence an Atwood number of $A=0.99$.  To avoid certain topological entanglement at the poles of the sphere, we excise very small disks centered at both poles.
Our mesh has 160x80 cells.


Given that $\rho_{\operatorname{in} }  \le \rho_{\operatorname{out} } $, the bubble rises and both the Rayleigh-Taylor and Kelvin-Helmholtz instabilities are
initiated.
As shown in Figure \ref{bubble_comparison}, the $z$-model simulations show many of the characteristic bubble shapes observed in the
high-speed photography from the experiments of Bhaga and Weber \cite{BhWe1981}.
For large enough bubbles, inertial effects dominate surface tension and viscosity, and our inviscid model  appears to capture the basic roll-up structure.
In this situation, an initially spherical bubble goes through several characteristic stages before a topological transition occurs, at which point our model no longer applies. Figure \ref{bubble_comparison} shows a comparison between the higher-order $z$-model and experiments for three of these characteristic regimes.

\section{Conclusions and Future Work}



We have derived a straightforward and efficient model for interface evolution in 3D RTI, and we have demonstrated its efficacy for simulating single-mode and random multi-mode RTI. It is robust enough to simulate spherical interface geometries, and fast enough to study viscous mixing through ensemble averages of randomly perturbed inviscid simulations. 
To our knowledge, this is the only three-dimensional interface model which has the ability to reproduce experimental results with the accuracy we have obtained. 
Computation times for our models are already much better than for equivalent three-dimensional Eulerian finite difference schemes, but with the addition of a fast multipole method, computations can be further sped up by orders of magnitude. All of our code is written in MATLAB, without any concern for speed, and is run on a 2.3GHz 8-core Intel i9 processor. Thus we expect that much better optimization is possible. 

As has been done in two dimensions by Ramani and Shkoller \cite{RaSh2020}, we intend to include the $z$-model in a three-dimensional multiscale compressible flow code, where the incompressible part of the flow is governed by the $z$-model. Other work will study the analytical properties of the medium- and higher-order $z$-models, investigate their asymptotic behaviors in certain well-understood regimes, perform more comprehensive comparisons to experimental data, and extend the interface model to a level-set method which can allow for topological transitions.

\subsection*{Acknowledgments} 
We would like to thank Jon Reisner at Los Alamos for many fruitful discussions on the subject of this paper.
Research reported in this publication was supported by NSF grant DMS-2007606 and   DTRA grant HDTRA11810022.  This research was also supported by
 Office of Defense Nuclear Nonproliferation Research and Development NA-24 and we note that the views of the authors do not necessarily reflect the views of the USG.  This work was also supported by the Laboratory Directed Research and Development Program of the Los Alamos National Laboratory, which
is under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy under
DOE Contracts W-7405-ENG-36 and LA-UR-10-04291.   This research was also supported by Office of Defense Nuclear Nonproliferation Research and Development NA-22.
 


\begin{bibdiv}
\begin{biblist}

\bib{AlCr1990}{article}{
      author={Albeverio, S.},
      author={Cruzeiro, A.-B.},
       title={Global Flows with Invariant (Gibbs) Measures for Euler and Navier-Stokes Two Dimensional Fluids},
        date={1990},
     journal={Commun. Math. Phys},
      volume={129},
       pages={431\ndash 444},
}

\bib{BaBe2004}{article}{
      author={Baker, G.~R.},
      author={Beale, J.~T.},
       title={Vortex blob methods applied to interfacial motion},
        date={2004},
     journal={Journal of Computational Physics},
      volume={196},
      number={1},
       pages={233\ndash 258},
}

\bib{BaMeOr1984}{article}{
      author={Baker, G.~R.},
      author={Meiron, D.~I.},
      author={Orszag, S.~A.},
       title={Boundary integral methods for axisymmetric and three-dimensional
  rayleigh-taylor instability problems},
        date={1984},
     journal={Physica D: Nonlinear Phenomena},
      volume={12},
      number={1-3},
       pages={19\ndash 31},
}

\bib{BaPh2006}{article}{
      author={Baker, G.~R.},
      author={Pham, L.~D.},
       title={A comparison of blob methods for vortex sheet roll-up},
        date={2006},
     journal={Journal of Fluid Mechanics},
      volume={547},
       pages={297\ndash 316},
}

\bib{BeMa1982}{article}{
      author={Beale, J.~T.},
      author={Majda, A.},
       title={Vortex methods. i. convergence in three dimensions},
        date={1982},
     journal={Mathematics of computation},
      volume={39},
      number={159},
       pages={1\ndash 27},
}

\bib{BeMa1985}{article}{
      author={Beale, J.~T.},
      author={Majda, A.},
       title={High order accurate vortex methods with explicit velocity
  kernels},
        date={1985},
     journal={Journal of Computational Physics},
      volume={58},
      number={2},
       pages={188\ndash 208},
}

\bib{BhWe1981}{article}{
      author={Bhaga, D.},
      author={Weber, M.~E.},
       title={Bubbles in viscous liquids: shapes, wakes and velocities},
        date={1981},
     journal={Journal of fluid Mechanics},
      volume={105},
       pages={61\ndash 85},
}

\bib{Bi1962}{inproceedings}{
      author={Birkhoff, G.},
       title={Helmholtz and taylor instability},
        date={1962},
   booktitle={Proc. symp. appl. math},
      volume={13},
       pages={55\ndash 76},
}

\bib{BiFi1959}{article}{
      author={Birkhoff, G.},
      author={Fisher, J.},
       title={Do vortex sheets roll up?},
        date={1959},
     journal={Rendiconti del Circolo matematico di Palermo},
      volume={8},
      number={1},
       pages={77\ndash 90},
}

\bib{CaDeFrGoRaReSh2020}{article}{
      author={Canfield, J.},
      author={Denissen, N.},
      author={Francois, M.},
      author={Gore, R.},
      author={Rauenzahn, R.},
      author={Reisner, J.},
      author={Shkoller, S.},
       title={A comparison of interface growth models applied to
  Rayleigh--Taylor and Richtmyer--Meshkov instabilities},
        date={2020},
     journal={Journal of Fluids Engineering},
      volume={142},
      number={12},
       pages={121108},
}

\bib{Ch1973}{article}{
  title={Numerical study of slightly viscous flow},
  author={Chorin, A.~J.},
  journal={Journal of fluid mechanics},
  volume={57},
  number={4},
  pages={785--796},
  year={1973},
  publisher={Cambridge University Press}
}

\bib{CoTa1973}{article}{
      author={Cole, R.~L.},
      author={Tankin, R.~S.},
       title={Experimental study of Taylor instability},
        date={1973},
     journal={The Physics of Fluids},
      volume={16},
      number={11},
       pages={1810\ndash 1815},
}

\bib{Da1967}{article}{
      author={Daly, B.~J.},
       title={Numerical study of two fluid rayleigh-taylor instability},
        date={1967},
     journal={The Physics of Fluids},
      volume={10},
      number={2},
       pages={297\ndash 307},
}

\bib{ChEmWa1960}{article}{
      author={Emmons, H.~W.},
      author={Chang, C.~T.},
      author={Watson, B.~C.},
       title={Taylor instability of finite surface waves},
        date={1960},
     journal={Journal of Fluid Mechanics},
      volume={7},
      number={2},
       pages={177\ndash 193},
}

\bib{Go2002}{article}{
      author={Goncharov, V.~N.},
       title={Analytical model of nonlinear, single-mode, classical
  rayleigh-taylor instability at arbitrary atwood numbers},
        date={2002},
     journal={Physical review letters},
      volume={88},
      number={13},
       pages={134502},
}

\bib{GrSh2017}{article}{
      author={Granero-Belinch{\'{o}}n, R.},
      author={Shkoller, S.},
       title={A model for rayleigh--taylor mixing and interface turnover},
        date={2017-01},
     journal={Multiscale Modeling {\&} Simulation},
      volume={15},
      number={1},
       pages={274\ndash 308},
         url={https://doi.org/10.1137/16m1083463},
}

\bib{Ha1978}{article}{
  title={Convergence of vortex methods for Euler’s equations},
  author={Hald, O. and Del Prete, V.~M.},
  journal={Mathematics of Computation},
  volume={32},
  number={143},
  pages={791--809},
  year={1978}
}

\bib{Ha1979}{article}{
  title={Convergence of vortex methods for Euler’s equations. II},
  author={Hald, O.~H.},
  journal={SIAM Journal on Numerical Analysis},
  volume={16},
  number={5},
  pages={726--755},
  year={1979},
  publisher={SIAM}
}

\bib{Ha1987}{article}{
  title={Convergence of vortex methods for Euler’s equations, III},
  author={Hald, Ole H},
  journal={SIAM journal on numerical analysis},
  volume={24},
  number={3},
  pages={538--582},
  year={1987},
  publisher={SIAM}
}

\bib{HiNi1981}{article}{
      author={Hirt, C.~W.},
      author={Nichols, B.~D.},
       title={Volume of fluid (vof) method for the dynamics of free
  boundaries},
        date={1981},
     journal={Journal of computational physics},
      volume={39},
      number={1},
       pages={201\ndash 225},
}

\bib{Kr1986b}{article}{
      author={Krasny, R.},
       title={Desingularization of periodic vortex sheet roll-up},
        date={1986},
     journal={Journal of Computational Physics},
      volume={65},
      number={2},
       pages={292\ndash 313},
}

\bib{Kr1986}{article}{
      author={Krasny, R.},
       title={A study of singularity formation in a vortex sheet by the
  point-vortex approximation},
        date={1986},
     journal={Journal of Fluid Mechanics},
      volume={167},
       pages={65\ndash 93},
}

\bib{LiWe2003}{article}{
      author={Liska, R.},
      author={Wendroff, B.},
       title={Comparison of several difference schemes on 1d and 2d test
  problems for the euler equations},
        date={2003},
     journal={SIAM Journal on Scientific Computing},
      volume={25},
      number={3},
       pages={995\ndash 1017},
}

\bib{RaReSh2019a}{article}{
      author={Ramani, R.},
      author={Reisner, J.},
      author={Shkoller, S.},
       title={A space-time smooth artificial viscosity method with wavelet
  noise indicator and shock collision scheme, {P}art 1: {T}he 1-{$D$} case},
        date={2019},
        ISSN={0021-9991},
     journal={J. Comput. Phys.},
      volume={387},
       pages={81\ndash 116},
         url={https://doi.org/10.1016/j.jcp.2019.02.049},
      review={\MR{3923419}},
}

\bib{RaReSh2019b}{article}{
      author={Ramani, R.},
      author={Reisner, J.},
      author={Shkoller, S.},
       title={A space-time smooth artificial viscosity method with wavelet
  noise indicator and shock collision scheme, {P}art 2: {T}he 2-{$D$} case},
        date={2019},
        ISSN={0021-9991},
     journal={J. Comput. Phys.},
      volume={387},
       pages={45\ndash 80},
         url={https://doi.org/10.1016/j.jcp.2019.02.048},
      review={\MR{3923418}},
}

\bib{RaSh2020}{article}{
      author={Ramani, R.},
      author={Shkoller, S.},
       title={A multiscale model for rayleigh-taylor and richtmyer-meshkov
  instabilities},
        date={2020-03},
     journal={Journal of Computational Physics},
      volume={405},
       pages={109177},
         url={https://doi.org/10.1016/j.jcp.2019.109177},
}

\bib{Re1984}{article}{
      author={Read, K.~I.},
       title={Experimental investigation of turbulent mixing by rayleigh-taylor
  instability},
        date={1984},
     journal={Physica D Nonlinear Phenomena},
      volume={12},
      number={1-3},
       pages={45\ndash 58},
}

\bib{Reisner}{article}{
	author={Reisner, J.},
	title={Private communication}
}

\bib{ReSeSh2013}{article}{
      author={Reisner, J.},
      author={Serencsa, J.},
      author={Shkoller, S.},
       title={A space-time smooth artificial viscosity method for nonlinear
  conservation laws},
        date={2013},
        ISSN={0021-9991},
     journal={J. Comput. Phys.},
      volume={235},
       pages={912\ndash 933},
         url={https://doi.org/10.1016/j.jcp.2012.08.027},
      review={\MR{3017628}},
}

\bib{Ro1931}{article}{
      author={Rosenhead, L.},
       title={The formation of vortices from a surface of discontinuity},
        date={1931},
     journal={Proceedings of the Royal Society of London. Series A, Containing
  Papers of a Mathematical and Physical Character},
      volume={134},
      number={823},
       pages={170\ndash 192},
}

\bib{Ro1956}{article}{
      author={Rott, N.},
       title={Diffraction of a weak shock with vortex generation},
        date={1956},
     journal={Journal of Fluid Mechanics},
      volume={1},
      number={1},
       pages={111\ndash 128},
}

\bib{LeRy1984}{article}{
      author={Ryskin, G.},
      author={Leal, L.~G.},
       title={Numerical solution of free-boundary problems in fluid mechanics.
  part 2. buoyancy-driven motion of a gas bubble through a quiescent liquid},
        date={1984},
     journal={Journal of Fluid Mechanics},
      volume={148},
       pages={19\ndash 35},
}

\bib{Ta1964}{article}{
      author={Takami, H.},
       title={A numerical experiment with discrete vortex approximation with
  reference to the rolling up of a vortex sheet. dept. of aeron. \& astron},
        date={1964},
     journal={Stanford Univ. Report SUDAER},
      volume={202},
}

\bib{Th2012}{article}{
  title={Using multiscale norms to quantify mixing and transport},
  author={Thiffeault, J.L.},
  journal={Nonlinearity},
  volume={25},
  number={2},
  pages={R1},
  year={2012},
  publisher={IOP Publishing}
}

\bib{TrBuEsJuAl2001}{article}{
      author={Tryggvason, G.},
      author={Bunner, B.},
      author={Esmaeeli, A.},
      author={Juric, D.},
      author={Al-Rawahi, N.},
      author={Tauber, W.},
      author={Han, J.},
      author={Nas, S.},
      author={Jan, Y.-J.},
       title={A front-tracking method for the computations of multiphase flow},
        date={2001-05},
     journal={Journal of Computational Physics},
      volume={169},
      number={2},
       pages={708\ndash 759},
         url={https://doi.org/10.1006/jcph.2001.6726},
}

\bib{TrUn1991}{article}{
      author={Unverdi, S.~O.},
      author={Tryggvason, Gr{\'e}tar},
       title={A front-tracking method for viscous, incompressible, multi-fluid
  flows},
        date={1992},
     journal={Journal of computational physics},
      volume={100},
      number={1},
       pages={25\ndash 37},
}

\bib{JaWi2007}{article}{
      author={Wilkinson, J.~P.},
      author={Jacobs, J.~W.},
       title={Experimental study of the single-mode three-dimensional
  Rayleigh-Taylor instability},
        date={2007-12},
     journal={Physics of Fluids},
      volume={19},
      number={12},
       pages={124102},
         url={https://doi.org/10.1063/1.2813548},
}

\bib{Yo1984}{article}{
      author={Youngs, D.~L.},
       title={Numerical simulation of turbulent mixing by Rayleigh-Taylor
  instability},
        date={1984},
     journal={Physica D: Nonlinear Phenomena},
      volume={12},
      number={1-3},
       pages={32\ndash 44},
}

\bib{Yo1989}{article}{
      author={Youngs, D.~L.},
       title={Modelling turbulent mixing by Rayleigh-Taylor instability},
        date={1989},
     journal={Physica D: Nonlinear Phenomena},
      volume={37},
      number={1-3},
       pages={270\ndash 287},
}

\end{biblist}
\end{bibdiv}



\end{document}