The onset of filamentation on vorticity interfaces in two-dimensional Euler flows

Abstract

Two-dimensional Euler flows, in the plane or on simple surfaces, possess a material invariant, namely the scalar vorticity normal to the surface. Consequently, flows with piecewise-uniform vorticity remain that way, and moreover evolve in a way which is entirely determined by the instantaneous shapes of the contours (interfaces) separating different regions of vorticity – this is known as ‘contour dynamics’. Unsteady vorticity contours or interfaces often grow in complexity (lengthen and fold), either as a result of vortex interactions (like mergers) or ‘filamentation’. In the latter, wave disturbances riding on a background, equilibrium contour shape appear to inevitably steepen and break, forming filaments, repeatedly– and perhaps endlessly. Here, we revisit the onset of filamentation. Building upon previous work and using a weakly nonlinear expansion to third order in wave amplitude, we derive a universal, parameter-free amplitude equation which applies (with a minor change) both to a straight interface and a circular patch in the plane, as well as circular vortex patches on the surface of a sphere. We show that this equation possesses a local, self-similar form describing the finite-time blow up of the wave slope (in a re-scaled long time proportional to the inverse square of the initial wave amplitude). We present numerical evidence for this self-similar blow-up solution, and for the conjecture that almost all initial conditions lead to finite-time blow up. In the full contour dynamics equations, this corresponds to the onset of filamentation.

1. Introduction

Ideal two-dimensional flows, governed by Euler’s equations, constitute an infinite-dimensional Hamiltonian system (Morrison 1998). Not only are there energy and momentum invariants (symmetry permitting), but each fluid particle conserves its scalar vorticity (the normal component of vorticity for general orientable surfaces, see Saffman 1995; Dritschel & Boatto 2015). Because such flows are incompressible, they evolve purely by vorticity re-arrangement, determined by a velocity field that depends linearly, but non-locally, on the vorticity distribution. Few exact solutions are known, and almost all of them are either steady or in relative equilibrium (see the discussion in Abrashkin & Yakubovich (1984), Aleman & Constantin (2012), Constantin & Martin (2017), Crowdy (2004), Krishnamurthy et al. (2021), Majda & Bertozzi (2002) and Stuart (1967)).

For piecewise-uniform vorticity distributions, the flow evolution depends only on the instantaneous shapes of the contours separating regions of uniform vorticity and on the vorticity jumps across them (Zabusky, Hughes & Roberts 1979; Dritschel 1989). The resulting system of equations is called ‘contour dynamics’. These equations also constitute an infinite-dimensional Hamiltonian system, manifest in practice by contours which deform, elongate and generally grow in complexity, apparently indefinitely (Dritschel 1988c ). Nonetheless, mathematically, contours of a certain regularity class (at least possessing a continuous tangent vector) preserve that regularity class for all time (Chemin 1993; Bertozzi & Constantin 1993).

The seemingly inevitable growth in complexity of these contours, or ‘vorticity interfaces’, was first investigated by Dritschel (1988c ), who studied the behaviour of small disturbances to circular contours (or ‘vortex patches’) and to straight contours, both of which are otherwise in equilibrium. Numerical simulations performed using contour dynamics, including curvature-controlled point redistribution and a regularisation called ‘surgery’ (Dritschel 1988b ), demonstrated that a range of initial disturbances having small wave slope gradually steepened and folded over, a process called ‘filamentation’. Moreover, after the first filament forms, further filaments are generated at a frequency equal to half of the vorticity jump across the interface. For a fascinating historical account going back to Lord Kelvin (Thomson 1880), see Craik (2012).

The present study focuses on the onset of filamentation, namely the wave-steepening stage and the approach to infinite wave slope. This was also studied mathematically in Dritschel (1988c ), who developed a weakly nonlinear theory describing the progressive steepening of an arbitrary disturbance. For this, the disturbance was expressed as a slowly varying amplitude, $\mathcal {A}$ , multiplied by the relatively fast linear oscillation, $\exp (\mathrm {i}\omega t/2)$ , with frequency equal to half of the vorticity jump ( $\omega /2$ ) across the interface. Then, by expanding the equations of contour dynamics to third order in $\mathcal {A}$ , and by requiring no secular growth, a cubically nonlinear evolution equation for $\mathcal {A}$ was derived. This was shown to accurately describe the onset of filamentation by direct comparison with full contour dynamics numerical simulations.

In the present paper, we revisit this equation and demonstrate that, by introducing an appropriately re-scaled slow time, the amplitude equation for $\mathcal {A}$ is parameter free. Hence, it applies equally well to disturbances propagating on circular contours in the plane and on the surface of a sphere. It also applies, with a minor modification (one less term), to disturbances propagating on straight contours in the plane. A companion paper (Constantin, Germain & Dritschel 2024) focuses on the mathematical structure of this equation, its conservation laws and symmetry properties, as well as exact solutions in special cases. Here, we study this equation numerically and, in particular, provide evidence for a finite-time singularity in the wave slope. This singularity appears to have a self-similar form in the diffusive similarity variable $(x-x_c)/(t_c-t)^{1/2}$ , where $x_c$ is the point of maximum wave slope and $t_c$ is the singularity time.

The paper is organised as follows. In § 2 we review the weakly nonlinear theory presented in Dritschel (1988c ), then show that the amplitude equation can be written in a parameter-free way in an appropriately re-scaled time variable. In § 3, we examine numerical simulations starting from a superposition of two low-wavenumber sinusoidal waves. Generically, for this and many other initial conditions, we find progressive wave steepening. Our results suggest a finite-time singularity, which is here analysed by fitting the wave shape to the form proposed in § 3.2. Our conclusions are offered in § 4.

2. Weakly nonlinear theory for vorticity interfaces

We briefly recap the analysis presented in Dritschel (1988c ), hereafter referred to as D88c. In appendix A therein, a cubic-order amplitude equation was derived for small wave-slope disturbances to a circular vortex patch on the surface of a sphere. This in fact also applies to a circular patch in the two-dimensional plane $\mathbb {R}^2$ in a certain limit, as well as to a straight interface in $\mathbb {R}^2$ , as explained below.

2.1. The sphere

In Cartesian coordinates, the equations of contour dynamics on the sphere closely resemble those on the plane. For a single contour $\mathcal {C}$ across which the vorticity jumps by $\omega$ (the difference between the vorticity to the left of the contour and that to the right), each point $\boldsymbol {x}\in \mathcal {C}$ evolves according to

