# Ephemeral Antibubbles: Spatiotemporal Evolution from Direct Numerical
Simulations

Nairita<EMAIL_ADDRESS>Los Alamos National Laboratory, NM 87545,
USA  Rashmi<EMAIL_ADDRESS>TIFR Center for
Interdisciplinary Sciences, Tata Institute of Fundamental Research,
Gopanpally, Hyderabad, 500046, India  Prasad<EMAIL_ADDRESS>TIFR Center for Interdisciplinary Sciences, Tata Institute of Fundamental
Research, Gopanpally, Hyderabad, 500046, India  Rahul
<EMAIL_ADDRESS>also at Jawaharlal Nehru Centre For Advanced
Scientific Research, Jakkur, Bangalore, India. Centre for Condensed Matter
Theory, Department of Physics, Indian Institute of Science, Bangalore 560012,
India

###### Abstract

Antibubbles, which consist of a shell of a low-density fluid inside a high-
density fluid, have several promising applications. We show, via extensive
direct numerical simulations (DNSs), in both two and three dimensions (2D and
3D), that the spatiotemporal evolution of antibubbles can be described
naturally by the coupled Cahn-Hilliard-Navier-Stokes (CHNS) equations for a
binary fluid. Our DNSs capture elegantly the gravity-induced thinning and
breakup of an antibubble via the time evolution of the Cahn-Hilliard scalar
order parameter field $\phi$, which varies continuously across interfaces, so
we do not have to enforce complicated boundary conditions at the moving
antibubble interfaces. To ensure that our results are robust, we supplement
our CHNS simulations with sharp-interface Volume-of-Fluid (VoF) DNSs. We track
the thickness of the antibubble and calculate the dependence of the lifetime
of an antibubble on several parameters; we show that our DNS results agree
with various experimental results; in particular, the velocity with which the
arms of the antibubble retract after breakup scales as $\sigma^{1/2}$, where
$\sigma$ is the surface tension.

###### pacs:

47.27.E-,47.27.eb,47.55.P-

##  I Introduction

Antibubbles, which comprise a shell of a low-density fluid inside a high-
density fluid, have been known for close to $90$ years, since the work of
Hughes and Hughes Hughes and Hughes (1932). In contrast to an ordinary bubble,
an antibubble has two surfaces, which trap a certain volume of fluid between
them. Therefore, the contact area of an antibubble is much larger than that of
a bubble with the same fluid volume; this property can be exploited for a
variety of chemical reactions. Furthermore, antibubbles are commonly used in
clinical diagnostic imaging, sonoporation (see, e.g., Ref. Kotopoulis  _et
al._ (2014)), as agents for ultrasound-guided drug delivery Johansen  _et al._
(2015a); Kotopoulisa  _et al._ (2015); Postema and Gilja (2007), and for
active leakage detection Johansen  _et al._ (2015b). Clearly, an antibubble is
unstable in the presence of gravity, and, if the inner core of the antibubble
is denser than the outer core, the antibubble rises under gravity; because of
hydrostatic pressure in the outer core, the fluid rises from the bottom to the
top, resulting in a thinning, and subsequent collapse, of the outer shell.
Although there have been a number of experimental investigations of the
spatiotemporal evolution of an antibubble Dorbolo  _et al._ (2006); Kim and
Vogel (2006); Dorbolo  _et al._ (2010, 2003); Bai  _et al._ (2016); Kalelkar
(2017); Chatzigiannakis  _et al._ (2021); Chatzigiannakis and Vermant (2020);
Song  _et al._ (2020) and drops Will  _et al._ (2021); Li  _et al._ (2021);
Soto  _et al._ (2020a, b); Prakash  _et al._ (2016) to name a few, over the
past few decades, theoretical studies of antibubble evolution have been
initiated only recently Scheid  _et al._ (2012); Zou  _et al._ (2013);
Sob’yanin (2015); Yang  _et al._ (2020); Vitry  _et al._ (2019); Song  _et
al._ (2020) and they have not, attempted, hitherto, to address the
spatiotemporal evolution of antibubbles in detail. The number of experimental
studies of the complete spatiotemporal evolution of antibubbles is also
limited, partly because great care has to be exercised to stabilise
antibubbles; often, surfactant molecules have to be introduced into the high-
density liquid phase for such stabilization.

We develop a natural, multiphase model for antibubbles and demonstrate how we
can use direct numerical simulations (DNSs) to follow the spatiotemporal
development of ephemeral, but beautiful, antibubbles. We show that the Cahn-
Hilliard-Navier-Stokes (CHNS) equations, which have been used to study a
variety of problems in binary-fluid flows Boffetta and Ecke (2012); Perlekar
_et al._ (2014); Scarbolo  _et al._ (2013); Pal  _et al._ (2016); Perlekar
_et al._ (2017); Gibbon  _et al._ (2016, 2018); Fan  _et al._ (2016, 2017,
2018), provide a minimal theoretical framework for studying the spatiotemporal
evolution of antibubbles; in addition to a velocity field, the CHNS system
employs a phase field $\phi$ that distinguishes between the two fluid phases.
We use the CHNS equations Boffetta and Ecke (2012); Perlekar  _et al._ (2014);
Pal  _et al._ (2016); Gibbon  _et al._ (2016); Perlekar  _et al._ (2017);
Gibbon  _et al._ (2018) to study antibubbles in two-dimensions (2D) and in
three dimensions (3D) by using extensive DNSs. To complement our CHNS results,
we also employ an alternative Volume-of-Fluid (VoF) numerical scheme that is a
sharp-interface method.

Our studies yield several interesting results that (a) provide explanations
for many experimental observations and (b) suggest new experimental studies:
Our results for the spatiotemporal evolution of antibubbles, in 2D and 3D,
show clearly how antibubble breakup occurs either because of gravity-induced
thinning or the puncturing of its bottom boundary; the the collapsing
antibubble then forms a rim that retracts with a velocity $v_{rim}$. We
uncover signatures of this collapse in Fourier-space spectra of the Cahn-
Hilliard field $\phi$, the velocity, and the vorticity. We show that
$v_{rim}\sim\sigma^{1/2}$, where $\sigma$ is the surface tension; and this
power-law exponent is independent of the kinematic viscosity $\nu$ of the
background fluid. We investigate the dependence of the scaled antibubble
lifetime $\tau_{1}/\tau_{g}$ on $\sigma$, $\nu$, and the scaled outer radius
of the antibubble $R_{0}/h_{0}$, in 2D; here, $\tau_{g}\equiv\sqrt{R_{0}/Ag}$,
with $A$ the Atwood number and $g$ the acceleration due to gravity, and
$h_{0}$ is the initial thickness of the antibubble shell. We compare our
results with experiments and earlier theoretical studies.

The remainder of this paper is organized as follows: In Sec. II we present the
models and the numerical methods that we use. In Sec. III we present the
results our DNSs in the following subsections:

  * •

subsection III.1 on the spatiotemporal evolution of antibubbles (2D DNSs);

  * •

subsection III.2 on the temporal evolution of Fourier spectra of $\phi$ and
the velocity and vorticity fields that are associated with antibubbles in our
2D DNSs;

  * •

subsection III.3 on the velocity of the retracting rim of the ruptured
antibubble in our 2D DNSs;

  * •

subsection III.4 on the signatures of this rupture in the time dependence of
the energy;

  * •

