2 Formulation of the problem

We study a two-dimensional incompressible fluid that fills a bounded domain ${\mathcal{D}}$ with mass density $\unicode[STIX]{x1D70C}=1$ . The boundary of the fluid domain, $\unicode[STIX]{x2202}{\mathcal{D}}(t),$ is a free surface given in the implicit form $F(x,y,t)=0$ . The fluid flow is potential and the velocity field is given by $\boldsymbol{v}(x,y,t)=\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}(x,y,t)$ and

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}=0. & & \displaystyle\end{eqnarray}$$

The non-stationary Bernoulli equation governs the evolution of the velocity potential, in particular, at $\unicode[STIX]{x2202}{\mathcal{D}}$ ,

(2.2) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D711}_{t}+{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})^{2}+p)|_{F(x,y,t)=0}=0, & & \displaystyle\end{eqnarray}$$

where $p=p(x,y,t)$ is the pressure, and the Bernoulli constant was absorbed in $\unicode[STIX]{x1D711}_{t}$ . We neglect the atmospheric pressure at the free surface and have $p=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}$ , where $\unicode[STIX]{x1D70E}$ is the surface tension coefficient and $\unicode[STIX]{x1D705}$ is the magnitude of the local curvature.

In general, a fluid particle on the free surface moves in both the tangential and the normal direction; however, it is only the motion in the normal direction that changes the shape of the fluid boundary. This motion is captured by the kinematic boundary condition

(2.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}t}+|\unicode[STIX]{x1D735}F|(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711},\boldsymbol{n})|_{F(x,y,t)}=0, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{n}=\boldsymbol{n}(x,y,t)$ is the outward unit normal to $\unicode[STIX]{x2202}{\mathcal{D}}$ at the point $(x,y)$ .

The total energy, ${\mathcal{H}},$ associated with the fluid flow is given by the sum of the kinetic energy, $T,$ and the potential energy, ${\mathcal{U}}$ :

(2.4a,b ) $$\begin{eqnarray}\displaystyle T=\frac{1}{2}\iint _{{\mathcal{D}}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D711})^{2}\,\text{d}x\,\text{d}y,\quad {\mathcal{U}}=\unicode[STIX]{x1D70E}\int _{\unicode[STIX]{x2202}D}\,\text{d}l. & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D70E}$ is the coefficient of surface tension and $\text{d}l$ is the elementary arclength along $\unicode[STIX]{x2202}{\mathcal{D}}$ .

Figure 1. Illustration of the conformal map: the fluid domain $z=x+\text{i}y$ is enclosed by the free surface $\unicode[STIX]{x2202}{\mathcal{D}}$ (a), and the conformal domain $w=u+\text{i}v$ is the periodic strip in the lower complex half-plane (b).

3 Conformal variables

Let $z(w)=x(w)+\text{i}y(w)$ be a conformal map to the fluid domain $(x,y)\in {\mathcal{D}}$ from a periodic strip $w=u+\text{i}v,-\unicode[STIX]{x03C0}\leqslant k_{0}u<\unicode[STIX]{x03C0},v\leqslant 0$ and $k_{0}$ is the base wavenumber for the parameterization of the circle. The illustration of the mapping is given in figure 1. The free surface $\unicode[STIX]{x2202}{\mathcal{D}}$ is mapped from $v=0$ and satisfies

(3.1) $$\begin{eqnarray}\displaystyle F(x(u,t),y(u,t),t)=0. & & \displaystyle\end{eqnarray}$$

The kinematic condition can be revealed from the observation of the time derivative of the implicit function $F$ , i.e.

(3.2) $$\begin{eqnarray}\displaystyle \frac{\text{d}F}{\text{d}t}=F_{t}+F_{x}x_{t}+F_{y}y_{t}=0, & & \displaystyle\end{eqnarray}$$

where the subscript denotes a partial derivative. After a change of coordinates from the $(x,y)$ -plane to the $(u,v)$ -plane, and applying the chain rule it becomes

(3.3) $$\begin{eqnarray}\displaystyle F_{t}+\frac{x_{t}x_{u}+y_{t}y_{u}}{|z_{u}|^{2}}F_{u}+\frac{x_{t}y_{u}-y_{t}x_{u}}{|z_{u}|^{2}}F_{v}=0, & & \displaystyle\end{eqnarray}$$