(2.1) \begin{equation} \frac {\textrm {d}\boldsymbol {x}}{\textrm {d} t}=-\frac {\omega }{4\pi }\oint _{\mathcal {C}}\log |\boldsymbol {x}'-\boldsymbol {x}|^2\,\textrm {d}\boldsymbol {x}' ;\end{equation}

(Dritschel 1988a ). In the case of the sphere, the points $\boldsymbol {x}$ and $\boldsymbol {x}'$ are three-dimensional but constrained to have magnitude $r=1$ , without loss of generality. In the case of the plane, $\boldsymbol {x}$ and $\boldsymbol {x}'$ are two-dimensional.

We now consider a circular vortex patch boundary lying at constant latitude $\phi =\phi _0$ , and separating vorticity $\omega _N$ to the north from $\omega _S$ to the south (the geometric centre of the patch lies on the $z$ axis between the north and south poles). Across the boundary or interface, the vorticity jumps by $\omega \equiv \omega _N-\omega _S$ , which here we take to be unity, again without loss of generality. (Note, this is not a model of rotating planetary flows, which have a continuous vorticity variation $\propto \sin \phi$ , but instead a mathematical model for studying general aspects of wave steepening on curved vorticity interfaces.)

It is convenient to consider displacements in the axial position or ‘height’ $z=\sin \phi$ of the vortex patch boundary, in particular displacements of the form

(2.2) \begin{equation} z(\theta ,t)=z_0-r_0^2\rho (\theta ,t) ,\end{equation}

where $\theta$ is the longitude coordinate, $t$ is time, $z_0=\sin \phi _0$ and $r_0=\cos \phi _0$ . Here, $\rho (\theta ,t)$ is the displacement function, considered small compared with unity (this form facilitates taking the planar limit, i.e. $\phi _0\to \pi /2$ ). We follow (A6) in D88c and take

(2.3) \begin{equation} \rho (\theta ,t)=a\eta (\theta ,t), \end{equation}

where $a\ll 1$ and $\eta ={\mathcal {O}}(1)$ . We assume here that the mean value of $\rho$ , and therefore $\eta$ , vanishes. (The mean value plays no role in the disturbance evolution, and corresponds to an ${\mathcal {O}}(a)$ shift of the mean latitude $\phi _0$ .)

Figure 1. Illustration of the linear evolution of disturbances to a vorticity interface, $\eta (\theta ,t)$ in (2.5), for $-\pi \leqslant \theta \leqslant \pi$ (shown as the abscissa). Here, $T=4\pi /\omega$ is the linear wave period and successive times shown are displaced downwards by an equal increment in the ordinate to avoid overlap. Note that the initial disturbance reverses after half a period, then recovers its initial form after one full period. After a quarter of a period, the solution turns from anti-symmetric to symmetric about its centre, and this also reverses after a further half-period.

Expanding (2.1) to ${\mathcal {O}}(a^3)$ , the equation satisfied by $\eta$ takes the form

(2.4) \begin{align} \frac {1}{\omega }&\frac {\partial \eta (\theta ,t)}{\partial t}-\frac {1}{4\pi }\int _0^{2\pi }\frac {\eta (\alpha ,t)\sin (\alpha -\theta )}{1-\cos (\alpha -\theta )}\textrm {d}\alpha =\frac {\partial G}{\partial \theta } \nonumber \\ &G(\theta ,t)=\frac {a}{4}z_0\eta ^2-\frac {a^2}{6}(z_0^2+1)\eta ^3 -\frac {a^2}{24\pi }\int _0^{2\pi }\frac {[\eta (\alpha ,t)-\eta (\theta ,t)]^3}{1-\cos (\alpha -\theta )}\textrm {d}\alpha \,, \end{align}

in a frame of reference rotating with the mean angular velocity $1/2\omega$ of the patch boundary (this is (A7) in D88c, for zero mean $\eta$ ). Without the nonlinear term ( $G=0$ ), the equation is linear in $\eta$ , and all solutions oscillate at the common frequency $1/2\omega$ , independent of the spatial structure of $\eta$ . Specifically, in linear theory, the general solution to (2.4) for $G=0$ is

(2.5) \begin{equation} \eta (\theta ,t)=\mathcal {A}(\theta )\textrm{e}^{{\frac {1}{2}}\mathrm {i}\omega t} + \mathrm {c.c.} \quad \textrm {with}\, \mathcal {A}(\theta )=\sum _{m=1}^{\infty }a_m \textrm{e}^{\mathrm {i} m\theta }, \end{equation}

where the coefficients $a_m$ are arbitrary complex constants, and ‘c.c.’ denotes complex conjugation. Denoting the real and imaginary parts of $\mathcal {A}$ by $\mathcal {A}_r$ and $\mathcal {A}_i$ , respectively, we can write this solution in the form

(2.6) \begin{align} \eta (\theta ,t)=2\left (\mathcal {A}_r(\theta )\cos \left({\frac {1}{2}}\omega t\right)- \mathcal {A}_i(\theta )\sin \left({\frac {1}{2}}\omega t\right)\right )\,. \end{align}

This is illustrated in figure 1 for an initially anti-symmetric disturbance. Note that, by a quarter of the linear wave period, the solution has become symmetric (this is $-2\mathcal {A}_i$ ), then changes back to the initial disturbance except reversed after another quarter period, and so on. The upshot is that the linear evolution is unsteady, evolving on a time scale inversely proportional to the vorticity jump $\omega$ .

For the nonlinear equation (2.4) with $G\neq 0$ , we follow the multiple-time-scale ansatz in (A8) of D88c and look for a solution of the form

(2.7) \begin{equation} \eta (\theta ,t)=\mathcal {A}(\theta ,t)\textrm{e}^{{\frac {1}{2}}\mathrm {i}\omega t} + \mathrm {c.c.} ,\end{equation}

where $\mathcal {A}(\theta ,t)=\mathcal {A}_0(\theta ,\tau )+a\mathcal {A}_1(\theta ,t,\tau )+\cdots$ , and where $\tau =\omega a^2 t$ is the slow time. Omitting the details, the equation for $\mathcal {A}_0$ (hereafter written simply as $\mathcal {A}$ ) is (A11) of D88c, rewritten here as

(2.8) \begin{equation} \frac {\partial \mathcal {A}}{\partial \tau }=\frac 12 \frac {\partial }{\partial \theta } \left (z_0^2 \mathcal {T}_1 + \mathcal {T}_2 - (z_0^2+1)|\mathcal {A}|^2\mathcal {A}\right ) ,\end{equation}

where