subsection III.5.1 on the dependence of the antibubble-rupture time on the
size of the antibubble and the surface tension;

  * •

subsection III.5.2 on the dependence of the antibubble-rupture time on the
kinematic viscosity;

  * •

subsection III.6 on the spatiotemporal evolution of antibubbles in our 3D
DNSs.

We present conclusions and a discussion of our results in Sec. IV.

##  II Model and Numerical Methods

An antibubble consists of a shell of light fluid of density $\rho_{1}$ inside
a background heavy fluid of density $\rho_{2}$; these fluids are immiscible
and incompressible. The following CHNS equations provide a natural theoretical
description for such a binary-fluid mixture Celani  _et al._ (2009); Scarbolo
and Soldati (2013); Yue  _et al._ (2004); Scarbolo  _et al._ (2011):

$\displaystyle\rho(\phi)D_{t}{\bm{u}}$ $\displaystyle=-\nabla
p+\varsigma(\phi)\nabla^{2}{\bm{u}}+\mathbf{F}_{\sigma}$
$\displaystyle+[\rho(\phi)-\displaystyle\frac{\rho_{1}+\rho_{2}}{2}]{\bm{g}}-\alpha\mathbf{u};$
(1) $\displaystyle D_{t}\phi$
$\displaystyle=\gamma\nabla^{2}{\mu}~{}{\rm{and}}~{}\nabla\cdot{\bm{u}}=0;$
(2)

In Eqs.(1)-(2), $D_{t}\equiv\left(\partial_{t}+{\bm{u}}\cdot\nabla\right)$ is
the convective derivative, ${\bm{u}}\equiv(u_{x},u_{y},u_{z})$ is the fluid
velocity, $\alpha$ is the coefficient of friction (we use this only in 2D),
${\bm{g}}$ is the constant acceleration due to gravity, which points
downwards, $\rho(\phi)\equiv\rho_{1}(1+\phi)/2+\rho_{2}(1-\phi)/2$ is the
density, $\varsigma(\phi)\equiv\nu_{1}(1+\phi)/2+\nu_{2}(1-\phi)/2$ is the
dynamic viscosity, $\phi({\bm{x}},t)$ is the order-parameter field at the
point ${\bm{x}}$ and time $t$ [with $\phi({\bm{x}},t)>0$ in the background
(majority) phase and $\phi({\bm{x}},t)<0$ in the antibubble-shell (minority)
phase]. For the illustrative results that we present, we choose $\rho(\phi)$
and $\varsigma(\phi)$ such that the kinematic viscosity $\nu$ is identical in
the two phases.

When we study 2D flows, we use the following stream-function-vorticity
formulation Boffetta and Ecke (2012); Pandit  _et al._ (2009):

$\displaystyle\left(\partial_{t}+{\bm{u}}\cdot\nabla\right)\omega$
$\displaystyle=$
$\displaystyle\nu\nabla^{2}\omega-\alpha\omega-\nabla\times(\phi\nabla\mu)-A\nabla\phi\times{\bm{g}};$
$\displaystyle\left(\partial_{t}+{\bm{u}}\cdot\nabla\right){\bf\phi}$
$\displaystyle=$
$\displaystyle\gamma\nabla^{2}{\mu}~{}{\rm{and}}~{}\nabla\cdot{\bm{u}}=0.$ (3)

Here, $\omega=(\nabla\times{\bm{u}}){\hat{\bm{e}}}_{z}$ is the vorticity. We
obtain ${\omega}({\bm{x}},t)$ and $\phi({\bm{x}},t)$ from our 2D DNS sup , and
from these we calculate the total kinetic energy
$E(t)=\langle|{\bm{u}}({\bm{x}},t)|^{2}\rangle_{{\bm{x}}}$ and the fluid-
energy dissipation rate
$\varepsilon(t)=\langle\nu|\omega({\bm{x}},t)|^{2}\rangle_{{\bm{x}}}$, where
$\langle\rangle_{{\bm{x}}}$ denotes the average over space, and
$\omega=\nabla\times{\bm{u}}$ is the fluid vorticity. We calculate the
lifetime of an antibubble falling under gravity in two ways: (a) from the
energy time series; and (b) from the minimum thickness at the South pole of
the antibubble.

At time $t=0$ we begin with the order-parameter profile Celani  _et al._
(2009)

$\displaystyle\phi({\bm{x}})$ $\displaystyle=$
$\displaystyle-\tanh\left[\frac{1}{\sqrt{2}\xi}\left(|{\bm{x}}-{\bm{x}}_{c}|-R_{0}\right)\right]$
(4)
$\displaystyle+\tanh\left[\frac{1}{\sqrt{2}\xi}\left(|{\bm{x}}-{\bm{x}}_{c}|-R_{1}\right)\right]-1.0,$

where $R_{0}$ and $R_{1}$ are, respectively, the initial magnitudes of the
outer and inner radii of the antibubble, whose centre is initially at
${\bm{x}}_{c}$; the interface width $\xi$ is measured by the dimensionless
Cahn number $Ch=\xi/L$. We use the initial thickness of the antibubble shell
$h_{0}\equiv R_{0}-R_{1}$ as a typical length scale.

In our CHNS studies we use the Boussinesq approximation, wherein
$\rho(\phi)D_{t}{\bm{u}}\approx\displaystyle\frac{\rho_{1}+\rho_{2}}{2}D_{t}{\bm{u}}$
and
$[\rho(\phi)-\displaystyle\frac{\rho_{1}+\rho_{2}}{2}]{\bm{g}}\approx-A\phi g$
in Eq.(1); the surface-tension force
$\mathbf{F}_{\sigma}\equiv-\phi\boldsymbol{\nabla}\mu$ and the chemical
potential $\mu({\bm{x}},t)$ follow from the Cahn-Hilliard free-energy
functional ${\mathcal{F}}$ (see the Supplemental Material sup ):

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\delta{\mathcal{F}}[\phi]/\delta\phi({\bm{x}},t);$
$\displaystyle{\mathcal{F}}[\phi]$ $\displaystyle=$
$\displaystyle\Lambda\int[(\phi^{2}-1)^{2}/(4\xi^{2})+|\nabla\phi|^{2}/2]d{\bm{x}};$
(5)

here, $\Lambda$ is the energy density with which the two phases mix in the
interfacial regime Celani  _et al._ (2009), $\xi$ sets the scale of the
interface width , $\sigma=\frac{2\sqrt{2}\Lambda}{3\xi}$ is the surface
tension, $\gamma$ is the mobility Yue  _et al._ (2004) of the binary-fluid
mixture; the Cahn number $Ch=\xi/L$, where $L$ is the linear size of our
simulation domain, is a dimensionless measure of the interfacial width; the
Schmidt number $Sc=\nu/D$. We assume that $\gamma$ and $\rho$ are independent
of $\phi$; we keep the diffusivity $D=\frac{\gamma\Lambda}{\xi^{2}}$ constant.
We study the retraction dynamics of an antibubble by varying the Atwood number
$A=(\rho_{2}-\rho_{1})/(\rho_{2}+\rho_{1})$ and the Bond number $Bo=Ag\rho
h_{0}^{2}/\sigma$. The Bond number is a convenient dimensionless ratio of body
forces (gravity), on the antibubble, and the surface or interfacial tension
$\sigma$. We express times in multiples of
$\tau_{g}\equiv\sqrt{R_{0}/Ag}=1.16$.