and the coefficient at $F_{v}$ gives the rate of change of the free surface $\unicode[STIX]{x2202}{\mathcal{D}}$ in the vertical direction in the $w$ -plane, which translates to the normal direction in the $(x,y)$ -plane. In conformal variables (2.3) can be written by noting that the unit normal $\boldsymbol{n}$ is given by

(3.4) $$\begin{eqnarray}\displaystyle \boldsymbol{n}=\frac{1}{|z_{u}|}(-y_{u},x_{u})^{T}, & & \displaystyle\end{eqnarray}$$

where we exploit the fact that $z(w)$ is subject to Cauchy–Riemann (CR) relations. A change of variables in (2.3) thus reveals that

(3.5) $$\begin{eqnarray}\displaystyle F_{t}+\frac{\unicode[STIX]{x1D713}_{v}}{|z_{u}|^{2}}F_{v}=0, & & \displaystyle\end{eqnarray}$$

where we have introduced a restriction of the velocity potential to the free surface:

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(u,t)=\unicode[STIX]{x1D711}(x,y,t)|_{F(x,y,t)=0}. & & \displaystyle\end{eqnarray}$$

By matching the coefficient at the partial derivative $F_{v}$ in (3.3) and (3.5), we discover the kinematic condition in the conformal domain,

(3.7) $$\begin{eqnarray}\displaystyle \frac{x_{t}y_{u}-y_{t}x_{u}}{|z_{u}|^{2}}=\frac{\unicode[STIX]{x1D713}_{v}}{|z_{u}|^{2}}, & & \displaystyle\end{eqnarray}$$

which ensures that the line $v=0$ maps to the free surface for all time. Having a derivative with respect to $v$ is a nuisance that can be alleviated by making use of the stream function, $\unicode[STIX]{x1D703}$ , and using CR relations together with Titchmarsh’s theorem (Titchmarsh Reference Titchmarsh1986) to find that $\unicode[STIX]{x1D713}_{v}=-{\hat{H}}\unicode[STIX]{x1D713}_{u}$ , where ${\hat{H}}$ denotes the Hilbert transform,

(3.8) $$\begin{eqnarray}\displaystyle {\hat{H}}f(u)=v.p.\int _{-\infty }^{\infty }\frac{f(u^{\prime })\,\text{d}u^{\prime }}{u^{\prime }-u}, & & \displaystyle\end{eqnarray}$$

where $v.p.$ denotes a Cauchy principal value integral.

The Bernoulli equation (2.2) determines the evolution of the velocity potential at the free surface. It is formulated in conformal variables by means of elementary calculus. The force of surface tension is proportional to the local curvature of $\unicode[STIX]{x2202}{\mathcal{D}}$ :

(3.9) $$\begin{eqnarray}\displaystyle p=-\unicode[STIX]{x1D70E}\frac{x_{u}y_{uu}-y_{u}x_{uu}}{|z_{u}|^{3}}. & & \displaystyle\end{eqnarray}$$

Under the conformal change of variables the surface potential $\unicode[STIX]{x1D713}(u,t)$ becomes a composite function,

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(u,t)=\unicode[STIX]{x1D711}(x(u,t),y(u,t),t), & & \displaystyle\end{eqnarray}$$

the time derivative of which together with (2.2), (3.7) and (3.9) reveal that

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{t}=-\frac{\unicode[STIX]{x1D713}_{u}^{2}-({\hat{H}}\unicode[STIX]{x1D713}_{u})^{2}}{2|z_{u}|^{2}}+\unicode[STIX]{x1D713}_{u}{\hat{H}}\left[\frac{{\hat{H}}\unicode[STIX]{x1D713}_{u}}{|z_{u}|^{2}}\right]-p, & & \displaystyle\end{eqnarray}$$

the dynamic boundary condition in the conformal domain.

4 The complex equations

In order to reveal the analytic structure of the problem at hand, it is convenient to introduce the complex potential $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D713}+\text{i}\unicode[STIX]{x1D703}$ . By the Titchmarsh theorem, the complex potential can also be written in the form

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D713}+\text{i}{\hat{H}}\unicode[STIX]{x1D713}=2\hat{P}\unicode[STIX]{x1D713}, & & \displaystyle\end{eqnarray}$$

