# Singularity formation of vortex sheets in two-dimensional Euler equations using the characteristic mapping method

Julius Bergmann, Thibault Maurel-Oujia, Xi-Yuan Yin, Jean-Christophe

Nave, Kai Schneider

► To cite this version:

Julius Bergmann, Thibault Maurel-Oujia, Xi-Yuan Yin, Jean-Christophe Nave, Kai Schneider. Singularity formation of vortex sheets in two-dimensional Euler equations using the characteristic mapping method. Physics of Fluids, 2024, 36 (12), pp.124126. 10.1063/5.0241214. hal-04929625

HAL Id: hal-04929625

https://hal.science/hal-04929625v1

Submitted on 4 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Singularity Formation in two-dimensional Euler Equations

## Singularity formation of vortex sheets in two-dimensional Euler equations using the Characteristic Mapping Method

J. Bergmann,¹,², a) T. Maurel-Oujia,¹ X.-Y. Yin,³ J.-C. Nave,⁴ and K. Schneider¹

¹) Institut de Mathématiques de Marseille, Aix-Marseille Université, CNRS, Marseille, France

²) Institut für Strömungsmechanik und Technische Akustik, Technische Universität Berlin, Berlin, Germany

³) LMFA-CNRS, Ecole Centrale de Lyon, Lyon, France

⁴) Department of Mathematics and Statistics, McGill University, Montréal, Québec, Canada

(Dated: 14 November 2024)

The goal of this numerical study is to get insight into singular solutions of the two-dimensional (2D) Euler equations for non-smooth initial data, in particular for vortex sheets. To this end high resolution computations of vortex layers in two-dimensional incompressible Euler flows are performed using the characteristic mapping method (CMM). This semi-Lagrangian method evolves the flow map using the gradient-augmented level set method (GALS). The semi-group structure of the flow map allows its decomposition into sub-maps (each over a finite time interval), and thus the precision can be controlled by choosing appropriate remapping times. Composing the flow map yields exponential resolution in linear time, a unique feature of CMM, and thus fine-scale flow structures can be resolved in great detail. Here the roll-up process of vortex layers is studied varying the thickness of the layer showing its impact on the growth of palinstrophy and possible blow up of absolute vorticity. The curvature of the vortex sheet shows a singular-like behavior. The self-similar structure of the vortex core is investigated in the vanishing thickness limit. Conclusions on the presence of possible singularities of two-dimensional Euler equations for non-smooth initial data are drawn by tracking them in the complex plane.

<footnote>a) Electronic mail: julius.bergmann .at. univ-amu.fr</footnote>
# Singularity Formation in two-dimensional Euler Equations

## I. INTRODUCTION

The emergence of coherent vorticity is ubiquitous in high Reynolds number turbulence, as for example in shear layers or in wall bounded flow in the form of thin boundary layers detaching and traveling into the bulk flow. The dynamics of thin vortex layers in the limit of vanishing viscosity is of tremendous interest for understanding the formation of singularities in the incompressible 2D Euler equations with non-smooth initial data. In contrast to 2D Euler with smooth initial conditions, for which results on global regularity and uniqueness of the solution are well known¹, it was shown by Székelyhidi² that for vortex sheet initial data infinitely many non-stationary weak solutions exist, which moreover conserve energy. For a review on the mathematics of turbulence, we refer e.g. to Majda and Bertozzi³, Bardos and Titi⁴. Extensive literature on vortex sheets and their dynamics is available, and an exhaustive review is beyond the scope of the current paper. The textbook of Saffman⁵ gives a nice overview and contains numerous references. Singularities in two-dimensional vortex layers have been studied in Caflisch et al.⁶ and their introduction gives some overview on the specific topic of the paper. The limit of thin shear layers has been investigated before (Baker and Shelley⁷, Baker and Pham⁸, Sohn⁹, Sohn¹⁰), while Pullin and Shen¹¹ studied recently vortex-sheet solutions.

In a recent work Caflisch et al.⁶ showed by means of direct numerical simulation using a classical pseudo-spectral method that 2D vortex layers may have complex singularities and after the roll up of the layer that the vortex cores may have unbounded vorticity in the limit of infinite Reynolds numbers. We recall that Caflisch et al.⁶ did almost exclusively viscous computations by solving Navier-Stokes using Fourier pseudo-spectral methods and coupling the vortex layer thickness with the viscosity, i.e. the inverse Reynolds number. Some computations for the inviscid case, i.e. for 2D Euler, were likewise presented in Caflisch et al.¹² and Caflisch et al.⁶. However, the inviscid computations were limited to minimum layer thickness values ($\delta = 0.0141$). Note that smaller values did not work due to the necessary resolution requirements, however they also used larger values (cf. table 3 in Caflisch et al.⁶). The motion of an inviscid vortex sheet, i.e. for vanishing thickness, is governed by the Birkhoff-Rott model equation¹³,¹⁴ which develops a singularity in finite time starting from smooth initial data¹⁵. Perturbations grow due to the Kelvin-Helmholtz instability and the vortex sheet does roll up. Regularized simulations using a vortex blob
# Singularity Formation in two-dimensional Euler Equations

method have been also performed in Caflisch et al.⁶, Baker and Pham⁸ and compared with Navier-Stokes computations.

The vanishing viscosity limit of 2D Navier-Stokes in the presence of boundaries was likewise studied in Nguyen Van Yen, Farge, and Schneider¹⁶, Nguyen Van Yen et al.¹⁷. Dipole-wall collisions were simulated and the existence of a Reynolds-independent energy dissipation rate was shown. In this context Prandtl's classical boundary layer argument was complemented, which states that both the boundary-layer thickness and dissipation rate are proportional to $Re^{-1/2}$. However for detaching boundary layers, Kato's scaling was shown to be more appropriate than Prandtl's scaling, which implies that the boundary layer scales with $Re^{-1}$. Some reviews of 2D flows with walls can be found in Clercx and Van Heijst¹⁸ and mathematical analysis of weak solutions of the 2D Navier-Stokes equations in bounded domains, in the vanishing viscosity limit, in Constantin et al.¹⁹.

State of the art for solving Navier-Stokes or Euler equations numerically with high precision are pseudo-spectral methods²⁰,²¹ which have been also extensively used for investigating nearly singular solutions of the 3D Euler equations²²,²³. A Cauchy-Lagrange method for computing 2D Euler flows was proposed in Podvigina, Zheligovsky, and Frisch²⁴. This semi-Lagrangian method exploits the time-analyticity of fluid particle trajectories and was shown to be more efficient than pseudo-spectral computations. However, detailed singularity studies have not been reported so far. An even more powerful tool for solving the incompressible Euler equations is the Characteristic Mapping Method (CMM). Evolving the flow map with a Gradient Augmented Level Set Method, developed in²⁵⁻²⁷, one can decompose the long time deformation into short time sub-maps due to the semi-group structure of the flow map. This yields a numerical scheme with exponential resolution in linear time developed for linear advection in Mercier, Yin, and Nave²⁸ and 2D Euler in²⁹, allowing to capture the exponential growth of vorticity gradients. The implementation of the method has global third order convergence in space and time and its efficiency has been demonstrated in comparison with spectral and Cauchy-Lagrange methods in²⁹. More recently an extension to 3D incompressible Euler flows has been proposed in³⁰. The compositional adaptivity of CMM is an essential feature which allows detailed insight into the small scales of the solution without using prohibitive numerical resolutions.

The goal of the present paper is to study inviscid vortex layers with regard to possible singularities using high resolution CMM computations. The motivation is triggered by the
## Singularity Formation in two-dimensional Euler Equations

vortex layer simulations of Caflisch et al.⁶ where, with a few exceptions, viscous results have been presented and where the layer thickness is coupled with the viscosity. As presented in Caflisch et al.¹², for too thick vortex layers the vortex merging process starts to interfere with the formation of the two vortex blobs with spiral arms. Here we compute solutions of 2D incompressible Euler flows and study the dynamics of these extremely thin vortex layers in the vanishing thickness limit and investigate possible singularities. The aim is to get some insights about possible singularities in 2D incompressible Euler for non-smooth initial data. Curvature and vortex strength of the vortex center-line are analyzed and a temporal and spatial normalization unveils the dynamics for vanishing thickness limit, as well as the investigation of singularities in the complex plane. The palinstrophy growth and energy spectra show and distinguish the impact of both the forming vortices and vortex merger process. Thanks to the high resolution capabilities of CMM we get insight into the fine scale structure of vortex cores and their dynamics in Euler flows.

The remainder of the paper is organized as follows. Set-up and initial conditions are discussed in section II. A short description of the characteristic mapping method for solving the 2D incompressible Euler equations is given in section III. Section IV gives an overview on the performed computations and numerical results are then presented in section V. A singularity analysis in the complex plane is performed in section VI. Finally, conclusions and perspectives for future work are given in section VII.

## II. GOVERNING EQUATIONS AND INITIAL CONDITION

We consider inviscid flow in a $2\pi$ periodic domain $\Omega = [-\pi, \pi] \times [-\pi, \pi]$ in the plane, governed by the 2D incompressible Euler equations. The starting point is the vorticity transport equation,