We solve the CHNS Eqs.(1)-(2) by using the pseudospectral method with periodic
boundary conditions; because of the cubic nonlinearity in the chemical
potential $\mu$, we use $N/2$-dealiasing Canuto  _et al._ (2012). We use the
exponential Adams-Bashforth method ETD2 Cox and Matthews (2002) for time
marching. For our 2D DNSs, we use computers with Graphics Processing Units
(e.g., the NVIDIA K80), which we program in CUDA
http://www.nvidia.com/object/cuda$\\_$home$\\_$new.html. . We use conventional
CPU-based computers for our 3D DNSs. We have $N$ collocation points in each
direction, and our domain length is $L=2\pi$. Our efficient code allows us to
explore the CHNS parameter space and carry out very long simulations that are
essential for our studies.

For the VoF method we note the following points: It is a sharp-interface
method, so it has $Ch=0$. The surface tension force for VoF is
$\mathbf{F}_{\sigma}\equiv\sigma\kappa\delta_{s}\mathbf{n}$, where $\kappa$ is
the curvature, $\delta_{s}$ is the Dirac delta function on the interface, and
$\mathbf{n}$ is the normal to the interface. We use Basilisk
http://basilisk.fr. ; Popinet (2018), an open-source solver for carrying out
both 2D and 3D axisymmetric VoF DNSs. Basilisk, which employs the Bell-
Collela-Glaz advection scheme Bell  _et al._ (1989) and the implicit viscosity
solver http://basilisk.fr. ; Popinet (2018), is parallelized over conventional
CPUs http://basilisk.fr. ; Popinet (2018). The VoF solver does not invoke the
Boussinesq approximation; and it can handle large density and viscosity
contrasts. The breakup of an interface is sensitive to the resolution in the
VoF, so we use an initial configuration in which the antibubble is already
punctured at the bottom; and we then investigate its spatiotemporal evolution.
This initial condition is similar to that used in experiments and in Ref.
Sob’yanin (2015). We use the following boundary conditions in our VoF
simulations: (a) In our 2D DNSs we employ periodic boundary conditions in all
directions; (b) in the 3D axisymmetric case we use an axisymmetric boundary
condition on the $z=0$ axis and the no-slip condition ${\bm{u}}=0$ at other
boundaries.

We list the parameters for some of our representative DNS runs in Tables 1, 2,
3, and 4. We also give Tables with the details of all our DNSs in the
Supplemental Material sup . As we show below, by utilising both the CHNS
framework and the VoF method, we can capture accurately the rupture and
retraction dynamics of an antibubble in an elegant and numerically efficient
way.

| $N$ | $R_{0}/\xi$ | $R_{1}/\xi$ | $R_{0}/h_{0}$ | $h_{0}/\xi$ | $D$ | $Ch$ | $Bo(\times 10^{-3})$ | $\nu$ | $Sc$
---|---|---|---|---|---|---|---|---|---|---
R1 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $36$ | $0.007$ | $1.75$
R2 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $12$ | $0.007$ | $1.75$
R3 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $8.9$ | $0.007$ | $1.75$
R4 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $5.0$ | $0.007$ | $1.75$
R5 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $4.0$ | $0.007$ | $1.75$
R6 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $4.5$ | $0.007$ | $1.75$
R7 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $3.3$ | $0.007$ | $1.75$
R8 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $2.7$ | $0.007$ | $1.75$
R9 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $2.4$ | $0.007$ | $1.75$
R10 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $2.1$ | $0.007$ | $1.75$
R11 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.9$ | $0.007$ | $1.75$
R12 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.7$ | $0.007$ | $1.75$
R13 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.5$ | $0.007$ | $1.75$
R14 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.4$ | $0.007$ | $1.75$
R15 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.3$ | $0.007$ | $1.75$
R16 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $1.1$ | $0.007$ | $1.75$
R17 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $0.8$ | $0.007$ | $1.75$
R18 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $0.73$ | $0.007$ | $1.75$
R19 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $0.56$ | $0.007$ | $1.75$
R20 |  $1024$ (2D) | $76.6$ | $69.7$ | $7$ | $6.96$ | $0.004$ | $0.0028$ | $0.45$ | $0.007$ | $1.75$
R21 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $12$ | $0.007$ | $1.75$
R22 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $5.0$ | $0.007$ | $1.75$
R23 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $4.0$ | $0.007$ | $1.75$
R24 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $3.3$ | $0.007$ | $1.75$
R25 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $3.1$ | $0.007$ | $1.75$
R26 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $2.7$ | $0.007$ | $1.75$
R27 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $2.4$ | $0.007$ | $1.75$
R28 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $1.9$ | $0.007$ | $1.75$
R29 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $1.7$ | $0.007$ | $1.75$
R30 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $1.5$ | $0.007$ | $1.75$
R31 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $1.3$ | $0.007$ | $1.75$
R32 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $1.1$ | $0.007$ | $1.75$
R33 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $0.8$ | $0.007$ | $1.75$
R34 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $0.73$ | $0.007$ | $1.75$
R35 |  $1024$ (2D) | $76.6$ | $69.7$ | $8$ | $6.96$ | $0.004$ | $0.0028$ | $0.45$ | $0.007$ | $1.75$

Table 1: The parameters $R_{0},R_{1},h_{0},\xi$, $\sigma$, $\nu,D,CH,Sc$ for
our CHNS DNS runs R1-R35. Initially, $R_{0}$ and $R_{1}$ are, respectively,
the outer and inner radii of the antibubble; $h_{0}\equiv(R_{0}-R_{1})$;
$\sigma$ is the surface tension; $\nu$ is the kinematic viscosity. In all the
runs, the Atwood number times the acceleration due to gravity $Ag=0.99$; the
number of collocation points is $N^{2}$; the diffusivity is $D$; the Cahn
number $Ch=\xi/L$, $\xi$ being the interface width. The typical values of
$R_{0}/h_{0}$ used in experiments are about $38$ Zou  _et al._ (2013). For all
our 2D runs, the friction coefficient $\alpha=0.001$.

| $N$ | $R_{0}/\xi$ | $R_{1}/\xi$ | $R_{0}/h_{0}$ | $h_{0}/\xi$ | $D$ | $Ch$ | $Bo(\times 10^{-3})$ | $\nu$ | $Sc$
---|---|---|---|---|---|---|---|---|---|---
R36 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $12$ | $0.007$ | $1.75$
R37 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $5.0$ | $0.007$ | $1.75$
R38 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $4.0$ | $0.007$ | $1.75$
R39 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $3.3$ | $0.007$ | $1.75$
R40 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $3.1$ | $0.007$ | $1.75$
R41 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $2.7$ | $0.007$ | $1.75$
R42 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $2.4$ | $0.007$ | $1.75$
R43 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $1.9$ | $0.007$ | $1.75$
R44 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $1.7$ | $0.007$ | $1.75$
R45 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $1.5$ | $0.007$ | $1.75$
R46 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $1.3$ | $0.007$ | $1.75$
R47 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $1.1$ | $0.007$ | $1.75$
R48 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $0.8$ | $0.007$ | $1.75$
R49 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $0.73$ | $0.007$ | $1.75$
R50 |  $1024$ (2D) | $76.6$ | $69.7$ | $9$ | $6.96$ | $0.004$ | $0.0028$ | $0.45$ | $0.007$ | $1.75$
R51 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.007$ | $1.75$
R52 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.01$ | $2.6$
R53 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.01$ | $3.9$
R54 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.02$ | $5.0$
R55 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.02$ | $6.8$
R56 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.035$ | $8.75$
R57 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.035$ | $10.0$
R58 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.05$ | $12.5$
R60 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.05$ | $14.0$
R61 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.05$ | $15.75$
R62 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.05$ | $17.5$
R63 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.08$ | $20$
R64 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.12$ | $30$
R65 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.18$ | $45$
R66 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.18$ | $58$
R67 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.27$ | $67.5$
R68 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.27$ | $80.0$
R69 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.4$ | $100$
R70 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.6$ | $151.5$
R71 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.01$ | $0.96$ | $241.7$
R72 |  $1024$ (2D) | $76.6$ | $69.7$ | $11$ | $6.96$ | $0.004$ | $0.0028$ | $0.1$ | $0.96$ | $241.7$