where $2\hat{P}=1+\text{i}{\hat{H}}$ and $\hat{P}$ is the projection operator. Both the complex functions $\unicode[STIX]{x1D6F7}$ and $z$ are analytic in the periodic strip in the lower complex half-plane. Unlike the problem of travelling waves on a fluid of finite or infinite depth, in our formulation the free surface $\unicode[STIX]{x2202}D$ is a closed curve in the $(x,y)$ -plane and, hence, the conformal map $z(u+\text{i}v)$ must be a periodic function of the variable $u$ , rather than an identity transformation plus a periodic correction. Therefore, $z(w)$ (as well as $\unicode[STIX]{x1D6F7}(w)$ ) is expanded as a Fourier series, and the analyticity in the periodic strip requires that only the non-positive Fourier coefficients are non-zero, i.e.

(4.2) $$\begin{eqnarray}\displaystyle z(w)=\hat{z}_{0}+\mathop{\sum }_{k_{0}\mathbb{N}}\hat{z}_{k}\text{e}^{-\text{i}kw}, & & \displaystyle\end{eqnarray}$$

where $\hat{z}_{k}$ denotes the Fourier coefficients of $z$ . As evident from this expansion, as $w\rightarrow -\text{i}\infty$ the derivative of the conformal map $z_{w}\rightarrow 0$ and $z(w\rightarrow -\text{i}\infty )=z_{0}$ . In other words, we define a new analytic function $\unicode[STIX]{x1D70C}(w)$ that satisfies

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}z_{w}=\text{e}^{-\text{i}k_{0}w} & & \displaystyle\end{eqnarray}$$

at every point in the strip, and $\unicode[STIX]{x1D70C}(-\text{i}\infty )=(-\text{i}k_{0}\hat{z}_{1})^{-1}\neq 0$ . We also introduce a new zero-mean, analytic function, $\unicode[STIX]{x1D708}=\text{i}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6F7}_{u}$ . When the kinematic condition is written in the complex form and multiplied by $|\unicode[STIX]{x1D70C}|^{2}$ , the choice of $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D708}$ becomes transparent, i.e.

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle z_{t}\bar{z}_{u}\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D70C}}-\bar{z}_{t}z_{u}\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D70C}}=\bar{\unicode[STIX]{x1D6F7}}_{u}\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D70C}}-\unicode[STIX]{x1D6F7}_{u}\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}z_{t}\text{e}^{\text{i}k_{0}w}-\bar{\unicode[STIX]{x1D70C}}\bar{z}_{t}\text{e}^{-\text{i}k_{0}w}=\text{i}(\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D708}}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D708}), & \displaystyle\end{eqnarray}$$

where $z_{t}(w\rightarrow -\infty )\rightarrow 0$ and the left-hand side is the difference of two complex analytic functions. We apply the projection operator $\hat{P}$ to have

(4.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}z_{t}\text{e}^{\text{i}k_{0}u}=\text{i}\hat{P}[\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D708}}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D708}]. & & \displaystyle\end{eqnarray}$$

After elementary calculation, we conclude that the new analytic function $\unicode[STIX]{x1D70C}$ satisfies the pseudo-differential equation

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{t}=\text{i}(U\unicode[STIX]{x1D70C}_{u}-U_{u}\unicode[STIX]{x1D70C})-k_{0}\unicode[STIX]{x1D70C}U, & & \displaystyle\end{eqnarray}$$

where $U$ , given by

(4.8) $$\begin{eqnarray}\displaystyle U=\hat{P}(\unicode[STIX]{x1D70C}\bar{\unicode[STIX]{x1D708}}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D708}), & & \displaystyle\end{eqnarray}$$

is the complex transport velocity. The equation for the complex potential is found by applying $2\hat{P}$ to (3.11):