$$ \partial_t \omega + (\mathbf{u} \cdot \nabla)\omega = 0 \qquad (1) $$

where the vorticity is defined as $\omega = \nabla \times \mathbf{u}$ and $\mathbf{u}$ is the incompressible velocity, satisfying $\nabla \cdot \mathbf{u} = 0$. The curl operator is invertible, and the velocity can be computed from the vorticity, $\mathbf{u} = -\nabla \times \Delta^{-1}\omega$ using the Biot-Savart operator. Equation 1 is completed with a suitable initial condition $\omega_0(\mathbf{x}) = \omega(\mathbf{x}, t=0)$, here a regularized vortex sheet and where $\mathbf{x} = (x, y)$ with $(x, y) \in \Omega$ as the cardinal directions. All variables have been non-dimensionalized,
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_2.png" /><figcaption>FIG. 1: Initial condition $\omega_0$ for $\delta \approx 0.032$ (a) and energy spectra for different $\delta$-values (b).</figcaption></figure>

similar to Caflisch et al.⁶, as follows,

$$x = x^* \frac{1}{\lambda}, \quad t = t^* \frac{\Gamma}{\lambda^2}, \quad u = u^* \frac{\lambda}{\Gamma}, \quad \omega = \omega^* \frac{\lambda^2}{\Gamma} \quad (2)$$

With the domain length $L_x = L_y = 2\pi$ we have for the spatial scaling $\lambda = L_x/(2\pi) = 1$. The circulation $\Gamma$ equals $2\pi$, as shown below. The regularized vortex sheet initial condition proposed in Caflisch et al.⁶ reads,

$$\omega_0(x, y) = \frac{1}{\sqrt{2\pi\delta}} \exp\left(-\frac{(y - \phi(x))^2}{2\delta^2}\right) \quad (3)$$

where $\delta$ is the thickness parameter and $\phi$ a perturbation function. The initial field corresponds to a vorticity line with Gaussian cross-section of thickness $\delta$ centered around a perturbation function $\phi(x)$ shown in figure 1. The profile for $\phi(x)$ is sinusoidal with $\phi(x) = \sin(x)/2$ and $(x, y) \in \Omega$. For the perturbation, the amplitude over wavelength ratio is $1/(4\pi)$. The position of the center line $\phi$ is the curve of maximum vorticity, which oscillates around $y = 0$. For $\delta \to 0$ this vortex layer converges towards a vortex sheet used in Caflisch et al.⁶ as initial condition for the Birkhoff-Rott equation.

The energy spectrum of the initial condition describes the initial frequencies present in the flow field. It is defined as:

$$E(k, t) = \frac{1}{2} \sum_{k-1/2 \le |k| < k+1/2} |\hat{u}(k, t)|^2 \quad (4)$$

Here $\hat{\cdot}$ denotes the 2D Fourier transform. The initial condition, a vortex line regularized with a Gaussian cross-section, does in the limit of $\delta \to 0$ approach a vortex line, i.e. a Dirac
Singularity Formation in two-dimensional Euler Equations

distribution. The entrophy spectrum thus exhibits a $k^0$ scaling and consequently the energy spectrum yields a $k^{-2}$ scaling, using the relation $E(k, t) = k^{-2}Z(k, t)$. For large $\delta$ values the regularization becomes more and more visible resulting in a faster, i.e. exponential, decay for large wave-numbers, as shown in figure 1b.

For viscous flow simulations, considered in Caflisch et al.⁶, $\delta$ was related to the Reynolds number $Re = \Gamma/\nu$ via $\delta = Re^{-1/2}$, with $\Gamma = \int \omega(\mathbf{x}, t = 0)d\mathbf{x}$ being the initial circulation of the vortex layer and $\nu$ being the kinematic viscosity. Note that in the present study $\nu$ vanishes. For the given initial condition, the initial circulation is

$$ 
\begin{aligned}
\int_{\Omega} \omega_0(\mathbf{x}) d\mathbf{x} &= \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} \frac{1}{\sqrt{2\pi}\delta} \exp\left(-\frac{(y - \phi(x))^2}{2\delta^2}\right) dxdy \\
&= \int_{-\pi}^{\pi} \int_{(-\pi - \phi(x))/2\delta}^{(\pi - \phi(x))/2\delta} \frac{1}{\sqrt{\pi}} \exp(-v^2) dv dx \\
&= \int_{-\pi}^{\pi} \frac{1}{2} \left[ \text{erf}(v) \right]^{(\pi - \phi(x))/2\delta}_{(-\pi - \phi(x))/2\delta} dx \\
&= 2\pi + o(\exp(-\delta^{-2})).
\end{aligned}
\quad (5) $$

The evaluation limits $(\pm\pi - \phi(x))/2\delta$ arise from the truncation of the Gaussian profile by the periodic box (on the upper and lower sides of which there is technically a discontinuity). Here $\text{erf}(z) = \frac{2}{\sqrt{\pi}} \int_{0}^{z} \exp(-v^2) dv$ is the error function which rapidly approaches 1 for $z \to \infty$, i.e. for $\delta \to 0$. Hence, the error function can be approximated by 1 for all investigated thickness values $\delta$ to a certainty far below machine precision, justified by the expansion $\text{erf}(x) = 1 - e^{-x^2} \left( \frac{1}{\sqrt{\pi}} \left( \frac{1}{x} - \frac{1}{2x^3} + \frac{3}{4x^5} - \frac{15}{8x^7} \right) + o(x^{-8}e^{-x^2}) \right)^{31}$. The term $[\text{erf}(v)]^{(\pi - \phi(x))/\delta}_{(-\pi - \phi(x))/\delta}$ can therefore be approximated to 2 with error of order $o(\exp(-(\pi - \frac{1}{2})^2\delta^{-2}))$ using that $|\phi(x)| \le \frac{1}{2}$.

For an incompressible and inviscid flow, energy $E(t)$ and entrophy $Z(t)$ remain constant over time, while the palinstrophy does increase super exponentially³². Those quantities are defined by

$$ E(t) = \frac{1}{2} \int_{\Omega} |\mathbf{u}(\mathbf{x}, t)|^2 d\mathbf{x}, \quad (6) $$

$$ Z(t) = \frac{1}{2} \int_{\Omega} |\omega(\mathbf{x}, t)|^2 d\mathbf{x}, \quad (7) $$

$$ P(t) = \frac{1}{2} \int_{\Omega} |\nabla \omega(\mathbf{x}, t)|^2 d\mathbf{x}. \quad (8) $$

With varying thickness $\delta$, the initial energy $E(t = 0)$ was found to increase linearly with decreasing vortex sheet thickness with a limit value estimated as 1.52 as $\delta \to 0$ (figure 2),
## Singularity Formation in two-dimensional Euler Equations

The initial enstrophy $Z(t = 0)$ is shown below to scale with $\delta^{-1}$ and the initial palinstrophy $P(t = 0)$ with $\delta^{-3}$, both values thus go to infinity with vanishing vortex sheet thickness. Similar to the circulation, we get for the initial enstrophy:

$$ 
\begin{aligned}
& \frac{1}{2} \int_{\Omega} |\omega_0(\mathbf{x})|^2 d\mathbf{x} \\
&= \frac{1}{2} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} \frac{1}{2\pi\delta^2} \exp\left(-\frac{(y-\phi(x))^2}{\delta^2}\right) dxdy \\
&= \frac{1}{2} \int_{-\pi}^{\pi} \int_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta} \frac{1}{2\pi\delta} \exp(-v^2) dvdx \\
&= \frac{1}{2} \int_{-\pi}^{\pi} \frac{1}{2\sqrt{\pi}\delta} \left[\text{erf}(v)\right]_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta} dx \\
&= \frac{\sqrt{\pi}}{2\delta} + o(\exp(-\delta^{-2})).
\end{aligned} 
\quad (9) $$

Respectively, the initial palinstrophy can be computed analytically with the gradient of the initial vorticity:

$$ \nabla\omega_0 = - \left( \frac{-\phi'(x)/\delta}{1/\delta} \right) \frac{y+\phi(x)}{\delta} \omega_0(x, y) \quad (10) $$

We thus obtain

$$ 
\begin{aligned}
& \int_{\Omega} |\nabla\omega_0|^2 d\mathbf{x} \\
&= \frac{1}{2\pi\delta^3} \int_{-\pi}^{\pi} \int_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta} (1+(\phi'(x))^2) v^2 \exp(-v^2) dv dx \\
&= \frac{1}{8\pi\delta^3} \int_{-\pi}^{\pi} \left(1+\frac{\cos^2(x)}{4}\right) \left[\sqrt{\pi}\operatorname{erf}(v) - 2v \exp(-v^2)\right]_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta} dx \\
&= \frac{9\sqrt{\pi}}{32\delta^3} + o(\exp(-\delta^{-2})).
\end{aligned} 
\quad (11) $$

We bound the size of this error using $\max\phi(x)$, by $\operatorname{erf}((\pi - 1/2)/\delta) = 1 - o(\exp(-(\pi - 1/2)/\delta)^2)$ and similarly for the $2u \exp(-u^2)$ term. This then gives the initial palinstrophy $P(t = 0) = \frac{9\sqrt{\pi}}{32\delta^3} + o(\exp(-\delta^{-2}))$.

The numerical values for normalized initial enstrophy $Z(t = 0) \cdot \delta = \frac{\sqrt{\pi}}{2} \approx 0.886$ and palinstrophy $P(t = 0) \cdot \delta^3 = \frac{9\sqrt{\pi}}{32} \approx 0.499$ encountered in the simulations are consistent with the analytically derived values up to machine precision ($10^{-16}$). The same observation applies for the initial circulation, matching the value $2\pi$ to machine precision.
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_3.png" /><figcaption>FIG. 2: Initial energy $E(t = 0)$ as a function of the sheet thickness $\delta$ showing a linear decrease with a value of $\approx 1.52$ for $\delta \to 0$.</figcaption></figure>

### III. CHARACTERISTIC MAPPING FOR 2D EULER

In this section, we briefly recall the characteristic mapping method (CMM) for solving the 2D incompressible Euler equations. For details and numerical implementation we refer the reader to Yin et al.²⁹. The numerical results in this paper are computed using a GPU CUDA implementation of the method. The open source code is available on GitHub³³.

The CMM for the 2D incompressible Euler equations is based on the 2D scalar vorticity formulation (equation 1) coupled with the CMM for linear transport. We recall that the characteristic map $\chi_B$ is the backward Lagrangian flow map and can be thought of as the

<figure><img src="image_4.png" /><figcaption>FIG. 3: Illustration of the pullback operation by the characteristic map. Figure (a) shows the initial vorticity, figure (b) overlays onto figure (a) the backward map, and figure (c) shows the time $t$ vorticity obtained by pullback using the relabelling symmetry.</figcaption></figure>
## Singularity Formation in two-dimensional Euler Equations

back-to-label operator for the advection of Lagrangian particles. Indeed, for any particle trajectory $\gamma(t)$ given by,

$$ \frac{d}{dt}\gamma(t) = u(\gamma(t), t) \qquad (12) $$

with initial condition $\gamma(0) = \gamma_0$, we have that $\chi_B$ satisfies,

$$ \gamma_0 = \chi_B(\gamma(t), t). \qquad (13) $$

It follows by the method of characteristics, that any scalar quantity $\phi$ advected by the same velocity $u$ satisfies $\phi(x, t) = \phi_0(\chi_B(x, t))$, this is known as the relabelling symmetry.

In the vorticity equations, $\omega$ is a scalar-valued quantity advected by $u$ and thus has the above relabelling symmetry. The velocity field $u$ is in turn obtained from the Biot-Savart law. The coupling of the vorticity equations with the characteristic map yields the following governing equations for the numerical method:

$$ \partial_t \chi_B + (u \cdot \nabla) \chi_B = 0 \qquad (14a) $$

$$ \omega(x, t) = \omega_0(\chi_B(x, t)) \qquad (14b) $$

$$ u = -\nabla \times \Delta^{-1} \omega \qquad (14c) $$

Numerically, the advection equation for the map (14a) is discretized through the Gradient-Augmented Level-Set method$^{25-27}$ which is a semi-Lagrangian method using Hermite cubic spatial interpolation with Runge-Kutta time-stepping schemes. The vorticity evaluation step (14b) is performed on a fixed fine grid by interpolation of the Hermite cubic representation of $\chi_B$ followed by a direct evaluation of the initial condition $\omega_0$. Figure 3 illustrates the computation of the vorticity by pullback on the initial condition. Lastly, the velocity field is defined through the Biot-Savart law (14c) as the curl of the Hermite interpolant of the stream function, which we obtain from a spectral solver for the Poisson equation using Fast Fourier Transforms. The time dependence of the velocity field is obtained from a Lagrange extrapolation in time using stream function data from the three most recent time steps.

This approach has several desirable numerical properties for the problem investigated here. Firstly, the use of the advection solution operator $\chi_B$ maintains a back-to-label symmetry of the vorticity field. This ensures a non-dissipative numerical method, since there exists a coordinate transform which transports the numerical solution to the exact one. This property is beneficial to the study of singularity formation, since artificial viscosity
## Singularity Formation in two-dimensional Euler Equations

is eliminated as a potential regularization mechanism preventing blow-up. Secondly, the characteristic maps benefit from the structure of the Lie-groups of volume-preserving diffeomorphisms. Therefore, a long-time map can be represented as the composition of multiple short-time sub-maps. For a time subdivision $0 < t_1 < t_2 < ... < t_{m-1} < T$, the sub-map decomposition is an adaptive and multi-resolution representation of the full map $\chi_B$ given by

$$ \chi_B(\cdot, T) = \chi_{[t_1,0]} \circ \chi_{[t_2,t_1]} \circ \cdots \circ \chi_{[T,t_{m-1}]} \quad (15) $$

where $\chi_{[t_i, t_{i-1}]}$ represents the backwards map in the time interval $[t_{i-1}, t_i]$. On each time sub-interval, the spatial deformation is comparatively small and can be well resolved on a coarser grid, the remapping step is triggered dynamically as the numerical resolution of the sub-map coarse grids are depleted. The full map obtained from the map composition in (15) retains all scales formed by each sub-map.

The movement of the vortex sheet center line, parametrized by $x(\theta, t)$, is tracked by following the trajectories of Lagrangian particles in the vortex core. The initial position corresponds to the perturbation function $\phi(x)$, i.e. $x(\theta, t=0) = (\theta, \phi(\theta) + \pi)$, $\theta \in [0, 2\pi]$. The initial vortex center curve is discretized using $N_P$ uniformly distributed sample particles with $\theta_n = n 2\pi/N_P$. A sufficiently high number of points $N_P \ge 10^6$ ensures that all dynamics of the material line is captured even under strong elongation effects.

Their time evolution under the numerical velocity field computed from the CMM is then tracked individually using standard explicit Runge-Kutta methods. We must note here that the velocity field obtained from the CMM is only $C^0$ in space since it is defined as the curl of the Hermite cubic interpolant of the stream function. The size of the discontinuities in higher derivatives depends on the grid size $N_\psi$ used for the stream function with the jumps in $i^{th}$ derivative scaling like $\|D^4\psi\|_\infty N_\psi^{-4+i}$. This introduces a negligible error in the $L^\infty$ error for the curve position and derivatives but can be a source of noise for the regularity analysis.

## IV. PERFORMED COMPUTATIONS

Several runs with successively decreasing $\delta$-values were executed on state-of-the art graphics cards, maximizing the available usable memory. All used parameters are summarized in
<figure><img src="image_5.png" /><figcaption>FIG. 4: Evolution of initial condition with $\delta \approx 0.032$ portraying the roll-up of the two vortices. The images are centered around 0 with extent $\pi/2$.</figcaption></figure>

table I. Two grid sizes were used: A coarser one for the description of the flow map $\chi$ and velocity $u$ and one for the initial vorticity $\omega_0$ defined for each sub-map as well as the discrete vorticity used in the Biot-Savart law. The merit of this lies in the enhanced smoothness of the velocity field, obviating the need for a highly detailed representation. Additionally, fine scales of the flow captured by the flow map are retained from the composition of the individual sub-maps. The time-step $\Delta t$ is set after a Courant-Friedrich-Lewis (CFL) number of 1/3 of the coarse map to neglect its influence, however the semi-Lagrangian method can easily deal with large CFL numbers > 1 as well. The incompressibility threshold of $\delta_{inc,b}$ ensures the volume preservation of the sub-maps, i.e. we monitor that $|\det \nabla \chi_B| \le \delta_{inc,B}$. The parameter of the local stencil size $\epsilon_m$ defines the distance at which spatial derivatives for the creation of the Hermite interpolants with the GALS method are computed, together with its corresponding map update stencil order. The filter size $k_{LP}$ defines the cut-off frequency for the low-pass filter and was disabled by setting it to high frequencies of negligible energy. The fluid and embedded particle time schemes are chosen to adapted versions of the Runge-Kutta scheme with improved efficiency for the CM method. At last, $N_P$ is the number of equidistant particles embedded in the vortex center line. Further in-depth analysis and explanation of all used parameters are reported in Yin et al.²⁹, Bergmann³⁴.

The parameters were specifically chosen to maximize the coarse grid size for the given GPU memory, which in return leads to lower growth of incompressibility error and therefore more accurate flow representation over time. With too low grid resolution, a simulation with small thickness $\delta$ leads to premature emergence of Kelvin-Helmholtz (KH) instabilities along
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_6.png" /><figcaption>FIG. 5: Close-up of upper vortex with overlayed center-line (a) and corresponding local palinstrophy $P = |\nabla\omega(\mathbf{x}, t)|^2$.</figcaption></figure>

the line of high vorticity. This was observed to be suppressed with increased grid size and therefore attributed to numerically induced errors. For the highest available grid settings on the NVIDIA A100 cards, a value of $\delta \approx 0.007$ was the lowest observed value where perturbations are bounded and do not dominate the flow behavior. In simulations with smaller

<figure><img src="image_7.png" /><figcaption>FIG. 6: Comparison between our results and Caflisch et al.⁶ with $\delta \approx 0.014$ and at $t = 2.85$. They show excellent agreement. Figure (b) from Caflisch et al., "Complex singularity analysis for vortex layer flows", J. Fluid Mech., 932, A21, 29, 2021 © Cambridge University Press, reproduced with permission.</figcaption></figure>
## Singularity Formation in two-dimensional Euler Equations

<table><caption>TABLE I: Simulation parameters for the different runs on Nvidia V100 and A100 machines.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>V100</th><th>A100</th></tr></thead><tbody><tr><td>Coarse & stream function grid</td><td>N<sub>coarse</sub> & N<sub>ψ</sub></td><td>8192</td><td>12288</td></tr><tr><td>Fine & sample grid</td><td>N<sub>fine</sub> & N<sub>ω</sub></td><td>12288</td><td>24576</td></tr><tr><td>Time-step</td><td>Δt<sub>fluid</sub></td><td>1/24576</td><td>1/36864</td></tr><tr><td>Incompressibility threshold</td><td>δ<sub>inc,b</sub></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Local stencil size for GALS method</td><td>ε<sub>m</sub></td><td>10<sup>-4</sup></td><td>5 · 10<sup>-5</sup></td></tr><tr><td>Low-pass filter size</td><td>k<sub>LP</sub></td><td>4096</td><td>12288</td></tr><tr><td>Fluid time integration scheme</td><td></td><td>RK4Mod</td><td>RK3Mod</td></tr><tr><td>GALS map update stencil</td><td></td><td>4th order</td><td>4th order</td></tr><tr><td>Number of particles</td><td>N<sub>P</sub></td><td>10<sup>6</sup></td><td>10<sup>7</sup></td></tr><tr><td>Particle time integration scheme</td><td></td><td>RK4Mod</td><td>RK3Mod</td></tr></tbody></table>

vortex sheet thickness, the flow is completely dominated by the emerging KH-vortices before any roll-up process starts to occur. In total, simulations for six different $\delta$-values have been performed successfully, ranging from $\delta \approx 0.007 = 1/\sqrt{2 \cdot 10^4}$ to $\delta \approx 0.045 = 1/\sqrt{5 \cdot 10^2}$. The values were placed with $1/\delta^2$ in order to ensure comparability to the Navier-Stokes results from Caflisch et al.⁶ by their overall dynamics, as they did for their inviscid results.

All runs were performed until the map-stack of emerging sub-maps depleted the CPU RAM of the available compute nodes. The computations were computed on the IDRIS supercomputer. Up to 43 convoluted maps for A100 runs and 117 convoluted maps for V100 runs were captured. Details on the executed run architecture and final observed time together with the amount of captured windings of the emerging vortex structures are given in table II. The settings ensure that energy and enstrophy are conserved to an error of $\approx 2 \cdot 10^{-5}$, see figure 7, which is supported by the error of the CM-method being not dissipative in nature for transported quantities.

Simulations with successively decreasing thickness values $\delta$ reveal the flow behavior for the limit $\delta \to 0$. An example for the vortex dynamics for different times and with $\delta \approx 0.032$ is shown in figure 4. Interestingly, two small vortices start to roll up along a common center each. These have a compressed, round-shaped vortex core of high vorticity and two elongated spiral arms. Both vortex cores, being of equal vorticity sign, attract each other and will for long times roll-up into one another. With decreasing sheet thickness, the initial roll-up process is triggered earlier, and the emerging structures are smaller in scale. Figure 5a shows that the vortex core is of slight elliptic shape, being stretched towards the second
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_8.png" /><figcaption>FIG. 7: Relative energy and enstrophy error to the initial condition $t = 0$ for simulations with different $\delta$-values over time. Both quantities are well conserved up to $2 \cdot 10^{-5}$.</figcaption></figure>

vortex due to the mutual attraction. This shape looks similar to those obtained by Baker and Shelley⁷ using layer interfaces for thin vortex layers with the Birkhoff-Rott model. The local palinstrophy captured in figure 5b, is very pronounced in the spiral arms and shows the spiralling structures also present within the vortex core itself.

The acquired results are comparable with inviscid flow fields reported by Caflisch et al.¹² and Caflisch et al.⁶. The depicted vortex in figure 6a is in excellent agreement with figure 14d of Caflisch et al.⁶, given in figure 6b. It is important to mention that the original results are depicted in non-equal aspect ratio, due to which the vortex core appears more elliptic.

<table><caption>TABLE II: Overview on the used vorticity layer thickness values $\delta$ including the obtained final time and the number of windings, for which computations were performed.</caption><thead><tr><th>δ ≈</th><th>1/δ²</th><th>Run on</th><th>Final time</th><th>Captured windings</th></tr></thead><tbody><tr><td>0.0045</td><td>5 × 10<sup>4</sup></td><td>A100</td><td>1.24</td><td>0</td></tr><tr><td>0.007</td><td>2 × 10<sup>4</sup></td><td>A100</td><td>2.29</td><td>3</td></tr><tr><td>0.01</td><td>1 × 10<sup>4</sup></td><td>V100</td><td>2.86</td><td>5</td></tr><tr><td>0.0141</td><td>5 × 10<sup>3</sup></td><td>V100</td><td>3.46</td><td>6</td></tr><tr><td>0.022</td><td>2 × 10<sup>3</sup></td><td>V100</td><td>4.37</td><td>6</td></tr><tr><td>0.032</td><td>1 × 10<sup>3</sup></td><td>V100</td><td>5.82</td><td>7</td></tr><tr><td>0.045</td><td>5 × 10<sup>2</sup></td><td>V100</td><td>7.59</td><td>7</td></tr></tbody></table>
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_9.png" /><figcaption>FIG. 8: Initial palinstrophy growth shows a $t^2$-scaling (a). At later times an exponential palinstrophy growth is found which is due to the formation of rolled-up vortices (b). The steepness of the slope increases for smaller $\delta$-values and exponential growth starts at earlier times.</figcaption></figure>

V. ANALYSIS OF THE RESULTS

In the following, the numerical results are analyzed and discussed. Comparison of global quantities as well as the material line uncover the different flow dynamics. Normalization in space and time with respective scaling laws give insight into dynamics for the non-smooth initial condition when $\delta \to 0$.

The palinstrophy $P$ over time, shown in figure 8, exhibits two different growth stages. At early times, all performed simulations show an algebraic growth following a $t^2$-scaling. This initial growth period is later on overtaken by an exponential growth once vortex structures emerge. The steepness of the growth increases with decreasing vortex sheet thickness $\delta$. For the slope $s$ of the exponential growth with $P(t) \propto \exp(st)$, an empirical relation of $s \approx \delta^{-0.77}$ was found with linear fitting. This means that in the limit of $\delta \to 0$ palinstrophy diverges. This is consistent with results in the literature for Navier-Stokes where exponential growth of palinstrophy was derived, see e.g. Lesieur$^{35}$. Early in the evolution a power-law behavior was predicted in Ayala and Protas$^{32}$ and estimates of the maximum palinstrophy growth were given. The lower limit of the investigated $\delta$-values, i.e. $\delta \approx 0.0045$, experiences artifacts in the form of numerically accelerated secondary Kelvin-Helmholtz instabilities along the
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_10.png" /><figcaption>FIG. 9: Arc length of the material curve for $\delta = 0.032$ with vortex center as vertical dotted line (a). Curvature $\kappa_c$ (b) and true vortex strength $\gamma_c$ (c) over arc length $s(\theta)$ around vortex center for $\delta = 0.014$.</figcaption></figure>

material line, resulting in the premature and irregular steep palinstrophy growth.

With the material line $x(\theta, t)$, given by the individual particle positions, we define following⁶ several quantities, the arc length $s(\theta)$, curvature $\kappa_c(\theta)$ and true vortex strength $\gamma_c(\theta)$. These are used to analyze the material line dynamics and are defined respectively as,

$$s(\theta) = \int_{0}^{2\pi} |x_{\theta}| d\theta, \qquad (16)$$

$$\kappa_c(\theta) = (x_{\theta}y_{\theta\theta} - y_{\theta}x_{\theta\theta})/(x_{\theta}^2 + y_{\theta}^2)^{3/2}, \qquad (17)$$

$$\gamma_c(\theta) = |x_{\theta}|^{-1}. \qquad (18)$$

with the notation $x_{\theta} = \frac{\partial x}{\partial \theta}$ and correspondingly for the other derivatives. Numerically the derivatives were computed using central 4th order finite differences over the particle positions for both the first and second derivative. Both the elongation of the spiral arms and contraction around the vortex core from figure 5 can be found in the evolution of the arc length of the material line in figure 9a. The overall arc length is both increased by the emerging two vortices and the global vortex merging process with the effect intensifying strongly over time especially around the vortex core position. The almost vertical lines correspond to strong stretching of the material line, while the almost horizontal part corresponds to extreme compression of the particles. Our results for the Euler case are comparable with those shown in Caflisch et al.⁶ figure 12b for Navier-Stokes.

The curvature shown in figure 9b showcases the strength of the winding process, forming
Singularity Formation in two-dimensional Euler Equations

two peaks of opposite sign around the vortex cores. These are situated at the edge of the coa-
lesced vortex centers from which at longer times the spiral arms start. Involving second order
derivatives the curvature experiences high-frequency oscillations from numerical derivation
of the particle positions. A Gaussian filter with standard deviation of 3 · 10⁻⁴ × Nₚ was
used to mitigate this effect. In the center of the developed vortex the true vortex strength
reaches a peak showcasing strong convergence towards the vortex core. Before any spiral
arms are formed all particles are compressed towards the vortex core (with γₑ > 1). Once
the vortices start to rotate, two spiral arms will form and elongation occurs (γₑ < 1) outside
the vortex core. Both quantities reach maximum peak values at different times, depicted by
the green curve in figure 9c, and then start to disperse.

Tracking the maximum value of curvature and vortex strength over time enables compar-
isons between different δ-values, shown in figure 10. Here, the curves have been matched
empirically to the given scaling laws. The curvature was found to empirically scale with
δ⁻⁰.⁹ and the true vortex strength with δ⁻⁰.³¹, for which an explanation was not found.
Nonetheless, this scaling supports the conjecture of Caflisch et al.⁶ that only for the limit
of δ → 0 the curvature as well as true vortex strength go towards infinity for finite time.
Interestingly, over time, both quantities show the same behavior with an initial steep increase

<figure><img src="image_11.png" /><figcaption>FIG. 10: Temporal re-scaling for the maximum values of the curvature κs and true vortex strength γs as well as temporal re-scaling for the vortex turnover with empirically determined coefficients. The critical times ts, given in figure (b), were found to be decreasing with lower δ-values.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 11: Temporal re-scaling for the vortex turnover (a) and palinstrophy growth (b) with empirically determined coefficients. The critical times $t_s$ are given in figure (a).</figcaption></figure>

to an overall maximum, decreasing again with oscillations. As the observed dispersion effect appears similar for different vortex sheet thicknesses $\delta$ it is interpreted as a flow behavior rather than a result of numerical artifacts. It can be explained by the flow evolution. At first the vorticity sheet close to the vortex core starts to roll up. An inner and an outer region are formed based on the true vortex strength: one close to the center of the vortex where the material line with surrounding region of vorticity is condensed and one outside where the material line is strongly elongated from the formation of the spiral arms due to the rotation of the vortex core. At later stages the vortex core is reformed into a circular shape (observable in figure 5a) and rotates with quasi-constant speed due to the induced velocity from the Biot-Savart law. This shape bears similarities to what Baker and Shelley$^7$ observed for vortex layer interfaces.

This is also visible in figure 11a. Here the occurrence of the vanishing gradient is plotted as a function of a re-scale time, where each occurrence marks when the gradient of the material line at the spiral center becomes zero, i.e. when either the x- or y-component of the derivative of the material line position $\frac{\partial x(\theta, t)}{\partial \theta}$ vanishes. As the vortex turnover over re-scaled time occurs in constant periods (linear growth in figure 11a), the vortex rotates with constant speed. It is similar to a solid body rotation and portrays a stable vortex core.

Additionally, the palinstrophy growth was also found to be re-scalable in time using $t_s$ and $\delta$, as shown in figure 11b.
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_13.png" /><figcaption>FIG. 12: Spatial re-scaling with different factors for the fourth time that $\partial x(\theta, t)/\partial\theta \approx 0$ in the $y$ direction for different $\delta$-values.</figcaption></figure>

In all graphs of figures 10 and 11, a temporal scaling has been used. As for different $\delta$-values the results are similar, those four criteria were used to match them and determine an empirical relation. The ansatz was taken from Caflisch et al.⁶ using $(t-t_s)\delta^{-1}$, where $t_s$ is the singularity time and $\delta^{-1}$ as re-scaling factor. The flow scales linearly with decreasing $\delta$-values. The critical times are not equal to $t_s = 1.505$, the value of the Birkhoff-Rott equation but were observed to be slightly larger. They are closely constant for larger $\delta$-values but changed once it was decreased, possibly going towards the critical time of the BR-equation for vanishing vortex sheet thickness $\delta$. The linear scaling with $\delta^{-1}$ in comparison to the scaling of $\delta^{-2/3}$ for viscous flow reported by Caflisch et al.⁶ results from different flow dynamics. Results by Baker and Pham⁸ for vortex blob methods for Euler flow report on the same re-scaling factor on time. Without viscosity to diffuse the flow, the initial vorticity field can only be transported and the whole dynamics of the condensed vortex core is an accumulation of vorticity present in the initial vortex sheets, which scales directly with the thickness $\delta$.

The material curves for different $\delta$-values have strong spatial self-similarity. Results obtained by Caflisch et al.⁶ for Navier-Stokes simulations and Baker and Pham⁸ for vortex blob methods show a scaling of $\delta^{-1}$, however, this did not match our computed results.
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_14.png" /><figcaption>FIG. 13: Energy spectra $E(k)$ defined in equation 4 for different $\delta$-values. Both start with an initial scaling of $k^{-2}$ with exponential decay for large $k$ as given in figure 1b.</figcaption></figure>

These show an empirical scaling of approximately $\delta^{-0.9}$, which is slightly different, as shown in figure 12. This similarity scaling is valid for all reported simulations and all times around the vortex cores and only breaks down with the vortex merging process for later times (observed for $\delta \approx 0.045$ for $t > 7$). There has been no reason found for the discrepancies to other literature, however differences might arise from the formation of the outer and inner region, which differs in shape to Baker and Pham$^8$ and might behave different to those of the viscous Navier-Stokes results of Caflisch et al.$^6$.

With the scaling relation of the thickness of the vortex core region, it reduces to a point size for the limit $\delta \to 0$. The point of maximum curvature, being at the interface between the core region and the start of the spiral arms, for the limits falls into the center point where the diverging true vortex strength compresses the vortex line into an individual point.

The characteristic of the initial energy spectra, defined in equation 4, have already been explained in section II. For larger $\delta$-values the initial profile with exponential decay experiences a shift of energy to finer scales (figure 13a). In comparison, the results with lower $\delta$-values with initially more pronounced $k^{-2}$ Dirac-scaling show only little shift in energy from larger to finer scales (figure 13b). However, for both simulations the two vortices have already emerged for the captured time and performed several windings, while the results for figure 13b did not yet observe any developed distortion from the vortex merging process.
# Singularity Formation in two-dimensional Euler Equations

Eventually, once this global process continues the energy spectra are expected to behave similarly as to that for larger $\delta$-values. The formation of the local vortices exhibits therefore only little impact on the energy spectra.

The two vortices which form are quite unstable once perturbed to sufficient degree. With ongoing global filamentation process of the two vortices, the shape of the vortex core becomes more and more elliptic (figure 14) towards each other. This brings an imbalance to the induced velocity field and secondary vortices form, breaking down the regularity. The individual parts of the spiral arm become compressed into filament structures and secondary vortices form.

<figure><img src="image_15.png" /><figcaption>FIG. 14: Material line (in black) and vorticity for $\delta = 0.032$ at different times. The initially more stable spiral arms experience distortions and filament structures start to emerge with ongoing global vortex merger process.</figcaption></figure>

# VI. SINGULARITY ANALYSIS

The complex singularity of an analytic function can be analyzed using the Fourier transform. The width of the analyticity strip can be obtained by considering the asymptotic behavior of the Fourier spectrum governed by Laplace's formula. To obtain information about singularities outside the analyticity strip the Borel-Polya-van der Hoeven (BPH) method has been proposed³⁶.
# Singularity Formation in two-dimensional Euler Equations

## A. Borel-Polya-van der Hoeven method

The Borel-Polya-van der Hoeven (BPH) method³⁶ is a numerical tool for finding complex singularities of single variable functions by combining the information on singularities obtained from Borel transforms for Taylor series³⁷ with numerical techniques for asymptotic interpolation³⁸. This method has been successfully applied to find singularities for the 1D Burgers equations³⁶ and has been used to analyze potential singularity formation in the 2D Euler vortex sheet problem⁶. The following is a summary of the method along with some small modifications implemented to adapt to our case.

Given a complex-valued function $f(Z)$ with $Z \in \mathbb{C}$ and with formal power series

$$f(Z) = \sum_{n=0}^{\infty} a_n Z^n, \qquad (19)$$

its Borel transform and Borel-Laplace transform are given respectively by

$$f^B(\xi) = \sum_{n=0}^{\infty} \frac{a_n}{n!} \xi^n \qquad (20a)$$

$$f^{BL}(Z) = \frac{1}{Z} f\left(\frac{1}{Z}\right) = \sum_{n=0}^{\infty} \frac{a_n}{Z^{n+1}}, \qquad (20b)$$

so named since $f^{BL}$ is formally the Laplace transform of $f^B$.

The BPH method is built on Pólya's theorem which is based on the observation that for $a_n = c^n$ for some complex number $c$, the Borel transform is the exponential

$$f^B(\xi) = \sum_{n=0}^{\infty} \frac{(c\xi)^n}{n!} = e^{c\xi}, \qquad (21)$$

while the Borel-Laplace transform is a simple pole at $c$:

$$f^{BL}(Z) = \sum_{n=0}^{\infty} \frac{c^n}{Z^{n+1}} = \frac{1}{Z-c}. \qquad (22)$$

Writing $f^B$ in polar coordinates using $\xi = re^{-i\theta}$ (reversed phase parametrization for convenience), and $c = |c|e^{i\phi}$, we have that

$$\ln(f^B(r, \theta)) = |c| r(\cos(\phi - \theta) + i \sin(\phi - \theta)), \qquad (23)$$

and therefore $\partial_r \ln(|f^B(r, \theta)|)$ is maximized at $\theta = \phi$ to value $|c|$ thereby revealing the position of the pole.
## Singularity Formation in two-dimensional Euler Equations

For $f^{BL}$ given by a linear combination of multiple isolated poles $f^{BL}(Z) = \frac{C_1}{Z-c_1} + \frac{C_2}{Z-c_2} + \cdots + \frac{C_m}{Z-c_m}$, the supporting function $\sigma(\theta)$ is given by

$$ \sigma(\theta) := \lim_{r \to \infty} \partial_r \ln(|f^B(r, \theta)|) \\ = \max_{j=1,2,\dots,m} |c_j| \cos(\phi_j - \theta). \qquad (24) $$

The curve $\sigma(\theta)e^{i\theta}$ then describes the convex hull of the set of singularities. This also means that only the singularities of $f^{BL}$ furthest from the origin are identified from the vertices of the convex hull. Since the Borel-Laplace transform $f \mapsto f^{BL}$ involves the change of variable $Z \mapsto 1/Z$, we have that the complex reciprocals of these singularities are the singularities of $f$ closest to the origin.

This method is then applied to estimate the location of complex plane singularities for a real valued periodic function $u(x)$ on $S^1 \sim [-\pi, \pi)$. We write $u$ in terms of Fourier series and make the analytic extension in some strip around the real axis

$$ u(z) = \sum_{k \in \mathbb{Z}} \hat{u}_k e^{ikz} \qquad (25) $$

which is decomposed as the sum of two functions

$$ u^{+}(z) = \sum_{k>0} \hat{u}_k e^{ikz}, \qquad (26a) $$

$$ u^{-}(z) = \sum_{k>0} \hat{u}_{-k} e^{-ikz}. \qquad (26b) $$

We note that $u^+$ is analytic in the upper half plane, i.e. poles are contained in the lower half plane. Similarly, all poles of $u^-$ are in the upper half plane. The change of variables $Z = e^{iz}$ and $Z = e^{-iz}$ are made to $u^+$ and $u^-$ respectively to yield

$$ f_1(Z) = u^{+}(-i \ln Z), \qquad (27a) $$

$$ f_2(Z) = u^{-}(i \ln Z), \qquad (27b) $$

where we use the $[-\pi, \pi)$ branch cut. We note now that both $f_1$ and $f_2$ are analytic inside the unit disk corresponding to the image of the upper and respectively lower half planes under the change of coordinate. Applying the BPH method on $f_1$ and $f_2$ then reveals the poles closest to the unit circle and hence closest to the real axis after the logarithmic change of coordinate.
## Singularity Formation in two-dimensional Euler Equations

The algorithm of the BPH method is summarized as follows. For a function $f$ given as a truncated power series in (19) with coefficients $a_n$, we compute its Borel transform $f^B(\xi)$ on a grid in polar coordinates given by $\xi_{j,k} = r_j e^{-i\theta_k}$. The grid points are given by

$$r_j = r_0 + \frac{jR}{M}, \qquad \theta_k = \frac{2\pi k}{K} \qquad (28)$$

for $j = 0, 1, 2, \dots, M-1$ and $k = 0, 1, 2, \dots, K-1$ with an appropriate choice of $r_0$ and $R$ to be discussed later.

At each grid point, the Borel transform is given by a $N$-truncated sum of (20). If we select the number of rays $K$ to be equal to the number of power series terms $N$, we get that

$$f^B(\xi_{j,k}) = \sum_{n=0}^{N-1} a_n \frac{r_j^n}{n!} e^{-i \frac{2\pi k n}{N}} = \sum_{n=0}^{N-1} a_n \frac{r_j^n}{n!} e^{-i \theta_k n}, \qquad (29)$$

which, for fixed $j$, is the expression for the discrete Fourier transforms of the sequence $\{b_n\}_{n=0}^{N-1}$ with

$$\begin{aligned} b_n &= a_n \frac{r_j^n}{n!} \\ &= a_n \left( r_j \exp \left( \frac{-\ln \Gamma(n+1)}{n} \right) \right)^n \\ &= \exp \left( n \ln r_j + \ln a_n - \ln \Gamma(n+1) \right). \end{aligned} \qquad (30)$$

FFT algorithms are used to compute the discrete Fourier transform in order to reduce round-off errors and improve speed and the log-Gamma function is used to avoid numerical overflow in the computation of the larger $n$ terms.

Then for each fixed $\theta_k$, $\sigma(\theta_k)$ can be estimated from finite difference on $\ln(|f^B(r, \theta)|)$. Higher order asymptotic interpolation could be used but would require high precision arithmetic. For double precision computations, we opted for simple finite difference. From the supporting function $\sigma(\theta)$, we can identify $S^{BL} = \{y_1^*, y_2^*, \dots, y_s^*\}$, the vertices of the convex hull of the singularities of $f^{BL}$ which we recall lie within the unit complex disk. The choice of the values of $r_0$ and $R$ needs to be adjusted carefully. Since along the ray of each $y_i^*$, the Borel transform is expected to grow as $f^B(\xi) \sim \exp(y_i^* \xi)$, in order to capture the exponential growth, one should pick $r_0$ larger than $1/|y_i^*|$. The maximum radius $R$ should be picked as large as possible until the effects of finite precision arithmetic start dominating.

Applying the inverse Borel-Laplace transform gives us the set of singularities $S = \{Z_1^*, Z_2^*, \dots, Z_s^*\}$ for $f$, where $Z_j^* = 1/y_j^*$, which are now the singularities of $f$ closest to
Singularity Formation in two-dimensional Euler Equations

<figure><img src="image_16.png" /><figcaption>FIG. 15: Evolution of the complex plane singularities for the true vortex strength (a) and curvature (b) for various values of $\delta$ from $t = 0.8$ to 1.6. The position of the singularity at time $t = 1.5$ is marked by the circles. The singularity positions are computed using the BPH method using radii in [1, 2.4]. The use of larger radii was not possible due to the presence of noise in the data.</figcaption></figure>

the unit circle. The locations of the singularities $z_i^*$ of $u^\pm$ are then obtained by performing the inverse transform $z_i^* = \mp i \ln Z$. This means that if one has an a priori estimate for the radius of analyticity $\delta$ of $u$, the value $r_0$ should be chosen larger than $e^\delta$.

We applied the BPH method to analyze the singularity formation in the vortex core curve $\gamma$. As in$^6$, we compute the closest complex singularities for the vortex strength function $\sigma$ and the curvature function $\kappa$ given by

$$ \sigma = |\dot{\gamma}|^{-1}, \qquad (31) $$

$$ \kappa = \frac{|\dot{\gamma} \times \ddot{\gamma}|}{|\dot{\gamma}|^3}. \qquad (32) $$

Since the CMM relies on Hermite cubic interpolation for the stream function which yields a $C^0$ velocity field with small discontinuities in the derivative (order of $h^3$ where here $h \sim 1e-9$), the evolved vortex curve $\gamma$ is also only $C^0$. The lack of regularity shows up as noise in the large wave-numbers of $\hat{\gamma}$. To perform the BPH analysis on this data, a smoothing kernel $\exp(-1 \times 10^{-7}k^4)$ in Fourier space and a truncation to the 10000 first wave-numbers is applied to the curve as pre-processing. The results are shown in figure 15, where the position of the complex singularities are drawn for times between 0.8 and 1.6. The general trajectories of the complex plane singularities before time $t = 1.5$ are in agreement with
# Singularity Formation in two-dimensional Euler Equations

the results found for the viscous case in⁶. We observe that all singularities move towards the real axis as time advances until $t = 1.5$ (marked with circles) after which they seem to move away in some cases. The meaning of this is not clear from these tests as the accuracy of the singularity analysis is very limited due to the limited available arithmetic precision. We cannot exclude the possibility that the BPH method is detecting another complex singularity which overtakes the one tracked in the above figures resulting in a transition to a new singularity which is closest to the real axis.

## VII. CONCLUSIONS

The flow of vortex layers governed by the incompressible 2D Euler equations has been computed for successively decreasing layer thickness using the characteristic mapping method. This semi-Lagrangian method features exponential resolution in linear time and thus allows to capture the exponential growth of the vorticity gradients. Our results agree with pseudo-spectral computations of Caflisch et al.⁶ without viscosity and go even beyond their reported results. In particular, the range of results for vanishing vortex sheets thickness was extended down to $\delta \approx 0.007$ and longer captured simulations enabled further analysis of the emerging vortex structures. Energy and enstrophy are conserved to a high degree thanks to the non-dissipative feature of the CM method. The palinstrophy shows at the beginning a super-exponential growth, which is intercepted and dominated by a stronger exponential growth once the two vortices form. The dynamics of the center-line of the vortex sheet can be described by measures of the arc length, curvature and vortex strength. The arc length and vortex strength showcase, that the emerging vortices develop two individual regions: one core region of very strong compression around the vortex centers and one of elongation in the spiral arms that form. The core region forms a peak in compression directly at the center, which increases with decreasing vortex sheet thickness $\delta$. For the limit of $\delta \to 0$ the vortex strength is expected to form a singularity at the center of rotation where it is compressed into it with vortex strength tending to infinity. The curvature displays similar singular behavior. Instead of forming maximum values directly at the vortex center, it forms two at the interface between the inner and outer region of the condensed vortex core. These are of opposite sign and describe the transition of the center line from the vortex blob formation to the spiral arms and additionally scale with
Singularity Formation in two-dimensional Euler Equations

vortex sheet thickness $\delta$, again going to infinity for the limit of vanishing $\delta$-values, further supporting the formation of singularity for non-smooth initial data suggested in Caflisch *et al.*$^6$. The investigated quantities show strong self-similarity in time and space. With the help of the curvature, vortex strength and turnover time, the re-scaling in time of the vortex dynamics was unveiled. In fact, after initial build up these vortices rotate with constant speed, similar to a solid body rotation. The maximum values for the curvature and vortex strength also do not observe monotonous growth over time, but after a steep increase, they oscillate around a common value. The temporal re-scaling was found to be consistent with other vortex sheets results of incompressible Euler equations$^8$. In space, the formed vortices show strong self-similarity even for many turnover events and over a large range of $\delta$-values. The observed scaling slightly differs from that found in Caflisch *et al.*$^6$ in the viscous case, i.e. for Navier-Stokes or$^8$ for vortex blob methods. The energy spectra show for small $\delta$ values a power law scaling with slope close to -2. Simulations for longer times show, that at some point the vortex merger starts to distort the round-shaped vortex blobs, this leads to instability in the roll-up process and secondary vortices start to form, eventually breaking down the spiral structure by continuous filamentation.

In future work we will apply CMM to compute fine scale structures of vorticity gradients in 2D Euler. The transport equation of the vorticity gradients contains a source term for vorticity gradient stretching, similar to the vortex stretching source term in 3D Euler. The latter can be solved likewise by CMM, similar to what has proposed in Yin, Schneider, and Nave$^{30}$.

Another challenging perspective is applying CMM to study Euler flows in 3D and to investigate numerically possible singularities. First results presenting low resolution computations for the Hou and Li$^{22}$, Kerr$^{39}$ initial condition can be found in Yin, Schneider, and Nave$^{30}$. High resolution 3D computations considering different flow configurations, like Hou and Li$^{22}$, Kerr$^{39}$ and more recently the one by Moffatt and Kimura$^{40,41}$, will be published in forthcoming work.
Singularity Formation in two-dimensional Euler Equations

## ACKNOWLEDGEMENTS

The authors acknowledge partial funding from the Agence Nationale de la Recherche (ANR), grant ANR-20-CE46-0010-01. This work was granted access to the HPC resources of IDRIS under the allocation 2022-91664 made by GENCI. Center de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. Part of this research was finalized while KS was visiting the Institute for Mathematical and Statistical Innovation (IMSI) at University of Chicago, which is supported by the National Science Foundation (Grant No. DMS-1929348). JCN Acknowledges partial support from the NSERC Discovery Program.

## Appendix

In the following we present numerical validation of the open access CMM cuda code³³ similar to what has been done in Yin et al.²⁹ for the Matlab implementation.

### Validation of convergence order in space and time for Cuda-code³³

The convergence tests presented in the following we use the 4-mode-flow with the same parameters as in Yin et al.²⁹ in order to make the results comparable. The initial condition is given by,

$$ \omega_0(x, y) = \cos(x) + \cos(y) + 0.6 \cos(2x) + 0.2 \cos(3x). \quad (33) $$

which develops quickly small scale features and is thus a severe test case for any numerical method, see e.g. Podvigina, Zheligovsky, and Frisch²⁴, Yin et al.²⁹. In table III, all numerical parameters for the reference simulations can be found for which the error analyses have been carried out. The grid for the flow map $\chi(N_{coarse})$, initial vorticity $\omega_0(N_{fine})$ and sampling of the vorticity for the FFT ($N_\omega$) were chosen of same size and the velocity field uses a grid with increased values to reduce the influence of the non-smoothness of the velocity-gradients. The time-step $\Delta t$ is set to a Courant-Friedrich-Lewis (CFL) number of 2. The size $\epsilon_m$ defines the stencil size for the GALS method. In the presented computational study no low-pass filtering and no remapping was used in order to capture the sub-map error
## Singularity Formation in two-dimensional Euler Equations

correctly. All simulations were run until a final time of $t = 1$. For further understanding of the impact of parameters, the readers are directed to$^{29,34}$.

<table><caption>TABLE III: Settings of the reference simulation$^{29}$.</caption><thead><tr><th>Name</th><th>Value</th><th>Name</th><th>Value</th></tr></thead><tbody><tr><td>$N_{\text{coarse}}$</td><td>1024</td><td>$N_{\text{fine}}$</td><td>1024</td></tr><tr><td>$N_{\psi}$</td><td>2048</td><td>$N_{\omega}$</td><td>1024</td></tr><tr><td>$h_{\text{fluid}}$</td><td>1/512</td><td></td><td></td></tr><tr><td>$\epsilon_m$</td><td>10<sup>-3</sup></td><td>Initial condition</td><td>4-mode-flow</td></tr><tr><td>Fluid time scheme</td><td>RK3</td><td>Map update stencil</td><td>4th order</td></tr></tbody></table>

The errors of four quantities were computed to examine the convergence order. Those are the flow map and vorticity error in $L_{\infty}$-norm and the energy and entrophy conservation error in $L_2$-norm. All the errors were evaluated by sampling the map on a uniform 2048$^2$-grid and computing the vorticity and velocity respectively on this map.

$$ \text{Map error} = \|\chi_{\text{ref}}(\cdot, t_n) - \chi(\cdot, t_n)\|_{\infty} \quad (34) $$

$$ \text{Vorticity error} = \|\omega_{\text{ref}}(\cdot, t_n) - \omega(\cdot, t_n)\|_{\infty} \quad (35) $$

$$ \text{Energy error} = \|\omega(\cdot, t_n)^2\|_2 - \|\omega(\cdot, t_0)^2\|_2 \quad (36) $$

$$ \text{Enstrophy error} = \|\mathbf{u}(\cdot, t_n)^2\|_2 - \|\mathbf{u}(\cdot, t_0)^2\|_2 \quad (37) $$

According to Yin et al.$^{29}$, an error bound for the characteristic map is given by

$$ \tilde{\mathcal{E}}^n = \mathcal{O} (\Delta x^2 \min(\Delta t, \Delta x^2 \Delta t^{-1}) + \Delta t^s + \Delta t^p) \quad (38) $$

Here $s$ and $p$ are the orders of the s-stage Runge-Kutta scheme and the used order of Lagrange time interpolation for the velocity, respectively. Values were set as $p = s$ to balance computation requirements with achieved numerical accuracy. Due to the bi-cubic spatial Hermite interpolation used, errors of order $\mathcal{O}(\Delta x^3)$ for the convergence order in space are expected (figure 16). One order is reduced to the theoretical 4th order as the velocity to advect the map using the GALS framework is sampled as a first order derivative from the stream function. The CUDA code is capable of including first to fourth order Runge-Kutta scheme. An example for convergence in time with third order scheme paired with third order Lagrange-interpolation of the velocity is depicted in figure 16. Again, all quantities converge with the expected order. The evolution of fluid particles used in the manuscript
Singularity Formation in two-dimensional Euler Equations

to track the material line of the vortex sheet, is computed using Lagrangian point particles, where the velocity is set to that of the fluid flow,

$$ \frac{d\boldsymbol{x}_p}{dt} = \boldsymbol{u} \qquad (39) $$

Here $\boldsymbol{x}_p$ is the position of the point particle and $\boldsymbol{u}$ the fluid velocity. For time integration we tested different numerical schemes (cf. figure 16c), for spatial interpolation of the velocity at the particle positions bi-cubic Hermite interpolation is applied. Similarly to the map error quantities, the particle position error can be defined as

$$ \text{Particle error} = \|\boldsymbol{x}_{p,ref}(\cdot, t_n) - \boldsymbol{x}_p(\cdot, t_n)\|_{\infty} \qquad (40) $$

The particles were initially scattered with uniform random distribution over the computational domain. Due to the volume preservation property of the incompressible flow, their distribution will remain uniform as well. All fluid parameters were kept similarly to the reference computation (table III) and 10⁶ particles were deployed. The tested time-stepping methods show the expected convergence orders. Merely fourth order deployed time-stepping schemes are reduced to third order, which is due to the third order convergence of the fluid velocity (figure 16c).

<figure><img src="image_17.png" /><figcaption>FIG. 16: Convergence errors from equation 34-37 in space and time for the flow variables and equation 40 for the fluid particles. Expected convergence rates are plotted with dotted lines.</figcaption></figure>
Singularity Formation in two-dimensional Euler Equations

REFERENCES

¹V. Yudovich, “Non-stationary flow of an ideal incompressible liquid,” USSR Comput. Math. & Math. Phys. **3**, 1407–1456 (1963).

²L. Székelyhidi, “Weak solutions to the incompressible Euler equations with vortex sheet initial data,” C. R. Math. **349**, 1063–1066 (2011).

³A. Majda and A. Bertozzi, *Vorticity and Incompressible Flow*, Cambridge Texts in Applied Mathematics (Cambridge University Press, 2001).

⁴C. Bardos and E. Titi, “Mathematics and turbulence: where do we stand?” J. of Turb. **14**, 42–76 (2013), https://doi.org/10.1080/14685248.2013.771838.

⁵P. G. Saffman, *Vortex dynamics* (Cambridge university press, 1995).

⁶R. Caflisch, F. Gargano, M. Sammartino, and V. Sciacca, “Complex singularity analysis for vortex layer flows,” J. Fluid Mech. **932**, A21 (2022).

⁷G. Baker and M. Shelley, “On the connection between thin vortex layers and vortex sheets,” J. Fluid Mech. **215**, 161–194 (1990).

⁸G. Baker and L. Pham, “A comparison of blob methods for vortex sheet roll-up,” J. Fluid Mech. **547**, 297–316 (2006).

⁹S.-I. Sohn, “Singularity formation and nonlinear evolution of a viscous vortex sheet model,” Phys. Fluids **25**, 014106 (2013), https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/1.4789460/13470616/014106_1_online.pdf.

¹⁰S.-I. Sohn, “Self-similar roll-up of a vortex sheet driven by a shear flow: Hyperbolic double spiral,” Phys. Fluids **28**, 064104 (2016), https://pubs.aip.org/aip/pof/article-pdf/doi/10.1063/1.4953780/15833154/064104_1_online.pdf.

¹¹D. Pullin and N. Shen, “On vortex-sheet evolution beyond singularity formation,” J. Fluid Mech. **976**, A17 (2023).

¹²R. Caflisch, F. Gargano, M. Sammartino, and V. Sciacca, “Regularized Euler-α motion of an infinite array of vortex sheets,” Bollettino Unione Mat. Ital. **10**, 113–141 (2017).

¹³G. Birkhoff, “Helmholtz and Taylor instability,” in *Proceedings of Symposia in Applied Mathematics*, Vol. XII (Amer. Math. Soc., 1962) pp. 55–76.

¹⁴N. Rott, “Diffraction of a weak shock with vortex generation,” J. Fluid Mech. **1**, 111–128 (1956).

¹⁵D. Moore, “The spontaneous appearance of a singularity in the shape of an evolving vortex
Singularity Formation in two-dimensional Euler Equations

sheet,” Proc. R. Soc. A: Math. Phys. Eng. Sci. 365, 105–119 (1979).

<footnote n="16">¹⁶R. Nguyen Van Yen, M. Farge, and K. Schneider, “Energy dissipating structures produced by walls in two-dimensional flows at vanishing viscosity,” Phys. Rev. Lett. 106, 184502 (2011).</footnote>

<footnote n="17">¹⁷N. Nguyen Van Yen, M. Waidmann, R. Klein, M. Farge, and K. Schneider, “Energy dissipation caused by boundary layer instability at vanishing viscosity,” J. Fluid Mech. 849, 676–717 (2018).</footnote>

<footnote n="18">¹⁸H. Clercx and G. Van Heijst, “Dissipation of coherent structures in confined two-dimensional turbulence,” Phys. Fluids 29 (2017).</footnote>

<footnote n="19">¹⁹P. Constantin, M. Lopes Filho, H. J. Nussenzveig Lopes, and V. Vicol, “Vorticity measures and the inviscid limit,” Archive for Rational Mechanics and Analysis 234, 575–593 (2019).</footnote>

<footnote n="20">²⁰M. Hussaini and T. Zang, “Spectral methods in fluid dynamics,” Annu. Rev. Fluid Mech. 19, 339–367 (1987), https://doi.org/10.1146/annurev.fl.19.010187.002011.</footnote>

<footnote n="21">²¹T. Ishihara, T. Gotoh, and Y. Kaneda, “Study of high-Reynolds number isotropic turbulence by direct numerical simulation,” Annu. Rev. Fluid Mech. 41, 165–180 (2009), https://doi.org/10.1146/annurev.fluid.010908.165203.</footnote>

<footnote n="22">²²T. Hou and R. Li, “Computing nearly singular solutions using pseudo-spectral methods,” J. Comp. Phys. 226, 379–397 (2007).</footnote>

<footnote n="23">²³J. Gibbon, M. Bustamante, and R. Kerr, “The three-dimensional Euler equations: singular or non-singular?” Nonlinearity 21, T123 (2008).</footnote>

<footnote n="24">²⁴O. Podvigina, V. Zheligovsky, and U. Frisch, “The Cauchy–Lagrangian method for numerical analysis of Euler flow,” J. Comput. Phys. 306, 320–342 (2016).</footnote>

<footnote n="25">²⁵J.-C. Nave, R. Rosales, and B. Seibold, “A gradient-augmented level set method with an optimally local, coherent advection scheme,” J. Comput. Phys. 229, 3802–3827 (2010).</footnote>

<footnote n="26">²⁶B. Seibold, R. R. Rosales, and J.-C. Nave, “Jet schemes for advection problems,” Discrete and Continuous Dynamical Systems - B 17, 1229–1259 (2012).</footnote>

<footnote n="27">²⁷P. Chidyagwai, J.-C. Nave, R. R. Rosales, and B. Seibold, “A comparative study of the efficiency of jet schemes,” arXiv preprint arXiv:1104.0542 (2011).</footnote>

<footnote n="28">²⁸O. Mercier, X.-Y. Yin, and J.-C. Nave, “The characteristic mapping method for the linear advection of arbitrary sets,” SIAM Journal on Scientific Computing 42, A1663–A1685 (2020).</footnote>

<footnote n="29">²⁹X.-Y. Yin, O. Mercier, B. Yadav, K. Schneider, and J.-C. Nave, “A characteristic mapping</footnote>
Singularity Formation in two-dimensional Euler Equations

method for the two-dimensional incompressible Euler equations," J. Comput. Phys. **424**,
109781 (2021).

<sup>30</sup>X.-Y. Yin, K. Schneider, and J.-C. Nave, “A characteristic mapping method for the three-dimensional incompressible Euler equations,” J. Comput. Phys. **477**, 111876 (2023).

<sup>31</sup>M. Abramowitz and I. Stegun, *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables*, 9th ed. (Dover, New York, 1964).

³²D. Ayala and B. Protas, “Maximum palinstrophy growth in 2d incompressible flows,” J. Fluid Mech. **742**, 340–367 (2014).

³³B. Yadav, T. Maurel Oujia, N. Saber, J. Bergmann, and P. Krah, “CMM cuda code,” https://github.com/CharacteristicMappingMethod/cmm-turbulence (2023).

³⁴J. Bergmann, *Investigation of mixing and particle transport in 2D incompressible Euler flows using the characteristic mapping method*, Master’s thesis, I2M, Aix-Marseille Université (2022).

³⁵M. Lesieur, “Two-dimensional turbulence,” Turbulence in Fluids: Fourth Revised and Enlarged Edition, 311–348 (2008).

³⁶W. Pauls and U. Frisch, “A Borel transform method for locating singularities of Taylor and Fourier series,” J. Stat. Phys. **127**, 1095–1119 (2007).

³⁷G. Pólya, “Untersuchungen über Lücken und Singularitäten von Potenzreihen,” Ann. Math. **34**, 731–777 (1933).

³⁸J. van der Hoeven, “Algorithms for asymptotic interpolation,” J. of Symb. Comp. (2006).

³⁹R. Kerr, “Evidence for a singularity of the three-dimensional, incompressible Euler equations,” Phys. Fluids A: Fluid Dyn. **5**, 1725–1746 (1993), https://pubs.aip.org/aip/pof/article-pdf/5/7/1725/12783847/1725_1_online.pdf.

⁴⁰H. Moffatt and Y. Kimura, “Towards a finite-time singularity of the Navier–Stokes equa-
tions part 1. derivation and analysis of dynamical system,” J. Fluid Mech. **861**, 930–967
(2019).

⁴¹H. Moffatt and Y. Kimura, “Towards a finite-time singularity of the Navier–Stokes equa-
tions. part 2. vortex reconnection and singularity evasion,” J. Fluid Mech. **870**, R1 (2019).