Table 2: The parameters $R_{0},R_{1},h_{0},\xi$, $\sigma$, $\nu,D,CH,Sc$ for
our CHNS DNS runs R36-R72. Initially, $R_{0}$ and $R_{1}$ are, respectively,
the outer and inner radii of the antibubble; $h_{0}\equiv(R_{0}-R_{1})$;
$\sigma$ is the surface tension; $\nu$ is the kinematic viscosity. In all the
runs, the Atwood number times the acceleration due to gravity $Ag=0.99$, the
number of collocation points is $N^{2}$; the diffusivity is $D$; the Cahn
number $Ch=\xi/L$, $\xi$ being the interface width. The typical values of
$R_{0}/h_{0}$ used in experiments are about $38$ Zou  _et al._ (2013). For all
our 2D runs, the friction coefficient $\alpha=0.001$.

| $N$ | $R_{0}/\xi$ | $R_{1}/\xi$ | $R_{0}/h_{0}$ | $h_{0}/\xi$ | $D$ | $Ch$ | $Bo(\times 10^{-3})$ | $\nu$ | $Sc$
---|---|---|---|---|---|---|---|---|---|---
Q1 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $1.3$ | $0.0016$ | $0.27$
Q2 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $0.63$ | $0.0016$ | $0.27$
Q3 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $0.4$ | $0.0016$ | $0.27$
Q4 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $0.3$ | $0.0016$ | $0.27$
Q5 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $0.25$ | $0.0016$ | $0.27$
Q6 |  $256$ (3D) | $24.0$ | $12.0$ | $1.2$ | $12.0$ | $0.006$ | $0.013$ | $0.21$ | $0.0016$ | $0.27$
Q7 |  $512$ (3D) | $14$ | $10.5$ | $4.0$ | $3.5$ | $0.006$ | $0.011$ | $2.7$ | $0.0016$ | $0.27$

Table 3: The parameters $R_{0},R_{1},h_{0},\xi$, $\sigma$, $\nu,D,CH,Sc$ for
our CHNS DNS runs Q1-Q7. Initially, $R_{0}$ and $R_{1}$ are, respectively, the
outer and inner radii of the antibubble; $h_{0}\equiv(R_{0}-R_{1})$; $\sigma$
is the surface tension; $\nu$ is the kinematic viscosity. In all the runs, the
Atwood number times the acceleration due to gravity $Ag=0.5$, the number of
collocation points is for Q1-Q6 is $N^{3}=256^{3}$ and for Q7 is
$N^{3}=512^{3}$: the diffusivity is $D$; the Cahn number $Ch=\xi/L$, with
$\xi$ being the interface width. The typical values of $R_{0}/h_{0}$ used in
experiments are about $38$ Zou  _et al._ (2013).

| $R_{0}$ | $R_{1}$ | $h_{0}$ | $A$ | $g$ | $\nu$ | $Bo(\times 10^{-3})$
---|---|---|---|---|---|---|---
PR1 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $0.9$
PR2 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $1.17$
PR3 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $1.82$
PR4 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $4.1$
PR5 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $8.2$
PR6 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $16.4$
APR1 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $0.75$
APR2 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $0.9$
APR3 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $1.17$
APR4 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $1.82$
APR5 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $2.05$
APR6 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $4.1$
APR7 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $8.2$
APR8 | $0.859$ | $0.739$ | $0.12$ | $0.9$ | $1.$ | $0.007$ | $16.4$
APR9 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $0.9$
APR10 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $1.82$
APR11 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $2.05$
APR12 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $4.1$
APR13 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $8.2$
APR14 | $0.859$ | $0.739$ | $0.12$ | $0.01$ | $100.$ | $0.007$ | $16.4$

Table 4:  The parameters $R_{0},R_{1},h_{0}$, $A$, $g$, $\nu$ and $Bo$ for our
VoF runs: 2D(${\tt{PR1}}-{\tt{PR6}}$); and 3D axisymmetric
(${\tt{APR1}}-{\tt{APR14}}$). For these runs $\rho_{2}=1$. For the 2D and
axisymmetric 3D runs we use a square box of area $L^{2}=16$ and discretize it
with $N^{2}=2048^{2}$ collocation points. For our 2D VoF DNSs we employ
periodic boundary conditions in all directions; in 3D we use an axisymmetric
boundary condition on the $z=0$ axis and the no-slip boundary condition
${\bm{u}}=0$ at other boundaries.

##  III Results

We now present the results of our DNSs in the following subsections: In
subsection III.1 we illustrate the spatiotemporal evolution of 2D antibubbles
and, in subsection III.2, we give the temporal evolution of Fourier spectra of
$\phi$ and the velocity and vorticity fields that are associated with these
antibubbles. In subsection III.3 we explore the dependence on $\sigma$ of the
velocity of the retracting rim of a ruptured 2D antibubble; and, in subsection
III.4, we examine the signatures of this rupture in the time dependence of the
energy. Furthermore, we elucidate the dependence of the antibubble-rupture
time on the size of the antibubble and the surface tension (subsection
III.5.1) and on the kinematic viscosity (subsection III.5.2). We then present
illustrative results on the spatiotemporal evolution of 3D antibubbles in
subsection III.6

Figure 1:  Pseudocolor plots of the CHNS $\phi$ field with antibubble radius
$R_{0}/h_{0}=7$ at time (a) $t=0.05\tau_{g}$, (b) $t=3.26\tau_{g}$, (c)
$t=11.77\tau_{g}$ for $\nu=0.007$ and $\sigma=16.6$ (i.e, $Bo=0.001$, from our
2D DNS run R16); pseudocolor plots of the corresponding $\omega$ field
(d)–(f). These plots show that, initially, the antibubble rises under gravity,
breaks, (a,d), then it retracts (b,e), becomes a circular droplet and rises
again (e,f). Pseudocolor plots of the volume fraction for the 2D VoF run PR5
at time (g) $t=0$, (h) $t=0.1\tau_{g}$, (i) $t=0.3\tau_{g}$; (j), (k), and (l)
pseudocolor plots of the vorticity for the 2D VoF PR5 run at the same times as
in (g), (h), and (i), respectively; in these 2D VoF runs, $\sigma=1.68$; these
plots show that the antibubble retracts and forms a rim. Time is measured in
units of $\tau_{g}$.