(2.9) \begin{equation} \mathcal {T}_1=\mathrm {i}\left (\mathcal {A}\frac {\partial }{\partial \theta }(\mathcal {W}-\bar {\mathcal {W}})-|\mathcal {A}|^2\frac {\partial \mathcal {A}}{\partial \theta }\right ), \end{equation}

and

(2.10) \begin{equation} \mathcal {T}_2=-\frac {1}{4\pi }\int _0^{2\pi } \frac {|\mathcal {A}(\alpha ,\tau )-\mathcal {A}(\theta ,\tau )|^2[\mathcal {A}(\alpha ,\tau )-\mathcal {A}(\theta ,\tau )]} {1-\cos (\alpha -\theta )}\,\textrm {d}{\alpha }\,. \end{equation}

Above, $\mathcal {W}$ is the part of $|\mathcal {A}|^2$ expressible in positive wavenumbers, i.e.

(2.11) \begin{align} \mathcal {W}=\sum _{m=1}^{\infty } w_m \textrm{e}^{\mathrm {i} m\theta } \,, \end{align}

and an overbar denotes complex conjugation (note that $|\mathcal {A}|^2=\mathcal {W}+\bar {\mathcal {W}}+\mathsf {\textit P}$ , where $\mathsf {\textit P}$ is a constant of the motion, see below). By explicit calculation starting with the general Fourier series

(2.12) \begin{align} \mathcal {A}(\theta ,\tau )=\sum _{m=1}^{\infty }a_m(\tau ) \textrm{e}^{\mathrm {i} m\theta } \,, \end{align}

remarkably, it can be shown that $\mathcal {T}_2=\mathcal {T}_1$ (Constantin et al. 2024). In this case, by introducing the re-scaled slow time $\tilde \tau ={({1}/{2})}(z_0^2+1)\tau$ , we have more simply

(2.13) \begin{align} \frac {\partial \mathcal {A}}{\partial \tilde \tau } &= \frac {\partial \mathcal {B}}{\partial \theta } \\ \mathcal {B}(\theta ,\tilde \tau )&=-\frac {1}{4\pi }\int _0^{2\pi } \frac {|\mathcal {A}(\alpha ,\tilde \tau )-\mathcal {A}(\theta ,\tilde \tau )|^2[\mathcal {A}(\alpha ,\tilde \tau )-\mathcal {A}(\theta ,\tilde \tau )]} {1-\cos (\alpha -\theta )} \,\textrm {d}{\alpha } -|\mathcal {A}(\theta ,\tilde \tau )|^2\mathcal {A}(\theta ,\tilde \tau ) \nonumber \,, \end{align}

where we have chosen $\mathcal {B}=\mathcal {T}_2-|\mathcal {A}|^2\mathcal {A}$ above.

Equation (2.13) models the onset of filamentation on circular vortex patches at any axial position or ‘height’ $z_0$ , including the planar limit $z_0\to 1$ . In particular, it does not explicitly depend on the values of the vorticity $\omega _N$ and $\omega _S$ to the north and south of the interface. Of course, the vorticity jump requires $\omega _N-\omega _S=\omega$ , and the vanishing of the mean vorticity over the sphere (by Stokes’ theorem) requires $(1-z_0)\omega _N+(1+z_0)\omega _S=0$ . Combining these relations gives

(2.14) \begin{equation} \omega _N=\tfrac 12(1+z_0)\omega \qquad \mathrm{and} \qquad \omega _S=-\tfrac 12(1-z_0)\omega \,, \end{equation}

showing that, in general, the average vorticity near the interface

(2.15) \begin{equation} \tfrac 12 (\omega _N+\omega _S)=z_0\omega \,, \end{equation}

depends on $z_0$ . This average vorticity induces a mean shear at the position of the patch boundary unless $z_0=0$ , yet this shear plays no role in the filamentation equation (2.13) above, apart from changing the definition of the dimensionless time $\tilde \tau$ .

There is also a spectral form of this equation. When deducing $\mathcal {T}_1=\mathcal {T}_2$ above, we found that the Fourier coefficients $a_m(\tilde \tau )$ of $\mathcal {A}(\tilde \theta ,\tilde \tau )$ evolve according to

(2.16) \begin{equation} \frac {\textrm {d} a_m}{\textrm {d}\tilde \tau }=\tfrac 12 \mathrm {i} m \sum _{\substack {n\gt 0,p\gt 0,\\n+p\gt m}} (n+p-|n-m|-|p-m|-2)a_n a_p \bar {a}_{n+p-m} \,. \end{equation}

(This corrects (A12) in D88c, where the term involving $z_0^2$ there is incorrect.) Notably, $n+p-|n-m|-|p-m|-2=2(m-1)$ whenever both $n\geqslant m$ and $p\geqslant m$ . In particular, this shows that $a_1$ is an invariant of the dynamics: $\textrm {d} a_1/\textrm {d}{\tilde \tau }=0$ . It can be shown that this invariant corresponds to conservation of the $x$ and $y$ components of the angular impulse vector on the sphere

(2.17) \begin{equation} {\boldsymbol {J}}={\frac {1}{2}}\omega \oint _{\mathcal {C}} {\boldsymbol {x}}\times \textrm {d}{\boldsymbol {x}} ;\end{equation}

(see Appendix B of Polvani & Dritschel 1993). Additionally, the ‘momentum’ $\mathsf {\textit P}$ , ‘mass’ $\mathsf {\textit M}$ and (kinetic) energy $\mathsf {\textit E}$ are all invariant (Constantin et al. 2024); the momentum and mass are simple sums over spectral coefficients

(2.18) \begin{equation} \mathsf {\textit P}=\sum _{m\gt 0} |a_m|^2 \quad \textrm {and}\quad \mathsf {\textit M}=\sum _{m\gt 0} |a_m|^2/m \,, \end{equation}

while the energy has a more complex form (see Constantin et al. 2024 for details).

2.2. The line

In an appropriate limit, (2.13) (or (2.16)) also applies to the onset of filamentation for weakly nonlinear disturbances to a straight interface, say $y=0$ , on $\mathbb {R}^2$ . Guided by the analysis in Appendix B of D88c, we substitute $\theta =-\kappa x$ , $\alpha =-\kappa x'$ and $\mathcal {A}=\kappa \hat {\mathcal {A}}$ into (2.13). Then, we take the limit $\kappa \to 0$ , assuming $|\alpha -\theta | \ll 1$ . After dividing by $\kappa$ and dropping the hat on $\hat {\mathcal {A}}$ , we find

