as gravitational, which is set to zero when we carry out capillary simulations. ${R}$ is the radius of the drop, $\textbf{e}_r$ is the radial unit vector, and $a$ and $\omega$ are the amplitude and frequency of the oscillatory forcing.

Material properties such as density or viscosity are defined in the entire domain: $$\left.\begin{array}{c}
\rho(\textbf{x},t) = \rho_{1} +\left(  \rho_{2} -\rho_{1}\right)I(\textbf{x},t)\\
\mu(\textbf{x},t) = \mu_{1} +\left( \mu_{2} -\mu_{1}\right)I(\textbf{x},t).
\end{array}\right.
\label{eq:rhomu}$$ The indicator function $I$ in is the Heaviside function, whose value is zero in one phase and one in the other phase. In our discrete numerical implementation $I$ is approximated by $I_{\rm num}$, which represents a smooth transition across 3 to 4 grid cells, as pioneered in the immersed boundary method of [@Pes1977]. In our method $I_{\rm num}$ is generated using a vector distance function computed directly from the tracked interface [@ShJu09].

The fluid variables $\textbf{u}$ and $P$ are calculated by a projection method [@Chorin-mcomp-1968]. The temporal scheme is first order, with implicit time integration used for the viscous terms. For spatial discretization we use the staggered-mesh marker-in-cell (MAC) method [@HW1965] on a uniform finite-difference grid with second-order essentially non-oscillatory (ENO) advection [@SO1989]. The pressure and distance function are located at cell centers 