###  III.1 Spatiotemporal evolution: 2D

In Fig. 1 we display, via pseudocolor plots of $\phi$ and $\omega$, the
spatiotemporal evolution of a 2D antibubble rising under gravity in our CHNS
DNS; blue and red indicate heavy and light fluids, respectively. Figure 1(a)
shows the pseudocolor plot of $\phi$ for the antibubble at a time $\simeq
1.07\times 0.05\tau_{g}$. Initially, the antibubble rises because of gravity,
breaks (Fig. 1(b)), then displays arms that retract, forms a disc-shaped
droplet, which rises under gravity (Fig. 1(c)), because the antibubble-shell
(red) fluid is lighter than the background (blue) fluid. Figures 1(d)-(f) show
pseudocolor plots of the vorticity field at the same times as their $\phi$
counterparts in Figs. 1(a)-(c). The spatiotemporal evolution of such a 2D
antibubble is shown in the Videos LABEL:sup_video1 and LABEL:sup_video2 in the
Supplemental Material sup . Note that, just after forming the disc-type shape,
the antibubble goes down in a direction opposite to that dictated by gravity
(Fig. 1(e)), because the surface tension is high, so, during the retraction of
the arms, when interfacial energy is released, the retracting droplet can be
pushed in this opposite direction. Eventually this backward thrust is damped
by viscosity and, finally, the (light) droplet rises under gravity. We begin
with $\omega=0$, everywhere. We see that, because of the
$\nabla\times\phi\nabla\mu$ term on the right-hand side (RHS) of Eq.(3), a
vorticity field is generated at the outer and inner surfaces of the
antibubble, i.e., wherever $\phi$ changes in space. The vortical regions that
are generated initially have a size that is comparable to the interface width.
These vortices then grow with time, until they are damped by viscosity.
Similarly, the gravity term is also significant wherever $\phi$ changes sign.
The competition between the viscous, surface-tension, and gravity terms in the
CHNS equations determines the breakup time of an antibubble, as we show below.

In Figs. 1(g)-(l) we use pseudocolor plots to illustrate antibubble retraction
in our VoF DNS in 2D (run PR1). Both our CHNS and VoF DNSs show that, as it
ruptures, an antibubble forms a rim, which then retracts.

(a)(b)(c)(d)(e)(f)(g)(h)(i)

Figure 2:  Log–log (base 10) plots of the spectra (Eq. 6) (a)–(c) $S(k,t)$,
(d)–(f) $E(k,t)$, and (g)–(i) $\Omega(k,t)$ versus $k$, with $k_{0}=1$ and
time in units of $\tau_{g}$.

###  III.2 Fourier spectra

Fourier-space spectra give us a complementary view of the spatiotemporal
evolution of the rupture of an antibubble (cf., the physical-space pseudocolor
plots in Figs. 1 (a)-(l)). In Fig. 2 we show, for our 2D CHNS DNSs, log-log
(base 10) plots of the spectra for the field $\phi$, the energy spectrum
$E(k,t)$, and the enstrophy spectrum $\Omega(k,t)$, at different times during
this evolution. These spectra are defined as follows:

$\displaystyle S(k,t)$ $\displaystyle\equiv$
$\displaystyle\sum\limits_{k-\frac{1}{2}\leq k^{\prime}\leq
k+\frac{1}{2}}\langle\hat{\phi}({\bm{k}}^{\prime},t)\hat{\phi}(-{\bm{k}}^{\prime},t)\rangle;$
$\displaystyle E(k,t)$ $\displaystyle\equiv$
$\displaystyle\displaystyle\frac{1}{2}\sum\limits_{k-\frac{1}{2}\leq
k^{\prime}\leq
k+\frac{1}{2}}\langle\hat{\bm{u}}({\bm{k}}^{\prime},t)\hat{\bm{u}}(-{\bm{k}}^{\prime},t)\rangle;$
$\displaystyle\Omega(k,t)$ $\displaystyle\equiv$
$\displaystyle\displaystyle\frac{1}{2}\sum\limits_{k-\frac{1}{2}\leq
k^{\prime}\leq
k+\frac{1}{2}}\langle\hat{{\omega}}({\bm{k}}^{\prime},t)\hat{{\omega}}(-{\bm{k}}^{\prime},t)\rangle;$
(6)

here, the circumflex denotes a spatial Fourier transform and $k$ and
$k^{\prime}$ are, respectively, the magnitudes of the wave vectors ${\bm{k}}$
and ${\bm{k}}^{\prime}$. At early times (Figs. 2 (a), (d), and (g)), these
spectra, especially $S(k,t)$, show oscillations with small and large periods,
which are, respectively, related inversely to the outer radius of the
antibubble and the thickness of its shell. As the antibubble ruptures (Figs. 2
(b), (e), and (h)), it loses its circular shape, and, in turn, the spectra
lose their oscillations. Finally, the antibubble is replaced by a single
droplet (Fig. 1(c)) of the light fluid, which rises under gravity; at this
stage, these spectra, especially $S(k,t)$, show oscillations with a small
period, which is related inversely to the radius of this rising droplet. These
Fourier-space spectra show that the breakup of an antibubble leads to
turbulence, insofar as $E(k,t)$, $S(k,t)$, and $\Omega(k,t)$ extend
significantly over several orders of magnitude of $k$. The Taylor-microscale
Reynolds number $Re_{\lambda}$ is $0.54$, $128.5$, and $163.4$ in Figs. 2 (a),
(b), and (c), respectively; here, $Re_{\lambda}=u_{rms}\lambda/\nu$, where the
root-mean-square velocity $u_{rms}(t)=2\sum\limits_{k}E(k,t)$ and the Taylor
microscale $\lambda\equiv[\sum k^{2}E(k,t)/\sum E(k,t)]^{-1/2}$.

Figure 3:  Plot of the retraction velocity $v_{rim}$ versus the theoretical
estimate $v_{th}=\sqrt{\sigma/\rho_{2}a}$ for our 2D CHNS runs R1-R17, 2D VoF
runs PR1-PR6, and 3D axisymmetric VoF runs APR1-APR14; red circles show data
from our CHNS runs; green, yellow, and blue markers indicate data from our VoF
runs. The black line represents the theoretical line. Here the surface tension
$\sigma=\displaystyle\frac{Ag\rho h_{0}^{2}}{Bo}$, where $Bo$ is the Bond
number.

(a)(b)(c)(d)(e)(f)

Figure 4: Plots of (a) the total kinetic energy $E(t)$ and (b) the time
derivative $dE(t)/dt$ versus $t/\tau_{g}$ for the 2D CHNS runs R16-R20; (c)
plots versus $Bo$ of the antibubble lifetimes $\tau_{1}$ and $\tau_{2}$
calculated, respectively, from the energy time series (blue line with circles)
and the minimum antibubble thickness (orange line with squares), for the 2D
CHNS runs R1-R16. (d) Plots versus $Bo$ (for $Sc=1.75$), of
$\tau_{1}/\tau_{g}$ for $R_{0}/h_{0}=7$ (DNS run R1-R20, blue line with
circles), $R_{0}/h_{0}=8$ (DNS run R21-R35, orange line with crosses),
$R_{0}/h_{0}=9$ (DNS run R36-R50, green line with squares); (e) log-log plot
of $\tau_{1}/\tau_{g}$ versus $Sc$, over three orders of magnitude and at a
fixed value of $Bo$ (runs R51-R71); (f) plots of $h_{min}$ versus
$t/\tau_{g}$, at $Bo=0.01$ and $A=0.01$ and $g=99$, for $Sc=1.75$ (R51, blue
line), $Sc=1.75$ (R52, orange line) $Sc=20$ (R63, green line), and $Sc=151.5$
(R70, red line).