(4.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{t}=\text{i}U\unicode[STIX]{x1D6F7}_{u}-B. & & \displaystyle\end{eqnarray}$$

Here $B$ is given by

(4.10) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle B=\hat{P}\left[\frac{|\unicode[STIX]{x1D6F7}_{u}|^{2}}{|z_{u}|^{2}}-2\unicode[STIX]{x1D70E}\frac{x_{u}y_{uu}-y_{u}x_{uu}}{|z_{u}|^{3}}\right],\\ \displaystyle B=\hat{P}[|\unicode[STIX]{x1D708}|^{2}+2\unicode[STIX]{x1D70E}k_{0}|\unicode[STIX]{x1D701}|^{2}+2\text{i}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D701}\bar{\unicode[STIX]{x1D701}}_{u}-\bar{\unicode[STIX]{x1D701}}\unicode[STIX]{x1D701}_{u})],\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D701}^{2}=\unicode[STIX]{x1D70C}$ . We omit the trivial details of the calculation and skip to the result, which is the equation satisfied by $\unicode[STIX]{x1D708}$ :

(4.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}_{t}=\text{i}(U\unicode[STIX]{x1D708}_{u}-B_{u}\unicode[STIX]{x1D70C})-k_{0}U\unicode[STIX]{x1D708}. & & \displaystyle\end{eqnarray}$$

The newly discovered equation (4.11) is analogous to that originally discovered in Dyachenko (Reference Dyachenko2001) for the infinite fluid domain, and when complemented with (4.7) forms a closed system suitable for numerical simulation.

5 The choice of the reference frame

Equations (4.7) and (4.11) are formulated for the derivative of a conformal map and, henceforth, describe only the motion of the fluid relative to the point of the fluid, $\hat{z}_{0}$ . The motion of this point, namely, $\hat{z}_{0}(t)=z(-\text{i}\infty )$ is not captured in (4.7) and (4.11) and has to be recovered elsewhere. This missing puzzle piece comes from the momentum conservation that implies the centre of mass of the fluid is an inertial reference frame, and we chose it to be placed at the origin, $z=0$ :

(5.1) $$\begin{eqnarray}\displaystyle \iint _{{\mathcal{D}}}z\,\text{d}x\,\text{d}y=\iint z|z_{u}|^{2}\,\text{d}u\,\text{d}v=0. & & \displaystyle\end{eqnarray}$$

With the total mass of the fluid being a constant of motion given by

(5.2) $$\begin{eqnarray}\displaystyle m=\iint _{{\mathcal{D}}}\,\text{d}x\,\text{d}y, & & \displaystyle\end{eqnarray}$$

(5.1) can be written as

(5.3) $$\begin{eqnarray}\displaystyle m\hat{z}_{0}=-\iint (z-\hat{z}_{0})|z_{u}|^{2}\,\text{d}u\,\text{d}v & & \displaystyle\end{eqnarray}$$

and is used to recover the zero Fourier mode of the conformal map, quite similar to the zero mean condition that is often imposed in the problem with infinite fluid domain. Together with the system (4.7), (4.11) and equation (5.3), we can fully describe the motion of the boundary of the fluid $\unicode[STIX]{x2202}{\mathcal{D}}$ .

Figure 2. (a) The size of the semiminor axis of the ellipse according to the theoretical result from Longuet-Higgins (Reference Longuet-Higgins1972) (dashed) and from the direct numerics (solid). The snapshots of the surface shape corresponding to the points labelled by $A$ , $B$ and $C$ are illustrated on the right panel. (b) The shape of the fluid boundary at time $t=0$ , $t=0.24$ (A), $t=0.48$ (B) and $t=1.00$ (C). The shape of the surface is found from numerical simulation of (4.7) and (4.11) with initial conditions in accordance with the Dirichlet ellipse (6.5), except for case $C$ where the aspect ratio of the ellipse is taken directly from the theory.

Equations (4.7) and (4.11) for the analytic functions $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D708}$ are equivalent to the equations for $R=1/z_{u}$ and $V=\text{i}\unicode[STIX]{x1D6F7}_{u}R$ , i.e. 

(5.4a,b ) $$\begin{eqnarray}\displaystyle R_{t}=\text{i}(UR_{u}-U_{u}R),\quad V_{t}=\text{i}(UV_{u}-B_{u}R), & & \displaystyle\end{eqnarray}$$

with