(2.19) \begin{align} \frac {\partial \mathcal {A}}{\partial \tilde \tau }&=\frac {\partial \mathcal {B}}{\partial x} \nonumber \\ \mathcal {B}(x,\tilde \tau )&=\frac {1}{4\pi }\int _0^{2\pi } \frac {|\mathcal {A}(x',\tilde \tau )-\mathcal {A}(x,\tilde \tau )|^2 [\mathcal {A}(x',\tilde \tau )-\mathcal {A}(x,\tilde \tau )]} {1-\cos (x'-x)}\,\textrm {d}{x'} , \end{align}

assuming periodicity in $x$ on the interval $[0,2\pi ]$ . Without periodicity, one would need to replace $1-\cos (x'-x)$ by $(1/2)(x'-x)^2$ and extend the integration over $x'$ to the whole real line (this is effectively (B18) in D88c). This equation has also been re-derived by Biello & Hunter (2010), but in the differential form involving $\mathcal {T}_1$ in (2.9) (with $\theta$ replaced by $-x$ ) – see Constantin et al. (2024) for further details.

The only important difference from the spherical case is the absence of the ‘curvature’ term $|\mathcal {A}|^2\mathcal {A}$ . There is also an unimportant sign change due to changing from spherical to planar Cartesian coordinates. The mode amplitudes (or Fourier coefficients) $a_k(\tilde \tau )$ in

(2.20) \begin{align} \mathcal {A}(x,\tilde \tau )=\sum _{k=1}^{\infty } a_k(\tilde \tau ) \textrm{e}^{\mathrm {i} k x} ,\end{align}

evolve according to

(2.21) \begin{equation} \frac {\textrm {d} a_k}{\textrm {d}\tilde \tau }=-\tfrac 12 \mathrm {i} k \sum _{\substack {n\gt 0,p\gt 0,\\n+p\gt k}} (n+p-|n-k|-|p-k|)a_n a_p \bar {a}_{n+p-k} \,. \end{equation}

Compared with (2.16), the ‘ $-2$ ’ in the bracket there is now missing — this came from the curvature term $|\mathcal {A}|^2\mathcal {A}$ in the circular case. Now, $a_1$ in general varies in time: none of the $a_k$ are invariant. However, momentum $\mathsf {\textit P}$ and mass $\mathsf {\textit M}$ as defined in (2.18) are still invariant, as is the energy $\mathsf {\textit E}$ (see Constantin et al. 2024).

2.3. Local, self-similar evolution

Motivated by the results in § 3 below, we seek a local asymptotic approximation of the above equations, (2.13) and (2.19), depending on a single similarity variable combining space and time. This asymptotic approximation describes, we conjecture, the moments before the finite-time singularity in wave slope observed in § 3. That is, the asymptotic solution is invariant in the similarity variable, but contains a wave-slope singularity in the original space and time variables.

The similarity variable, $\xi$ , has the same form as that which features in the fundamental solution of the heat conduction equation, namely

(2.22) \begin{equation} \xi =\frac {\theta -\theta _{{max}}(\tilde \tau )}{\sqrt {\delta }} \quad \textrm {with}\, \delta =1-\tilde \tau /\tilde \tau _c, \end{equation}

for the case of the sphere (or replace $\theta$ by $x$ for the case of the line). We next seek an approximate solution of the form

(2.23) \begin{equation} \mathcal {A}(\theta ,\tilde \tau )=\delta ^{\mathrm {i} \mu } \psi (\xi ) \,, \end{equation}

for some constant $\mu$ . This is the only self-similar form consistent with the equations of motion: the imaginary exponent of $\delta$ in the prefactor guarantees that $|\mathcal {A}|_{max }$ does not vary significantly near the singularity time, as observed in the results below. Moreover, this form leads to $|\partial \mathcal {A}/\partial \theta |_{max }\propto 1/\sqrt {\tilde \tau _c-\tilde \tau }$ , again as observed.

An equation for $\psi (\xi )$ can be obtained as follows. Take $\theta =\theta _{max }+\delta ^{{1}/{2}}\xi$ and $\alpha =\theta _{max }+\delta ^{{1}/{2}}\xi '$ in (2.13). Then, assuming $\delta \ll 1$ while $\xi '-\xi =\mathcal {O}(1)$ , one can show that $\mathcal {B}=\delta ^{\mathrm {i}\mu -{{1}/{2}}}b(\xi )$ to leading order, where

(2.24) \begin{equation} b(\xi )=-\frac {1}{2\pi }\int _{-\infty }^{\infty } \frac {|\psi (\xi ')-\psi (\xi )|^2 [\psi (\xi ')-\psi (\xi )]} {(\xi '-\xi )^2}\,\textrm {d}{\xi '}; \end{equation}

(assuming sufficiently fast far-field decay of the integrand). Notably, the curvature term $|\mathcal {A}|^2\mathcal {A}$ in (2.13) is $\mathcal {O}(\delta ^{{{1}/{2}}})$ smaller, so is neglected here. (This curvature term is absent for the case of the line.) The right-hand side of the evolution equation (2.13), namely $\partial \mathcal {B}/\partial \theta$ , is thus $\delta ^{\mathrm {i}\mu -1}\textrm {d}{b}/\textrm {d}{\xi }$ to leading order. The left-hand side evaluates to

(2.25) \begin{equation} \mathrm {i}\mu \delta ^{\mathrm {i}\mu -1}\psi + \delta ^{\mathrm {i}\mu }\frac {\textrm {d}\psi }{\textrm {d}\xi }\frac {\partial \xi }{\partial \tilde \tau } \,. \end{equation}

However, since $\xi =(\theta -\theta _{max })\delta ^{-{\frac {1}{2}}}$ , we have

(2.26) \begin{equation} \frac {\partial \xi }{\partial \tilde \tau }=-\delta ^{-{\frac {1}{2}}}\frac {\textrm {d}\theta _{max }}{\textrm {d}\tilde \tau }+\delta ^{-1}\frac {\xi }{2\tilde \tau _c} \,. \end{equation}

For $\delta \ll 1$ , we can neglect the $\mathcal {O}(\delta ^{-{\frac {1}{2}}})$ term, implying that the left-hand side of (2.13) is

(2.27) \begin{equation} \delta ^{\mathrm {i}\mu -1}\left (\mathrm {i}\mu \psi +\frac {\xi }{2\tilde \tau _c}\frac {\textrm {d}\psi }{\textrm {d}\xi }\right ), \end{equation}

to leading order. This has the same $\delta ^{\mathrm {i}\mu -1}$ prefactor as the right-hand side, so upon cancelling this prefactor we obtain an equation entirely in terms of $\xi$

(2.28) \begin{equation} \mathrm {i}\mu \psi +\frac {\xi }{2\tilde \tau _c}\frac {\textrm {d}\psi }{\textrm {d}\xi }= -\frac {1}{2\pi }\frac {\textrm {d}}{\textrm {d}\xi }\int _{-\infty }^{\infty } \frac {|\psi (\xi ')-\psi (\xi )|^2 [\psi (\xi ')-\psi (\xi )]} {(\xi '-\xi )^2}\,\textrm {d}{\xi '} \,. \end{equation}

We conjecture that this equation describes the local blow up of wave slope on vorticity interfaces, and that almost all initial conditions are attracted to this blow-up solution in finite time. For the case of the line, we arrive at the same equation, except that now the complex conjugate of $\psi$ satisfies (2.28).

Note, we can further re-scale the similarity variable, using $X=\xi /\sqrt {2\tilde \tau _c}$ , so that the solution $\psi (X)$ depends on only one free real constant, $\lambda =2\tilde \tau _c\mu$

(2.29) \begin{equation} \mathrm {i}\lambda \psi +X\frac {\textrm {d}\psi }{\textrm {d}{X}}= -\frac {1}{2\pi }\frac {\textrm {d}}{\textrm {d}{X}}\int _{-\infty }^{\infty } \frac {|\psi (X')-\psi (X)|^2 [\psi (X')-\psi (X)]} {(X'-X)^2}\,\textrm {d}{X'} \,. \end{equation}

Supplemented by the requirement that the solution be normalisable, i.e.

(2.30) \begin{align} \int _{-\infty }^{\infty }|\psi (X)|^2\,\textrm {d}{X}=1\, ,\end{align}

this can be viewed as an eigen-problem, with $\lambda$ (real) serving as the eigenvalue.

3. Results

3.1. Numerical approach

For simulating the evolution of the weakly nonlinear equations (2.13) or (2.19), we have found that the most robust approach is to solve for the spectral coefficients directly, using (2.16) or (2.21), with the sums over $n$ and $p$ truncated to ensure $n$ , $p$ and $n+p-m$ all lie in the range $[1,M]$ (substitute $k$ for $m$ in the case of the line). This is more computationally intensive than using the differential or integral form of (2.13) or (2.19), but it is stable over the times of interest. However, numerical stability requires the introduction of a damping term of the form $-\nu \,(m/M)^p a_m$ on the right-hand side of (2.16) or (2.21), where

(3.1) \begin{equation} \nu (\tilde \tau )=\frac {\varepsilon \,r_{max }(\tilde \tau )}{\alpha } \,. \end{equation}

By careful experimentation, we have found that choosing $p=6$ and $\varepsilon =10^{-2}$ ensures that the upper end of the variance spectrum $|a_m|^2$ remains small compared with lower and mid-range over the entire simulation period. In effect, the highest wavenumber is damped by the factor $e^{-\varepsilon }$ every time step when the time step is limited by $r_{max }$ .

We employ a standard fourth-order Runge–Kutta time integration scheme, with time-step adaptation. The time step is limited to

(3.2) \begin{equation} \Delta \tilde \tau =\min \left (\alpha /r_{max },\,\beta /M,\,\tilde \tau _{save}-\tilde \tau \right ) \quad \textrm {with}\, r_{max }=|\partial ^2{\mathcal {A}_r}/\partial {\tilde \tau }\partial {\theta }|_{max } , \end{equation}

where $\tilde \tau _ {save}$ is the next data save time and $\mathcal {A}_r$ is the real part of $\mathcal {A}$ . Here, we have chosen $\alpha =0.2$ (or $\alpha =0.1$ in the case of the line whose evolution is faster), but results with half this value are almost indistinguishable. The time step $\alpha /r_{\it max }$ attempts to resolve the increasingly high frequencies on the approach to any singularity. We have also chosen $\beta =2048\times 10^{-3}$ ; the time step $\beta /M$ acts as a Courant–Friedrichs–Lewy (CFL) constraint. These criteria for limiting the time step were deduced by trial and error for a wide range of initial conditions.

3.2. The sphere

We focus on one example in this case which nonetheless typifies the behaviour seen in many others. The amplitude $\mathcal {A}$ is initialised with only two non-zero spectral coefficients, specifically $a_2=\sqrt {3}/2$ and $a_3=0.25\textrm {e}^{\mathrm {i}\pi /3}$ . Then, the ‘momentum’ $\mathsf {\textit P}=|a_2|^2+|a_3|^2$ equals $1$ , and $\mathsf {\textit P}=\sum _m |a_m|^2$ remains equal to $1$ for all time $\tilde \tau$ . We have performed simulations with mode truncations $M=256$ , $512$ , $1024$ and $2048$ . All give similar results, but the highest resolution allows a better resolution of the approach to a singularity in wave slope.

Figure 2. Time evolution of the wave amplitude $\mathcal {A}$ (a, with real and imaginary parts in blue and red, respectively) for a circular vorticity interface on a sphere, together with the corresponding power spectrum $|a_m|^2$ (b) for 5 selected times $\tilde \tau$ (increasing downwards). Initially, just $a_2$ and $a_3$ are non-zero.

The evolution of $\mathcal {A}$ and of its spectrum are shown at a few selected times in figure 2. The left column shows the real and imaginary parts of the amplitude $\mathcal {A}(\theta ,\tilde \tau )$ , while the right column shows the spectrum $|a_m(\tilde \tau )|^2$ . In general, the waves comprising the interface propagate to the left, but they also steepen. In the spectrum, this is associated with the development of a power-law form $|a_m|^2 \propto m^{q}$ , with $q\approx -3$ , over a large range of wavenumbers $m$ . At the final reliable time shown ( $\tilde \tau =0.355$ ), the amplitude $\mathcal {A}$ exhibits a near discontinuous variation, which is nonetheless marginally resolved (as seen, e.g., in a zoom at this time, discussed below).

Figure 3. Time evolution of various diagnostics, as labelled, from $\tilde \tau =0$ to $0.355$ (the last reliable time). In the figure for $\textrm {d}\theta _{max }/\textrm {d}\tilde \tau$ , 1-2-1 averages (replacing data values, say $f_i$ , by $(f_{i-1}+2f_i+f_{i+1})/4$ ) were repeated 64 times to remove most of the noise occurring around the maximum near $\tilde \tau =0.32$ (endpoint values were replaced by linear extrapolation of the averaged interior data points). This noise arises from the imprecision in locating $\theta _{{max}}$ .

We next provide evidence that the solution is approaching a finite-time singularity in the wave slope $|\partial \mathcal {A}/\partial \theta |_{\it max }$ . This wave slope is determined by finding the $\theta$ grid point having the largest value of $|\partial \mathcal {A}/\partial \theta |$ , then fitting a parabola through this value and the values on either side of this grid point. The maximum in this parabola occurs at $\theta =\theta _{ {max}}$ , which varies with $\tilde \tau$ . The same procedure is also used to determine the largest value of $|\mathcal {A}|$ , to better understand the nature of the solution near the conjectured singularity. These diagnostics are shown in figure 3. One sees a dramatic rise in the maximum wave slope (upper right panel) near the conjectured singularity time, whereas other diagnostics are tending to finite values. In particular, the maximum amplitude $|\mathcal {A}|_{{max}}$ hardly varies over the entire evolution. Also, the location of maximum wave slope moves at a moderate speed, between $-10$ and $-8$ approximately.

Figure 4. Time evolution of the maximum wave slope $s(\tilde \tau )=|\partial \mathcal {A}/\partial \theta |_{max }$ together with a fit to $\sqrt {c/(\tilde \tau _c-\tilde \tau )}$ (a), and the function $f(\tilde \tau )=s^2(\tilde \tau _c-\tilde \tau )-c$ (b) which would be zero for a perfect fit.

The maximum wave slope $s(\tilde \tau )=|\partial \mathcal {A}/\partial \theta |_{max }$ appears to exhibit a square-root singularity, i.e. $s\sim \sqrt {c/(\tilde \tau _c-\tilde \tau )}$ for $\tilde \tau _c\approx 0.35748$ , as shown in figure 4. The singularity time $\tilde \tau _c$ and the constant $c$ were determined by a least-squares fit over the period $0.345\leqslant \tilde \tau \leqslant 0.355$ . Specifically, $\tilde \tau _c$ and $c$ were determined by minimising $\sum w(\tilde \tau _j)f^2(\tilde \tau _j)$ for discrete times $\tilde \tau _j$ , where $f=s^2(\tilde \tau _c-\tilde \tau )-c$ and $w=s$ , a weight chosen to favour data closer to the singularity time. The root-mean-square (r.m.s.) error in the fit over this period is only $0.01189$ (note: $c\approx 0.72935$ ).

Figure 5. Zoom of a portion of the curve $\mathcal {A}_r(\theta ,\tilde \tau )$ at $\tilde \tau \approx 0.3553927$ (blue) compared with a contour dynamics simulation (black) at the equivalent time, here $362$ linear wave periods (see the text for details).

To appreciate how well the weakly nonlinear theory captures the onset of filamentation, it is compared in figure 5, at a time just beyond the last time shown in figure 2, with the results of a full contour dynamics simulation solving (2.1) numerically. The contour dynamics simulation was initialised with $z_0=0$ (i.e. the equator) and $\rho (\theta ,0)=2a\mathcal {A}_r(\theta ,0)$ in (2.2), with amplitude $a=1/40$ . The vorticity jump across the interface was taken to be $\omega =4\pi$ so that one unit of time corresponds to one full linear wave oscillation. The numerical parameters are now standard and may be found in Dritschel & Fontane (2010), apart from the large-scale length $L$ used to control the point resolution (here, we take $L=6.25\pi /M\approx 0.009587$ , giving 3936 initial contour nodes). The agreement in figure 5 is excellent, apart from a small discrepancy near the point of maximum wave slope. At earlier times, the solutions overlap within the plotted line width.

Figure 6. A portion of the interface in the full contour dynamics simulation at successive linear wave periods, from $t=362$ at the top to $t=376$ at the bottom. Here, $\rho (\theta ,t)$ is plotted over the range $1.85\leqslant \theta \leqslant 2$ . Note that, between the times shown, the interface exhibits a relatively fast oscillation over the linear wave period, analogous to that illustrated in figure 1.

Shortly after this time, the waves on the interface do break in the contour dynamics simulation, as shown in figure 6. Filamentation starts between $t=366$ and $367$ (note the interface shape changes greatly between each period, due to the underlying linear oscillation). In terms of the long time scale, filamentation occurs between $\tilde \tau =0.3593$ and $0.3603$ , a little later than the time $\tilde \tau _c\approx 0.35748$ indicated by the singularity fit in the weakly nonlinear theory.

The behaviour seen in this example appears to be common to any (multiple-harmonic) initial condition: inevitably the maximum wave slope appears to blow up in finite time, and this is associated with the spectrum filling out, possibly tending to the form $|a_m|^2 \sim m^{-3}$ at the singularity time. Moreover, the blow up is local, occurring over a small range in $\theta$ . This observation motivated the derivation of the self-similar equation (2.28), whose solution would have the same form. However, we have been unable to solve (2.28) directly, due to its nonlinear and non-local character. In lieu, we next provide evidence that the observed behaviour in figure 2 is consistent with the self-similar form proposed in (2.23). To this end, starting from a guess for the constant $\mu$ , we define

(3.3) \begin{equation} \psi (\xi )=\frac {\int _{\tilde \tau _1}^{\tilde \tau _2} \delta ^{-\mathrm {i}\mu -{\frac {1}{2}}}\mathcal {A}(\theta _{max }+\delta ^{{\frac {1}{2}}}\xi ,\tilde \tau )\textrm {d}\tilde \tau } {\int _{\tilde \tau _1}^{\tilde \tau _2}\delta ^{-{\frac {1}{2}}}\textrm {d}\tilde \tau }, \end{equation}

as a guess for the form of $\psi$ (here integrals over time are discretised as sums). We use $\delta ^{-{\frac {1}{2}}}$ to preferentially weight the data closer to the singularity time. We then seek the constant $\mu$ which minimises

(3.4) \begin{equation} H(\mu )=\frac {\int _{\tilde \tau _1}^{\tilde \tau _2} \delta ^{-{\frac {1}{2}}}\int _{-\xi _{max }}^{\xi _{max }}w(\xi )|\delta ^{-\mathrm {i}\mu }\mathcal {A}(\theta _{max }+\delta ^{{\frac {1}{2}}}\xi ,\tilde \tau )-\psi (\xi )|^2\textrm {d}\xi \textrm {d}\tilde \tau } {\int _{\tilde \tau _1}^{\tilde \tau _2}\delta ^{-{\frac {1}{2}}}\textrm {d}\tilde \tau \int _{-\xi _{max }}^{\xi _{max }}w(\xi )\textrm {d}\xi } ,\end{equation}

where $w(\xi )=1+\cos (\pi \xi /\xi _{max })$ is a spatial weight favouring values in the centre of the interval (and vanishing at the endpoints). Again, spatial integration is done by discrete sums (here over 1024 intervals in $\xi$ ). In practice, $H(\mu )$ exhibits a single, simple minimum at the target value of $\mu$ .

Figure 7. The estimated self-similar solution $\psi (\xi )$ (real part in cyan, imaginary part in magenta), together with scaled numerical profiles (see the text) at times $\tilde \tau =0.345$ , $0.347$ , $0.349$ , $0.351$ , $0.353$ and $0.355$ (real part in blue, imaginary part in red, with fading backwards in time).

We have chosen the time interval between $\tilde \tau _1=0.345$ and $\tilde \tau _2=0.355$ (including both endpoints) to do the analysis, because this interval provides the best fit to the observed wave-slope growth in figure 4. We have also chosen $\xi _{max }=1$ to limit periodicity effects coming from the range in $\theta$ when sampling $\mathcal {A}$ in (3.3) and (3.4). With these choices, we find $\mu \approx 0.2307$ where $H(\mu )\approx 0.005163$ . For this value of $\mu$ , we can construct $\psi (\xi )$ from (3.3) and compare this with the scaled numerical profiles, specifically with $\delta ^{-\mathrm {i}\mu }\mathcal {A}(\theta _{max }+\delta ^{{{1}/{2}}}\xi ,\tilde \tau )$ , at selected times in the sampling period. This is done in figure 7. While the collapse is not perfect, it is suggestive of the existence of a self-similar solution, $\psi (\xi )$ . The neglect of terms of $\mathcal {O}(\delta ^{{\frac {1}{2}}})$ likely results in some of the scatter observed. For instance, the neglected term in (2.26) involving $\textrm {d}\theta _{max }/\textrm {d}\tilde \tau$ is questionable, since this derivative is nearly $-8$ just before the singularity time, whereas $\delta ^{{\frac {1}{2}}}$ varies from $0.1117$ to $0.04983$ over the time period analysed. This neglected term is not necessarily small. Nonetheless, the results in figure 7 at least provide a hint to the possible form of the solution to (2.28).

3.3. The line

Equation (2.19) governing the onset of filamentation on the periodic line is nearly identical to that on the sphere, (2.13), except for the curvature term $|\mathcal {A}|^2\mathcal {A}$ in the definition of $\mathcal {B}$ , and a sign difference. In fact, the complex conjugate of $\mathcal {A}$ for the line satisfies the same equation as $\mathcal {A}$ for the sphere, after dropping the curvature term. Hence, to compare the sphere and line most closely, we use the same initial condition (after complex conjugation) in spectral form, namely $a_2(0)=\sqrt {3}/2$ and $a_3(0)=0.25\textrm {e}^{-\mathrm {i}\pi /3}$ , with all other coefficients zero.

Figure 8. Time evolution of the wave amplitude $\mathcal {A}$ (a, with real and imaginary parts in blue and red respectively) for a vorticity interface on the periodic line, together with the corresponding power spectrum $|a_k|^2$ (b) for 5 selected times $\tilde \tau$ (increasing downwards). Initially, just $a_2$ and $a_3$ are non-zero. Compare with figure 2 for the circular case.

The evolutions of $\mathcal {A}$ and of its spectrum $|a_k(\tilde \tau )|^2$ are shown at a few selected times in figure 8. As for the sphere, the waves generally propagate to the left and steepen. The spectrum develops a power-law form $\propto k^{q}$ , with the exponent $q$ tending to $-3$ at the latest time. The same behaviour was found for the sphere. At the final reliable time shown ( $\tilde \tau =0.173$ ), the amplitude $\mathcal {A}$ exhibits a near discontinuous variation. However, the details differ somewhat from those seen for the sphere in figure 2, for instance the ‘kink’ now appears in the real part of $\mathcal {A}$ rather than the imaginary part. Also, the evolution is nearly twice as fast. Finally, unlike for the sphere, the lowest harmonic $a_1$ varies in time.

Figure 9. Time evolution of various diagnostics, as labelled, from $\tilde \tau =0$ to $0.173$ (the last reliable time). In the figure for $\textrm {d}{x}_{\it max }/ {\rm d}{t}$ , 1-2-1 averages were repeated 8 times to remove most of the noise occurring around the maximum near $\tilde \tau =0.16$ . Compare with figure 3 for the circular case.

We next examine the evolution of the maximum wave slope $|\partial \mathcal {A}/\partial {x}|_{max }$ in figure 9, along with other diagnostics. There is again an indication that the system is approaching a finite-time singularity in wave slope. Near the final time shown, the maximum wave amplitude $|\mathcal {A}|_{max }$ hardly varies. The location of maximum wave slope, $x_c(\tilde \tau )$ , moves at a moderate speed, even faster than in the case of the sphere (indeed almost twice as fast).

Figure 10. Time evolution of the maximum wave slope $s(\tilde \tau )=|\partial \mathcal {A}/\partial {x}|_{max }$ together with a fit to $\sqrt {c/(\tilde \tau _c-\tilde \tau )}$ (a), and the function $f(\tilde \tau )=s^2(\tilde \tau _c-\tilde \tau )-c$ (b) which would be zero for a perfect fit. Compare with figure 4 for the circular case.

Like in the case of the sphere, the maximum wave slope $s(\tilde \tau )=|\partial \mathcal {A}/\partial {x}|_{max }$ appears to exhibit a square-root singularity $s\sim \sqrt {c/(\tilde \tau _c-\tilde \tau )}$ , now for $\tilde \tau _c\approx 0.17478$ , as shown in figure 10. The singularity time $\tilde \tau _c$ and the constant $c$ were determined by a least-squares fit over the period $0.163\leqslant \tilde \tau \leqslant 0.173$ , using the same approach used for the sphere. The r.m.s. error in the fit over this period is only $0.00788$ (note: $c\approx 0.65997$ ).

Figure 11. The estimated self-similar solution $\psi (\xi )$ (real part in cyan, imaginary part in magenta), together with scaled numerical profiles (see the text) at times $\tilde \tau =0.163$ , $0.165$ , $0.167$ , $0.169$ , $0.171$ and $0.173$ (real part in blue, imaginary part in red, with fading backwards in time). Compare with figure 7 for the circular case.

We next provide evidence for a self-similar solution of the form $\mathcal {A}=\delta ^{\mathrm {i} \mu } \psi (\xi )$ , where $\xi$ and $\delta$ have the same definitions given in (2.22) except $\theta$ is replaced by $x$ , and the complex conjugate of $\psi$ satisfies (2.28). As for the sphere, we fit the scaled numerical data over a time period extending from $\tilde \tau _1=0.163$ to $\tilde \tau _1=0.173$ to estimate the form of $\psi (\xi )$ and the exponent $\mu$ . That fit gives $\mu \approx 0.3150$ (cf. $\mu \approx 0.2307$ for the sphere), and $H(\mu )\approx 0.005163$ in (3.4). The estimated form of $\psi$ along with the scaled numerical data $\delta ^{-\mathrm {i}\mu }\mathcal {A}(x_{max }+\delta ^{{\frac {1}{2}}}\xi ,\tilde \tau )$ are shown in figure 11. The form of $\psi$ compares well with the spherical case in figure 7, though there are some differences. These differences, and the scatter in the fit, may result from the impact of higher-order terms, formally $\mathcal {O}(\delta ^{{\frac {1}{2}}})$ smaller but in practice potentially comparable due to the computational difficulty in getting closer to the singularity time. Nonetheless, these results hint at the likely form of the solution to the self-similar equation (2.28).

3.4. Stability of travelling waves

Figure 12. Time evolution of various diagnostics for a simulation starting from a weakly perturbed travelling wave having $\epsilon =0.1$ (see the text). See figure 3 for the analogous diagnostics in the case of a larger initial perturbation having $\epsilon =0.5$ .

Figure 13. As in figure 12 but for a slightly larger disturbance amplitude, here $\epsilon =0.12$ .

All single-mode disturbances, which take the form

(3.5) \begin{align} \mathcal {A}(\theta ,\tilde \tau )=a_k \textrm {e}^{\mathrm {i} k\theta -k(k-1)|a_k|^2\tilde \tau }\,, \end{align}

are exact solutions of (2.13) or (2.16), and translate in $\theta$ at speed $-(k-1)|a_k|^2$ (Constantin et al. 2024). (In the case of the line, the $k-1$ factor is replaced by $k$ .) Here, we briefly examine the stability of a few of these solutions numerically. We initialise as before by choosing another mode $p\gt k$ with initial amplitude $a_p(0)=\epsilon \textrm{e}^{\mathrm {i} \pi /3}$ , then chose $a_k(0)=\sqrt {1-\epsilon ^2}$ so that $\mathsf {\textit P}=\sum _m |a_m|^2=1$ , without loss of generality. Again, we take $k=2$ and $p=3$ , but take $\epsilon =0.1$ , which is 5 times smaller than in the previous example. (Other values are discussed below.)

Figure 14. Late-time evolution of the relative errors in momentum $\mathsf {\textit P}$ and mass $\mathsf {\textit M}$ . A subscript ‘0’ refers to their initial values. Prior to $\tilde \tau \approx 34$ , errors are $\mathcal {O}(10^{-11})$ but likely smaller because the spectral coefficients $a_m$ were only saved to 11 digit accuracy.

Figure 12 shows various diagnostics including the maximum wave slope $|\partial \mathcal {A}/\partial \theta |_{max }$ . Comparing with the case with $\epsilon =0.5$ in figure 3, we see that smaller $\epsilon$ delays the onset of filamentation, which appears to occur after a series of growing oscillations; this behaviour is found also for $\epsilon =0.11$ and $0.12$ , see e.g. Figure 13, whereas $\epsilon =0.13$ and larger $\epsilon$ exhibit a single diverging peak in wave slope. The oscillations correspond to the almost periodic amplification and reduction of the tail of the power spectrum, with each amplification being larger than the last. The wave form $\mathcal {A}$ begins to develop a kink and a near discontinuity by $\tilde \tau =37.2$ for $\epsilon =0.1$ . The recovery of $|\partial \mathcal {A}/\partial \theta |_{max }$ after this time is likely to be spurious – significant errors in the conserved momentum $\mathsf {\textit P}$ and mass $\mathsf {\textit M}$ develop especially around $\tilde \tau =37$ , when the final peak occurs – see figure 14. The numerical filtering at high wavenumbers $m$ arrests the growth in wave slope, which likely diverges soon after $\tilde \tau =37$ in reality. Notably, the best-fit spectral slope $q$ climbs through $-3$ at $\tilde \tau =37.06$ , when strong numerical dissipation occurs. All evidence points to a wave-slope singularity in finite time; it appears that even weakly perturbed translating waves are eventually subject to filamentation.

4. Discussion

This paper has revisited a generic aspect of vorticity interfaces, namely the tendency for unsteady disturbances to steepen and ‘break’, resulting in ‘filamentation’ (Dritschel 1988c ). Specifically, we have studied an equation first derived in Dritschel (1988c ) describing the weakly nonlinear development of shallow (small wave slope) disturbances to circular and linear interfaces. That equation is shown to have a simpler, universal form in a re-scaled slow-time variable (see also Constantin et al. 2024 for details of its mathematical structure, invariants and some exact solutions).

We have studied the onset of filamentation on both a circular interface (applicable to interfaces on the sphere or on the plane) and a periodic linear interface. In both cases, we find generically that unsteady initial conditions tend to steepen, apparently exhibiting an inverse square root time singularity in wave slope at one location. This is the manifestation of filamentation discussed in Dritschel (1988c ), who showed that the results compare well with full contour dynamics simulations up to the point of filamentation (this is confirmed here with more accurate simulations). The significance of this finding is that vortex patch boundaries inevitably and endlessly grow in complexity due to filamentation. Almost no initial condition can avoid this fate – only some sort of regularisation can prevent or limit filamentation.

Motivated by the numerical findings, we have developed a theory for the self-similar evolution of the interface close to the singularity in wave slope. This is a cubically nonlinear, nonlocal equation for a complex function $\psi (\xi )$ in a similarity variable $\xi$ , of the form that appears in the heat equation. The equation has an eigenvalue $\mu$ related to the time-scaling factor with exponent $\mathrm {i} \mu$ used to define $\psi$ . We have not been able to solve this equation, but have provided an estimate of the form of $\psi$ by fitting numerical data close to the singularity time. Future work will target finding numerical solutions of this self-similar equation and, more generally, rigourously proving the self-similar equation is an attractor for almost all initial conditions.