###  III.3 Rim-retraction velocity

When an antibubble ruptures, the surface-tension energy (the interfacial free
energy in the CHNS description) is converted to the kinetic energy of the
antibubble. The rate of change of the former is $\sigma d{\mathcal{S}}(t)/dt$,
where ${\mathcal{S}}(t)$ is the outer perimeter of the antibubble in our 2D
DNSs; and the kinetic energy that is released is $E_{M}=\int_{S_{f}}\left(\rho
v^{2}({\bm{x}},t)\right)d{\bm{x}}$, where $S_{f}$ is the area of the majority
phase surrounding the minority phase; this causes the film to retract. Recall
that Figs. 1(a)-(l) show that, in both our CHNS and VoF DNSs, as it ruptures,
an antibubble forms a rim, which then retracts. When this rim retracts, the
outer perimeter of the antibubble reduces by a length $R_{f}$ on one side
after the rupture; and, if $a$ is the thickness of the vanished film, then
$E_{M}=\rho v_{rim}^{2}aR_{f}$, where $v_{rim}$ is the rim-retraction
velocity. If we assume, furthermore, that the rim retracts with a constant
velocity, then $dv_{rim}/dt=0$, so $dR_{f}/dt=v_{rim}$. Finally, we find

$\sigma d{\mathcal{S}}(t)/dt\propto
dE_{M}/dt,~{}\Rightarrow~{}v_{rim}=\left(\sigma/\rho\right)^{1/2}.$ (7)

This dependence of $v_{rim}$ on $\sigma$ matches the experimental observations
in Ref. Scheid  _et al._ (2012). Moreover, this dependence is the same in both
our CHNS and VoF DNSs, even though, in the former, antibubble breakup occurs
because of gravity-induced thinning, whereas, in the latter, this breakup is
induced by puncturing the antibubble at its bottom boundary. The puncturing
initial condition has also been used in the experiments and a theory of
antibubble collapse Sob’yanin (2015); Zou  _et al._ (2013).

The plot in Fig. 3(b) shows that rim velocities, which we obtain from 2D CHNS
and VoF DNSs and our 3D axisymmetric VoF runs, for both high and low $A$, are
in excellent agreement with the theoretical prediction
$v_{th}\sim\sqrt{\sigma/\rho_{2}a}$ Sob’yanin (2015), where $a$ is the radius
of the rim ($\sqrt{h_{0}R_{0}/\pi})$.

###  III.4 Energy time series

When an antibubble bursts, the surface tension energy is converted into the
kinetic energy of the fluid; this yields a spike in the fluid-energy time
series (see Fig. 4(a)). Therefore, we identify the breakup time $\tau_{1}$ as
the instant at which $dE(t)/dt$ displays a maximum. We also define the breakup
time $\tau_{2}$ at which the antibubble-shell thickness $h_{min}$, at the
lower end of the antibubble, vanishes. In Fig. 4(b) we show plots of both
$\tau_{1}$ and $\tau_{2}$ versus $Bo$ for an antibubble with $R_{0}/h_{0}=11$
and $\nu=0.967$. This plot shows that both our estimates of the collapse
times, $\tau_{1}$ and $\tau_{2}$, agree with each other. (For details, see the
Supplemental Material sup .)

###  III.5 Antibubble-breakup times

The antibubble-breakup times, which we have defined above, depend on the
initial size of the antibubble, the surface tension, and the kinematic
viscosity. We explore these dependences below via our 2D CHNS DNSs.

####  III.5.1 Dependences on Size and Surface Tension

The size- and surface-tension dependences of the antibubble-breakup time
follow from the plots of $\tau_{1}/\tau_{g}$ versus $Bo$, which we give in
Fig. 4(d), for (i) $R_{0}/h_{0}=9$ (DNS runs R13-R16, pink line with
triangles), (ii) $R_{0}/h_{0}=8$ (DNS run R17-R20, red line with squares), and
(iii) $R_{0}/h_{0}=9$ (DNS run R21-R24, blue line with circles). In all these
plots, the viscosity is fixed and the Schmidt number $Sc=1.75$. Although the
scaling of $\tau_{1}$ by $\tau_{g}=\sqrt{R/Ag}$ makes all the curves collapse
on top of each, to a large degree, there is a small, but noticeable,
difference between these curves.

####  III.5.2 Dependence on the kinematic viscosity

In Fig. 4(e), we present a log-log plot of $\tau_{1}/\tau_{g}$ versus $Sc$,
over three orders of magnitude and at a fixed value of $Bo$; clearly,
$\tau_{1}/\tau_{g}$ increases with $Sc$, with low-$Sc$ and high-$Sc$
asymptotes, which are consistent with $Sc^{\frac{1}{3}}$ and
$Sc^{\frac{2}{3}}$.

In Fig. 4(f), we plot $h_{min}$ versus $t/\tau_{g}$, at $Bo=1.85\times
10^{-4}$ and $A=0.01$ and $g=99$, for $Sc=1.75$ (blue line with triangles),
$Sc=66$ (red line with diamonds), and $Sc=245$ (yellow line with crosses). We
note that $h_{min}$ falls rapidly, for $Sc=1.75$, but falls slowly, when
$Sc=5.75$. This observation agrees with experiments Scheid  _et al._ (2012).

###  III.6 Three-dimensional antibubbles

We now present illustrative results for 3D antibubbles evolving under gravity
from our CHNS and VoF runs Q1-Q7, in Table 3, and APR1-APR14, in Table 4,
respectively. In Fig. 5 (a) we give, for comparison, an image of an antibubble
from an experiment (courtesy of C. Kalelkar from Ref. Kalelkar (2017)). In
Figs. 5 (b) and (c) we show psuedocolor plots of two-dimensional sections of
$\phi$ (with $y=\pi$ in our 3D CHNS run R2) at (b) $t=0.1\tau_{g}$, and (c)
$t=0.3\tau_{g}$. These plots show clearly the how the South pole of the
antibubble becomes thinner, with the passage of time, while a dome develops at
its North pole. We also find that, at least for the parameters we use, the
initially spherical antibubble remains axisymmetric. Therefore, we design our
3D VoF simulations to be axisymmetric. In Figs. 5 (d)-(e) we show the rupture
of an antibubble in our 3D axisymmetric APR7 VoF run at (d)$t=0$, (e)
$t=0.1\tau_{g}$, and (f) $t=0.3\tau_{g}$.

As we did in our 2D CHNS studies, we identify the antibubble-breakup time
$\tau_{1}$ as the time when $dE/dt$ reaches its maximum value. In Fig. 5(f) we
plot $\tau_{1}/\tau_{g}$ versus $Bo$ for our 3D CHNS runs Q1-Q6; we find that
this plot is qualitatively similar to its 2D CHNS counterpart in Fig. 4(b). We
recall that the plot in Fig. 3(b) also shows that rim velocities, which we
obtain from 3D axisymmetric VoF runs, are in excellent agreement with the
theoretical prediction Sob’yanin (2015).