(5.5) $$\begin{eqnarray}\displaystyle & \displaystyle U=\hat{P}[V\bar{R}+\bar{V}R], & \displaystyle\end{eqnarray}$$

(5.6) $$\begin{eqnarray}\displaystyle & \displaystyle B=\hat{P}[|V^{2}|+2\text{i}\unicode[STIX]{x1D70E}(Q\bar{Q}_{u}-\bar{Q}Q_{u})], & \displaystyle\end{eqnarray}$$

where $Q^{2}=R$ .

There is, however, an important caveat, namely the function $R$ is no longer analytic at $w\rightarrow -\text{i}\infty$ , but instead it contains a term in its Fourier series expansion with a positive wavenumber, $k_{0}$ :

(5.7) $$\begin{eqnarray}\displaystyle R(w)=\text{e}^{\text{i}k_{0}w}\left(\hat{\unicode[STIX]{x1D70C}}_{0}+\mathop{\sum }_{k_{0}\mathbb{N}}\hat{\unicode[STIX]{x1D70C}}_{k}\text{e}^{-\text{i}kw}\right)=\unicode[STIX]{x1D70C}\text{e}^{\text{i}k_{0}w}. & & \displaystyle\end{eqnarray}$$

6 Numerical experiments

In the remainder of the paper we demonstrate the simulations of (4.7) and (4.11) and compare the results with available exact solutions. The simulations are performed on a uniform grid in the $u$ -variable using the Runge–Kutta fourth-order timestepping scheme. The spatial derivative, $\unicode[STIX]{x2202}_{u}$ , and projection operator, $\hat{P}$ , are applied as Fourier multipliers to the coefficients of Fourier series of the respective functions.

In the first simulation we compare the solutions of (4.7) and (4.11) with the Dirichlet ellipse (see Longuet-Higgins (Reference Longuet-Higgins1972) for details). The complex potential, $\unicode[STIX]{x1D6F7}$ , is a quadratic function of the conformal map $z$ for Dirichlet solutions:

(6.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}={\textstyle \frac{1}{2}}Az^{2}+\int f\,\text{d}t. & & \displaystyle\end{eqnarray}$$

Here $A=A(t)$ and $\int f\,\text{d}t$ is the Bernoulli constant. The surface shape $\unicode[STIX]{x2202}{\mathcal{D}}$ is given by

(6.2) $$\begin{eqnarray}\displaystyle \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1, & & \displaystyle\end{eqnarray}$$

where $a=a(t)$ and $b=b(t)$ are determined from $A$ as follows:

(6.3) $$\begin{eqnarray}\displaystyle a^{2}=\frac{1}{b^{2}}=\frac{A^{2}}{1-\sqrt{1-A^{2}}} & & \displaystyle\end{eqnarray}$$

and $A$ satisfies an ordinary differential equation (ODE)

(6.4) $$\begin{eqnarray}\displaystyle \frac{\text{d}A}{\text{d}t}=A^{2}\sqrt{1-A^{4}}. & & \displaystyle\end{eqnarray}$$

This ODE is solved numerically and is compared with the solution of (4.7) and (4.11), and agreement is demonstrated in figure 2(a). We measure the size of the semiminor axis $b(t)$ as obtained from both simulations and plot the result in figure 2(b). The initial data for the simulation is given by

(6.5a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}(t=0)=1,\quad \unicode[STIX]{x1D708}(t=0)=-\sqrt{2}\text{e}^{-\text{i}u}, & & \displaystyle\end{eqnarray}$$

and in the original variables yields

(6.6a,b ) $$\begin{eqnarray}\displaystyle z(t=0)=\text{e}^{-\text{i}w},\quad \unicode[STIX]{x1D6F7}(t=0)=\frac{\sqrt{2}}{2}z^{2}. & & \displaystyle\end{eqnarray}$$