Figure 5: (a) Image from an antibubble experiment Kalelkar (2017); psuedocolor
plot of order parameter showing a section at $y=\pi$ from 3D simulations of
antibubbles from run R2 at (b) $t=0.1\tau_{g}$ and (c) $t=0.3\tau_{g}$;
Pseudocolor plots of $\phi$ from our 3D axisymmetric VoF run APR7 run at
(d)$t=0$, (e) $t=0.1\tau_{g}$ (f), and $t=0.3\tau_{g}$ showing majority (blue)
and minority (red) phases; (g) plot of $\tau_{1}/\tau_{g}$ versus $Bo$ from
our 3D CHNS runs Q1-Q6.

##  IV Conclusions

An experimental study Dorbolo  _et al._ (2006) has the title Vita brevis of
antibubbles. We have shown how to use the CHNS system for studying the
spatiotemporal evolution of this short life of antibubbles, in both 2D and 3D.
Our DNSs of the CHNS system allow us to study, numerically and theoretically,
the collapse or breakup of an antibubble because of the gravity-induced
thinning at its South pole. This occurs in experiments via the drainage of air
from the lower end of the antibubble. In some experiments, where antibubbles
are stabilised, say by the addition of surfactants, antibubble rupture is
precipitated by piercing its shell; we have studied such rupture via the VoF
method. By considering Fourier-space spectra, we have shown that the breakup
of an antibubble leads to turbulence, insofar as these spectra have
significant weight over several orders of magnitude of wavenumbers. Our DNSs
have allowed us to study the dependence of the antibubble lifetime on the
surface tension, which is related inversely to the Bond number $Bo$, the
kinematic viscosity, which is related directly to the Schmidt number $Sc$, and
the ratio $R_{0}/h_{0}$ of the initial antibubble radius $R_{0}$ and its
thickness $h_{0}$. We have also shown how the antibubble-rim-retraction
velocity depends on the surface tension, in both CHNS and VoF DNSs.

The dependence of the antibubble lifetime on the surface tension has been
studied earlier Scheid  _et al._ (2012), through experiments and numerical
simulation of the balance equations obtained from lubrication theory. A recent
study Yang  _et al._ (2020) uses the Allen-Cahn-Navier-Stokes equations, with
order-parameter conservation enforced numerically via a Lagrange multiplier.
Our work extends significantly theoretical and DNS studies of antibubbles, by
using the CHNS system, in which order-parameter conservation is built in
manifestly, and the VoF method that can be used fruitfully if interfaces are
very thin and non-Boussinesq effects are present. Our results agree with
earlier results, where they exist. The number of studies on antibubbles is
limited partly because great care has to be exercised to stabilize antibubbles
in terrestrial (as opposed to zero-gravity) experiments. Often a surfactant
has to be introduced into the high-density liquid phase for such
stabilization. We hope our detailed study of the spatiotemporal evolution of
antibubbles, from their initiation to their rupture, and of their effect on
the background fluid, will lead to more experimental studies of the properties
of these ephemeral, but beautiful, inverted bubbles.

###### Acknowledgements.

We thank Chirag Kalelkar for introducing us to experimental studies of
antibubbles and for the image in Fig. 5 (a) from his experiments Kalelkar
(2017). We thank Nadia Bihari Padhan and Sumantra Sarkar for discussions. NP
thanks Los Alamos National Laboratory for support; PP and RR acknowledge
support from intramural funds at TIFR Hyderabad from the Department of Atomic
Energy (DAE), India and DST (India) Project No. ECR/2018/001135; RP thanks
CSIR and DST (India) for support and SERC (IISc) for computational resources.

## References

  * Hughes and Hughes (1932) W. Hughes and A.R. Hughes, “Liquid drops on the same liquid surface,” Nature 129, 59 (1932).
  * Kotopoulis  _et al._ (2014) S. Kotopoulis, A. Delalande, M. Popa, V. Mamaeva, G. Dimcevski, O.H. Gilja, M. Postema, B. Tore Gjertsen, and E. McCormack, “Sonoporation-enhanced chemotherapy significantly reduces primary tumour burden in an orthotopic pancreatic cancer xenograft,” Molecular imaging and biology 16, 53–62 (2014).
  * Johansen  _et al._ (2015a) K. Johansen, S. Kotopoulis, A.T. Poortinga, and M. Postema, “Nonlinear echoes from encapsulated antibubbles,” Physics Procedia 70, 1079–1082 (2015a).
  * Kotopoulisa  _et al._ (2015) S Kotopoulisa, K Johansenb, OH Giljaa, AT Poortingad, and M Postemab, “Acoustically active antibubbles,” Dynamics 2, 9 (2015).
  * Postema and Gilja (2007) Michiel Postema and Odd H Gilja, “Ultrasound-directed drug delivery,” Current pharmaceutical biotechnology 8, 355–361 (2007).
  * Johansen  _et al._ (2015b) Kristoffer Johansen, Spiros Kotopoulis, and Michiel Postema, “Ultrasonically driven antibubbles encapsulated by newtonian fluids for active leakage detection,” (2015).
  * Dorbolo  _et al._ (2006) Stéphane Dorbolo, Nicolas Vandewalle, Etienne Reyssat, and David Quéré, “Vita brevis of antibubbles,” Europhysics News 37, 24–25 (2006).
  * Kim and Vogel (2006) P Geon Kim and Jerusha Vogel, “Antibubbles: factors that affect their stability,” Colloids and Surfaces A: Physicochemical and Engineering Aspects 289, 237–244 (2006).
  * Dorbolo  _et al._ (2010) Stéphane Dorbolo, Denis Terwagne, R Delhalle, J Dujardin, N Huet, N Vandewalle, and N Denkov, “Antibubble lifetime: influence of the bulk viscosity and of the surface modulus of the mixture,” Colloids and Surfaces A: Physicochemical and Engineering Aspects 365, 43–45 (2010).
  * Dorbolo  _et al._ (2003) S. Dorbolo, H. Caps, and N. Vandewalle, “Fluid instabilities in the birth and death of antibubbles,” New J. of Phys. 5, 161 (2003).
  * Bai  _et al._ (2016) Lixin Bai, Weilin Xu, Pengfei Wu, Weijun Lin, Chao Li, and Delong Xu, “Formation of antibubbles and multilayer antibubbles,” Colloids and Surfaces A: Physicochemical and Engineering Aspects 509, 334–340 (2016).
  * Kalelkar (2017) Chirag Kalelkar, “The inveterate tinkerer,” Resonance 22, 955–960 (2017).
  * Chatzigiannakis  _et al._ (2021) Emmanouil Chatzigiannakis, Nick Jaensson, and Jan Vermant, “Thin liquid films: where hydrodynamics, capillarity, surface stresses, and intermolecular forces meet,” Current Opinion in Colloid & Interface Science , 101441 (2021).
  * Chatzigiannakis and Vermant (2020) Emmanouil Chatzigiannakis and Jan Vermant, “Breakup of thin liquid films: From stochastic to deterministic,” Physical Review Letters 125, 158001 (2020).
  * Song  _et al._ (2020) Youngsup Song, Lenan Zhang, and Evelyn N Wang, “Criteria for antibubble formation from drop pairs impinging on a free surface,” Physical Review Fluids 5, 123601 (2020).
  * Will  _et al._ (2021) Jelle B Will, Varghese Mathai, Sander G Huisman, Detlef Lohse, Chao Sun, and Dominik Krug, “Kinematics and dynamics of freely rising spheroids at high reynolds numbers,” Journal of Fluid Mechanics 912 (2021).
  * Li  _et al._ (2021) Yanshen Li, Christian Diddens, Andrea Prosperetti, and Detlef Lohse, “Marangoni instability of a drop in a stably stratified liquid,” Physical Review Letters 126, 124502 (2021).
  * Soto  _et al._ (2020a) Álvaro Moreno Soto, Pablo Peñas, Guillaume Lajoinie, Detlef Lohse, and Devaraj van der Meer, “Ultrasound-enhanced mass transfer during single-bubble diffusive growth,” Physical Review Fluids 5, 063605 (2020a).
  * Soto  _et al._ (2020b) Álvaro Moreno Soto, Detlef Lohse, and Devaraj van der Meer, “Diffusive growth of successive bubbles in confinement,” Journal of fluid mechanics 882 (2020b).
  * Prakash  _et al._ (2016) Vivek N Prakash, J Martínez Mercado, Leen van Wijngaarden, Ernesto Mancilla, Yoshiyuki Tagawa, Detlef Lohse, and Chao Sun, “Energy spectra in turbulent bubbly flows,” Journal of fluid mechanics 791, 174–190 (2016).
  * Scheid  _et al._ (2012) B. Scheid, S. Dorbolo, L.R. Arriaga, and E. Rio, “Antibubble dynamics: The drainage of an air film with viscous interfaces,” Phys. Rev. Lett. 109, 264502 (2012).
  * Zou  _et al._ (2013) J. Zou, C. Ji, B.G. Yuan, X.D. Ruan, and X. Fu, “Collapse of an antibubble,” Phys. Rev. E 87, 061002 (2013).
  * Sob’yanin (2015) D.N. Sob’yanin, “Theory of the antibubble collapse,” Phys. Rev. Lett. 114, 104501 (2015).
  * Yang  _et al._ (2020) Junxiang Yang, Yibao Li, Darae Jeong, and Junseok Kim, “Mathematical modeling and simulation of antibubble dynamics,” NUMERICAL MATHEMATICS-THEORY METHODS AND APPLICATIONS 13, 81–98 (2020).
  * Vitry  _et al._ (2019) Youen Vitry, Stéphane Dorbolo, Jan Vermant, and Benoit Scheid, “Controlling the lifetime of antibubbles,” Advances in colloid and interface science 270, 73–86 (2019).
  * Boffetta and Ecke (2012) Guido Boffetta and Robert E Ecke, “Two-dimensional turbulence,” Annual Review of Fluid Mechanics 44, 427–451 (2012).
  * Perlekar  _et al._ (2014) Prasad Perlekar, Roberto Benzi, Herman JH Clercx, David R Nelson, and Federico Toschi, “Spinodal decomposition in homogeneous and isotropic turbulence,” Physical review letters 112, 014502 (2014).
  * Scarbolo  _et al._ (2013) Luca Scarbolo, Dafne Molin, Prasad Perlekar, Mauro Sbragaglia, Alfredo Soldati, and Federico Toschi, “Unified framework for a side-by-side comparison of different multicomponent algorithms: Lattice boltzmann vs. phase field model,” Journal of Computational Physics 234, 263–279 (2013).
  * Pal  _et al._ (2016) N. Pal, P. Perlekar, A. Gupta, and R. Pandit, “Binary-fluid turbulence: Signatures of multifractal droplet dynamics and dissipation reduction,” Phys. Rev. E 93, 063115 (2016).
  * Perlekar  _et al._ (2017) Prasad Perlekar, Nairita Pal, and Rahul Pandit, “Two-dimensional turbulence in symmetric binary-fluid mixtures: Coarsening arrest by the inverse cascade,” Scientific Reports 7, 44589 (2017).
  * Gibbon  _et al._ (2016) John D Gibbon, Nairita Pal, Anupam Gupta, and Rahul Pandit, “Regularity criterion for solutions of the three-dimensional cahn-hilliard-navier-stokes equations and associated computations,” Physical Review E 94, 063103 (2016).
  * Gibbon  _et al._ (2018) J.D. Gibbon, A. Gupta, N. Pal, and R. Pandit, “The role of bkm-type theorems in 3d euler, navier–stokes and cahn–hilliard–navier–stokes analysis,” Physica D 376-377, 60–68 (2018).
  * Fan  _et al._ (2016) Xiang Fan, PH Diamond, L Chacón, and Hui Li, “Cascades and spectra of a turbulent spinodal decomposition in two-dimensional symmetric binary liquid mixtures,” Physical Review Fluids 1, 054403 (2016).
  * Fan  _et al._ (2017) Xiang Fan, PH Diamond, and Luis Chacón, “Formation and evolution of target patterns in cahn-hilliard flows,” Physical Review E 96, 041101 (2017).
  * Fan  _et al._ (2018) Xiang Fan, PH Diamond, and L Chacón, “Chns: A case study of turbulence in elastic media,” Physics of Plasmas 25, 055702 (2018).
  * Celani  _et al._ (2009) A. Celani, A. Mazzino, P. Muratore-Ginanneschi, and L. Vozella, “Phase-field model for the rayleigh–taylor instability of immiscible fluids,” J. Fluid Mech. 622, 115–134 (2009).
  * Scarbolo and Soldati (2013) L. Scarbolo and A. Soldati, “Turbulence modulation across the interface of a large deformable drop,” J. Turbulence 14, 27–43 (2013).
  * Yue  _et al._ (2004) P. Yue, J.J. Feng, C. Liu, and J. Shen, “A diffuse-interface method for simulating two-phase flows of complex fluids,” J. Fluid Mech. 515, 293–317 (2004).
  * Scarbolo  _et al._ (2011) L. Scarbolo, D. Molin, and A. Soldati, “Phase field model simulation of droplet deformation and breakup in wall bounded turbulence,” in  _APS Division of Fluid Dynamics Meeting Abstracts_ (2011).
  * Pandit  _et al._ (2009) Rahul Pandit, Prasad Perlekar, and Samriddhi Sankar Ray, “Statistical properties of turbulence: an overview,” Pramana 73, 157 (2009).
  * (41) “Attached supplementary material.” .
  * Canuto  _et al._ (2012) C. Canuto, M.Y. Hussaini, A. Quarteroni, A. Thomas Jr, _et al._ ,  _Spectral methods in fluid dynamics_ (Springer Science & Business Media, 2012).
  * Cox and Matthews (2002) S.M. Cox and P.C. Matthews, “Exponential time differencing for stiff systems,” J. Comp. Phys. 176, 430–455 (2002).
  * (44) http://www.nvidia.com/object/cuda$\\_$home$\\_$new.html., .
  * (45) http://basilisk.fr., .
  * Popinet (2018) Stéphane Popinet, “Numerical models of surface tension,” Annu. Rev. Fluid Mech. 50, 1–28 (2018).
  * Bell  _et al._ (1989) John B Bell, Phillip Colella, and Harland M Glaz, “A second-order projection method for the incompressible navier-stokes equations,” Journal of Computational Physics 85, 257–283 (1989).