Figure 3. Numerical simulation of the motion of a droplet with surface tension coefficient $\unicode[STIX]{x1D70E}=2$ : (a) The initial shape of the surface (blue) and the shape of the surface at time $t=0.725$ (red) and $t=1.450$ (green) in a simulation with initial kinetic energy $K=0.25\unicode[STIX]{x03C0}$ . The simulation stops when the spectrum of the free surface has become too wide and is not resolved on a uniform grid with $8192$ Fourier modes (see figure 4). (b) The results of simulations with $c^{2}=0.4,1.2$ and $2.0$ (blue, yellow and green, respectively). The oscillatory motion is exhibited for the blue and yellow curves, but for the green curve (kinetic energy $K=0.25\unicode[STIX]{x03C0}$ ) the restoring force from surface tension is insufficient, and the fluid droplet forms a narrowing neck and the droplet breaking is conjectured. We observe that the period of oscillations increases when more kinetic energy is available, and, thus, the restoring force becomes weaker as the droplet extends from its equilibrium shape. Although simulation for $K=0.25\unicode[STIX]{x03C0}$ stops at time $t=1.450$ , we conjecture that a capillary instability will occur and result in subsequent breaking of the droplet.

In the second set of simulations, we solve (4.7) and (4.11) in the presence of surface tension with $\unicode[STIX]{x1D70E}=2$ and take the initial data for the conformal map and the complex potential to be

(6.7a,b ) $$\begin{eqnarray}\displaystyle z(t=0)=\text{e}^{-\text{i}w},\quad \unicode[STIX]{x1D6F7}(t=0)=\frac{c}{2}\text{e}^{-2\text{i}w}, & & \displaystyle\end{eqnarray}$$

where the values of $c^{2}$ are $0.4$ , $1.2$ and $2.0$ . At the initial time the surface shape is a unit disc and the complex potential is a quadratic function of $z$ , and if this simulation was run at $\unicode[STIX]{x1D70E}=0$ , then the result would be a Dirichlet ellipse solution. In contrast to the $\unicode[STIX]{x1D70E}=0$ case, the surface shape exhibits oscillatory motion that is quite similar to standing waves in an infinite fluid, yet surface tension provides the restoring force instead of gravity. In the $c^{2}=2$ case, corresponding to the initial kinetic energy $K=0.25\unicode[STIX]{x03C0}$ , the restoring force is insufficient to lead to droplet oscillations, instead the droplet shape develops a narrowing neck and is reminiscent of fragmentation to smaller droplets (see figure 3).

We measure the accuracy of the simulations by verifying conservation of the fluid mass, $m$ , and the Hamiltonian, ${\mathcal{H}}$ , i.e.

(6.8) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=-\frac{1}{2}\int \unicode[STIX]{x1D713}{\hat{H}}\unicode[STIX]{x1D713}_{u}\,\text{d}u+\unicode[STIX]{x1D70E}\int |z_{u}|\,\text{d}u, & & \displaystyle\end{eqnarray}$$

which are conserved to $13$ digits of precision. In addition, we track the magnitude of Fourier coefficients $|\hat{\unicode[STIX]{x1D701}}_{k}|$ and $|\hat{\unicode[STIX]{x1D708}}_{k}|$ and ensure that it is resolved, illustration of the Fourier spectrum is presented in figure 4.

Figure 4. Numerical simulation of the motion of a droplet with the surface tension coefficient $\unicode[STIX]{x1D70E}=2$ : the Fourier coefficients $|\hat{\unicode[STIX]{x1D701}}_{k}|$ (solid) and $|\hat{\unicode[STIX]{x1D708}}_{k}|$ (dashed) computed at times $t=0$ (blue), $t=0.725$ (red) and $t=1.425$ (green).

We observe that given enough kinetic energy the restoring force provided by surface tension becomes incapable of supporting oscillatory motion (see figure 3(b)), and it becomes energy efficient to break the droplet into two smaller droplets that carry away half the kinetic energy each, because the characteristic volume of the resulting droplets is halved and the perimeter is decreased by $\sqrt{2}$ , the smaller droplets are more stable to further splitting than the initial droplet. However, when the energy is sufficiently large subsequent fragmentation of the smaller droplets may occur.

The actual process of droplet breaking is not tractable with the present techniques of conformal mapping, because at the instant of breaking the topology of the fluid domain changes. Furthermore, as the neck between the two cores extends and its walls close-in, the number of Fourier modes required to resolve the solution grows rapidly, suggesting that singularities in the analytic continuation of $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D708}$ are approaching the real axis $v=0$ .