Singularity formation of vortex sheets in 2D Euler equations using the characteristic mapping method

Abstract

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 2D 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 non-uniqueness of weak solutions of 2D Euler for non-smooth initial data are drawn and the presence of flow singularities is revealed tracking them in the complex plane.

1 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 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 (Yudovich, 1963), it was shown by Székelyhidi (2011) 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 to (Majda & Bertozzi, 2001; Bardos & Titi, 2013) and for 2D vortex layers we refer to the detailed introduction of Caflisch et al. (2022).

In a recent work Caflisch et al. (2022) 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. (2022) 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. (2017) and Caflisch et al. (2022). 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. (2022)). The motion of an inviscid vortex sheet, i.e. for vanishing thickness, is governed by the Birkhoff–Rott model equation (Birkhoff, 1962; Rott, 1956) which develops a singularity in finite time starting from smooth initial data (Moore & Stuart, 1979). Perturbations grow due to the Kelvin–Helmholtz instability and the vortex sheet does roll up. Regularized simulations using a vortex blob method have been also performed in Caflisch et al. (2022) 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 et al. (2011, 2018). 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 & Van Heijst (2017) and mathematical analysis of weak solutions of the 2D Navier–Stokes equations in bounded domains, in the vanishing viscosity limit, in Constantin et al. (2019).

State of the art for solving Navier–Stokes or Euler equations numerically with high precision are pseudo-spectral methods (Hussaini & Zang, 1987; Ishihara et al., 2009) which have been also extensively used for investigating nearly singular solutions of the 3D Euler equations (Hou & Li, 2007; Gibbon et al., 2008). A Cauchy–Lagrange method for computing 2D Euler flows was proposed in Podvigina et al. (2016). 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 Characterisitic Mapping Method (CMM). Evolving the flow map with a Gradient Augmented Level Set Method, developed in (Nave et al., 2010; Seibold et al., 2012; Chidyagwai et al., 2011), 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 et al. (2020) and 2D Euler in (Yin et al., 2021), 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 (Yin et al., 2021). More recently an extension to 3D incompressible Euler flows has been proposed in (Yin et al., 2023). 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 revisit the vortex layer computations of Caflisch et al. (2022) in the inviscid case in more depth using CMM and to extend the range to smaller $\delta$-values and longer times. As presented in Caflisch et al. (2017), 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 into the non-uniqueness of week solutions and about possible singularities. 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 2. A short description of the characteristic mapping method for solving the 2D incompressible Euler equations is given in section 3. Section 4 gives an overview on the performed computations and numerical results are then presented in section 5. A singularity analysis in the complex plane is performed in section 6. Finally, conclusions and perspectives for future work are given in section 7.

2 Governing equations and initial condition

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

$$\partial_{t}\omega+({\bm{u}}\cdot\nabla)\omega=0$$ (1)

where the vorticity is defined as $\omega=\nabla\times{\bm{u}}$ and $\bm{u}$ is the incompressible velocity, satisfying $\nabla\cdot\bm{u}=0$. The curl operator is invertible and the velocity can be computed from the vorticity, $\bm{u}=-\nabla\times\Delta^{-1}\omega$ using the Biot–Savart operator. Equation 1 is completed with a suitable initial condition $\omega_{0}({\bm{x}})=\omega({\bm{x}},t=0)$, here a regularized vortex sheet and where $\bm{x}=(x,y)$ with $(x,y)\in\Omega$ as the cardinal directions. All variables have been non-dimensionalized, similar to Caflisch et al. (2022), as follows,

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. (2022) reads,

$$\omega_{0}(x,y)=\frac{1}{\sqrt{2\pi}\delta}\exp\left(-\frac{(y-\pi-\phi(x))^{2}}{2\delta^{2}}\right)$$ (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$. The position of the center line $\phi$ is the curve of maximum vorticity, which oscillates around $y=\pi$. For $\delta\rightarrow 0$ this vortex layer converges towards a vortex sheet used in Caflisch et al. (2022) as initial condition for the Birkhoff–Rott equation.

(a) Initial condition $\omega_{0}$ for $\delta\approx 0.032$

(b) Energy spectrum for different $\delta$-values

Figure 1: Initial condition for $\delta\approx 0.032$ and energy spectra for different $\delta$-values.

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\leq|{\bm{k}}|<k+1/2}|\widehat{\bm{u}}({\bm{k},t})|^{2}$$ (4)

Here $\widehat{\cdot}$ denotes the 2D Fourier transform. The initial condition, a vortex line regularized with a Gaussian cross section, does in the limit of $\delta\rightarrow 0$ approach a vortex line, i.e. a Dirac distribution. The enstrophy 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 illustrated in figure 0(b). For viscous flow simulations, considered in Caflisch et al. (2022), $\delta$ was related via $\delta=Re^{-1/2}$, where $Re=\Gamma/\nu$ is the circulation based Reynolds number with $\Gamma=\int\omega({\bm{x}},t=0)d{\bm{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 re-shifted to $\Omega_{s}=[-\pi,\pi]\times[-\pi,\pi]$, the initial circulation is

$$\int_{\Omega_{s}}\omega_{0}(\bm{x})d\bm{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$$ (5)

$$=\int_{-\pi}^{\pi}\int_{(-\pi-\phi(x))/2\delta}^{(\pi-\phi(x))/2\delta}\frac{1}{\sqrt{\pi}}\exp{\left(-v^{2}\right)}dvdx$$

$$=\int_{-\pi}^{\pi}\frac{1}{2}\Bigl{[}{\rm{erf}}(v)\Bigr{]}_{(-\pi-\phi(x))/2\delta}^{(\pi-\phi(x))/2\delta}dx=2\pi+o(\exp(-\delta^{-2})).$$

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 ${\rm erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}\exp(-v^{2})dv$ is the error function which rapidly approaches $1$ for $z\rightarrow\infty$, i.e. for $\delta\rightarrow 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 ${\rm erf}(x)=1-e^{-x^{2}}\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}})$ (Abramowitz & Stegun, 1964). The term $[{\rm{erf}}(v)]_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta}$ can therefore approximated to 2 with error of order $o(\exp(-(\pi-\frac{1}{2})^{2}\delta^{-2}))$ using that $|\phi(x)|\leq\frac{1}{2}$.

For incompressible and inviscid flow, energy $E(t)$ and enstrophy $Z(t)$ remain constant over time while the palinstrophy does increase super exponentially (Ayala & Protas, 2014). Those quantities are defined via

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

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

$$P(t) =\frac{1}{2}\int_{\Omega}|\nabla\omega(\bm{x},t)|^{2}d\bm{x}\,.$$ (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 of $\approx 1.52$ for $\delta\rightarrow 0$ (figure 2), 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, for the shifted domain $\Omega_{s}$ we get for the initial enstrophy:

$$\frac{1}{2}\int_{\Omega_{s}}|\omega_{0}(\bm{x})|^{2}d\bm{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$$ (9)

$$=\frac{1}{2}\int_{-\pi}^{\pi}\int_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta}\frac{1}{2\pi\delta}\exp{\left(-v^{2}\right)}dvdx$$

$$=\frac{1}{2}\int_{-\pi}^{\pi}\frac{1}{2\sqrt{\pi}\delta}\Bigl{[}{\rm{erf}}(v)\Bigr{]}_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta}dx=\frac{\sqrt{\pi}}{2\delta}+o(\exp(-\delta^{-2})).$$

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

$$\nabla\omega_{0}=-\left(\begin{matrix}-\phi^{\prime}(x)/\delta\\ 1/\delta\end{matrix}\right)\frac{y-\phi(x)}{\delta}\omega_{0}(x,y)$$ (10)

We thus obtain

$$\int_{\Omega_{S}}|\nabla\omega_{0}|^{2}d\bm{x} =\frac{1}{2\pi\delta^{3}}\int_{-\pi}^{\pi}\int_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta}\left(1+(\phi^{\prime}(x))^{2}\right)v^{2}\exp(-v^{2})dvdx$$ (11)

$$=\frac{1}{8\pi\delta^{3}}\int_{-\pi}^{\pi}\left(1+\frac{\cos^{2}(x)}{4}\right)\Bigl{[}\sqrt{\pi}{\rm{erf}}(v)-2v\exp(-v^{2})\Bigr{]}_{(-\pi-\phi(x))/\delta}^{(\pi-\phi(x))/\delta}dx$$

$$=\frac{9\sqrt{\pi}}{32\delta^{3}}+o(\exp(-\delta^{-2})).$$

We bound the size of this error using $\max\phi(x)$, by $\rm{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 applies for the initial circulation, matching the value $2\pi$.

Figure 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\rightarrow 0$.

3 Characteristic mapping for 2D Euler

(a)

(b)

(c)

Figure 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.

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. (2021). The numerical results in this paper are computed using a GPU Cuda implementation of the method. The open source code is available on GitHub (Yadav et al., 2023).

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 back-to-label operator for the advection of Lagrangian particles. Indeed, for any particle trajectory $\gamma(t)$ given by,

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

with initial condition $\gamma(0)=\gamma_{0}$, we have that $\chi_{B}$ satisfies,

$$\gamma_{0}=\chi_{B}(\gamma(t),t).$$ (13)

It follows by the method of characteristics, that any scalar quantity $\phi$ advected by the same velocity $\bm{u}$ satisfies $\phi(\bm{x},t)=\phi_{0}(\chi_{B}(\bm{x},t))$, this is known as the relabelling symmetry.

In the vorticity equations, $\omega$ is a scalar-valued quantity advected by $\bm{u}$ and thus has the above relabelling symmetry. The velocity field $\bm{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}+({\bm{u}}\cdot\nabla)\chi_{B}=0$$ (14a)

$$\omega({\bm{x}},t)=\omega_{0}(\chi_{B}({\bm{x}},t))$$ (14b)

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

Numerically, the advection equation for the map (14a) is discretized through the Gradient-Augmented Level-Set method (Nave et al., 2010; Seibold et al., 2012; Chidyagwai et al., 2011) 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 is eliminated as 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}<\ldots<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\dots\circ\chi_{[T,t_{m-1}]}$$ (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 $\bm{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. $\bm{x}(\theta,t=0)=(\theta,\phi(\theta)+\pi),\quad\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}\geq 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.

4 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 table 1. Two grid sizes were used: A coarser one for the description of the flow map $\chi$ and velocity $\bm{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 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}|\leq\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 at the vortex center line. Further in-depth analysis and explanation of all used parameters are reported in Yin et al. (2021); Bergmann (2022).

(a) $t=5$ and $\delta\approx 0.032$

(b) Close-up with overlayed center line

(c) Close-up of local palinstrophy

Figure 4: Developed vortices in global field and close-up of upper vortex with overlayed center-line (a,b). After an initial phase two symmetric vortices start to form with center with condensed region of high vorticity and spiral arms. Figure (c) shows corresponding local palinstrophy $|\nabla\omega(\bm{x},t)|^{2}$.

(a) Results obtained with CMM

(b) Results from Caflisch et al. (2022)

Figure 5: Comparison between our results and Caflisch et al. (2022) with $\delta\approx 0.014$ and at $t=2.85$. The CMM results were adapted to the shifted region $\Omega_{s}$ for comparison purposes. They show excellent agreement. Figure adapted from Caflisch et al. (2022), printed with their permission.

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 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 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. (2022) by their overall dynamics, as they did for their inviscid results.

Parameter	V100	A100	Parameter	V100	A100
$N_{\text{coarse}}$ & $N_{\psi}$	$8192$	$12288$	$N_{\text{fine}}$ & $N_{\omega}$	$12288$	$24576$
$\Delta t_{\text{fluid}}$	$1/24576$	$1/36864$	$\delta_{inc,b}$	$10^{-3}$	$10^{-3}$
$\epsilon_{m}$	$10^{-4}$	$5\cdot 10^{-5}$	$k_{LP}$	$4096$	$12288$
Fluid time scheme	RK4Mod	RK3Mod	Map update stencil	4th order	4th order
$N_{P}$	$10^{6}$	$10^{7}$	Particle time scheme	RK4Mod	RK3Mod

Table 1: Simulation parameters for the different runs on Nvidia V100 and A100 machines. Given are the resolution of coarse and fine grids $N_{coarse}$ and $N_{fine}$, time-step $\Delta t_{\text{fluid}}$, incompressibility threshold $\delta_{inc,b}$, local stencil size $\epsilon_{m}$ for the GALS method with its stencil order, low-pass filter size $k_{LP}$, time scheme order for fluid and particle advection and number of particles.

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 2. The settings ensure that energy and enstrophy are conserved to an error of $\approx 2\cdot 10^{-5}$, see figure 6, which is supported by the error of the CM-method being not dissipative in nature for transported quantities.

$\delta\approx$	0.0045	0.007	0.01	0.0141	0.022	0.032	0.045
$1/\delta^{2}$	$5\times 10^{4}$	$2\times 10^{4}$	$1\times 10^{4}$	$5\times 10^{3}$	$2\times 10^{3}$	$1\times 10^{3}$	$5\times 10^{2}$
Run on	A100	A100	V100	V100	V100	V100	V100
Final time	1.24	2.29	2.86	3.46	4.37	5.82	7.59
Captured windings	0	3	5	6	6	7	7

Table 2: Overview on the used vorticity layer thickness values $\delta$ including the obtained final time and the number of windings, for which computations were performed.

Simulations with successively decreasing thickness values $\delta$ reveal the flow behavior for the limit $\delta\rightarrow 0$. An example for the flow dynamics at $t=5$ and with $\delta\approx 0.032$ is depicted 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. Eventually, both vortices start to attract each other and will eventually merge into a single stable vortex for long times. With decreasing sheet thickness the roll-up process is triggered earlier and the emerging structures are smaller in scale. Figure 3(b) shows that the vortex core is round with distortion towards the second vortex due to the mutual attraction. The local palinstrophy captured in figure 3(c), is very pronounced in the spiral arms and shows the spiraling structures also present within the vortex core itself.
The acquired results are comparable with inviscid flow fields reported by Caflisch et al. (2017) and Caflisch et al. (2022). The depicted vortex in figure 4(a) is in excellent agreement with figure 14d of Caflisch et al. (2022), given in figure 4(b). 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.

(a) Relative energy error $\frac{E(t)}{E(t=0)}-1$

(b) Relative enstrophy error $\frac{Z(t)}{Z(t=0)}-1$

Figure 6: 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}$.

(a) Scaling of palinstrophy growth at early times shown in log-log.

(b) Palinstrophy growth at later time shown in log-lin.

Figure 7: 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.

5 Analysis of the results

In the following, the numerical results will be 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\rightarrow 0$.
The palinstrophy $P$ over time, shown in figure 7, 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 a strong 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)$, a relation of $s\approx\delta^{-0.77}$ was found. This means that in the limit of $\delta\rightarrow 0$ palinstrophy diveres. This is consistent with results in the literature for Navier–Stokes where exponential growth of palinstrophy was derived, see e.g. Lesieur (2008). For small times a powerlaw behavior was predicted in Ayala & Protas (2014) and estimates of the maximum palinstrophy growth were given. The lower limit of the investigated $\delta$-values, i.e. for $\delta\approx 0.0045$, experiences artifacts in form of numerically accelerated secondary Kelvin–Helmholtz instabilities along the material line, resulting in the premature and irregular steep palinstrophy growth.

With the material line $\bm{x}(\theta,t)$, given by the individual particle positions, we define following (Caflisch et al., 2022) 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}|\bm{x}_{\theta}|d\theta\,,$$ (16)

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

$$\gamma_{c}(\theta) =|\bm{x}_{\theta}|^{-1}\,.$$ (18)

with the notation $\bm{x}_{\theta}=\frac{\partial\bm{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 4 can be found in the evolution of the arc length of the material line in figure 7(a). The overall arc length is both increased by the emerging two vortices as well as 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 correspond to extreme compression of the particles. Our results for the Euler case are much more pronounced but comparable with those shown in Caflisch et al. (2022) figure 12b for Navier–Stokes.

(a) $s(\theta)$ for $\delta\approx 0.032$

(b) $\kappa_{c}$ for $\delta\approx 0.014$

(c) $\gamma_{c}$ for $\delta\approx 0.014$

Figure 8: 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}$ (b) over arc length $s(\theta)$ around vortex center for $\delta=0.014$.

The curvature shown in figure 7(b) showcases the strength of the winding process, forming two peaks of opposite sign around the vortex cores. These are situated at the edge of the coalesced 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\cdot 10^{-4}\times N_{p}$ 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 $\gamma_{c}>1$). Once the vortices start to rotate, two spiral arms will form and elongation occurs ($\gamma_{c}<1$) outside the vortex core. Both quantities reach maximum peak values at different times, depicted by the green curve in figure  7(c), and then start to disperse.

(a) Curvature $\kappa_{c}(\theta)$

(b) True vortex strength $\gamma_{c}(\theta)$

Figure 9: Temporal re-scaling for the maximum values of the curvature $\kappa_{s}$ and true vortex strength $\gamma_{s}$ as well as temporal re-scaling for the vortex turnover with empirically determined coefficients. The critical times $t_{s}$, given in figure (b), were found to be decreasing with lower $\delta$-values.

(a) Vortex turnover

(b) Rescaled palinstrophy growth in log-log.

Figure 10: 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).

Tracking the maximum value of curvature and vortex strength over time enables comparisons between different $\delta$-values, shown in figure 9. Here, the curves have been matched empirically for the given scaling laws. The curvature was found to empirically scale with $\delta^{-0.9}$ and the true vortex strength with $\delta^{-0.31}$, for which an explanation was not found. Nonetheless, this scaling supports the conjecture of Caflisch et al. (2022) that only for the limit of $\delta\rightarrow 0$ the curvature as well as true vortex strength go towards infinity for finite time. For this limit the vortex core coalesces to an individual point with infinite curvature.
Interestingly, over time both quantities show the same behavior with an initial steep increase 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. Two regions are formed based on the true vortex strength: One close to the center of the vortex were 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 some point the vortex core is reformed into a circular shape and rotates with quasi-constant speed due to the resulting velocity from the Biot–Savart law.

(a) Re-scaled by $\delta^{-0.9}$

(b) Re-scaled by $\delta^{-1}$, as done in (Caflisch et al., 2022)

Figure 11: Spatial re-scaling with different factors for the fourth time that $\partial\bm{x}(\theta,t)/\partial\theta\approx 0$ in the $y$ direction for different $\delta$-values.

This is also visible in figure 9(a). Here the occurrence of the vanishing gradient is plotted over re-scaled 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{\bm{x}}(\theta,t)}{\partial\theta}$ vanishes. As the vortex turnover over re-scaled time occurs in constant periods (linear growth in figure 9(a)), 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 9(b). In all graphs of figures 9 and 10, 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. (2022) 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. (2022) results from different flow dynamics. With no viscosity 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$.

(a) $\delta=0.045$

(b) $\delta=0.01$

Figure 12: 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 0(b).

The material curves for different $\delta$-values have strong spatial self-similarity. However, while Caflisch et al. (2022) reported a scaling of $\delta^{-1}$, this did not match our computed results. These show an empirical scaling of $\delta^{-0.9}$, which is slightly different, as shown in figure 11. The 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$).
The characteristic of the initial energy spectra, defined in equation 4, have already been explained in section 2. For larger $\delta$-values the initial profile with exponential decay experiences a shift of energy to finer scales (figure 11(a)). 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 11(b)). However, for both simulations the two vortices have already emerged for the captured time and performed several windings, while the results for figure 11(b) did not yet observe any developed distortion from the vortex merging process. 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 merging process with increasing influence on the flow field, the shape of the vortex core becomes more and more elliptic (figure 13). This brings an imbalance to the 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.

(a) $t=5$

(b) $t=6.2$

(c) $t=7.4$

Figure 13: 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.

6 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 (Pauls & Frisch, 2007).

6.1 Borel–Polya–van der Hoeven method

The Borel–Polya–van der Hoeven (BPH) method (Pauls & Frisch, 2007) 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 (Pólya, 1933) with numerical techniques for asymptotic interpolation (van der Hoeven, 2006). This method has been successfully applied to find singularities for the 1D Burgers equations (Pauls & Frisch, 2007) and has been used to analyze potential singularity formation in the 2D Euler vortex sheet problem (Caflisch et al., 2022). The following is a summary of the method along with some small modifications implemented to adapt to our case.

Given a complex function $f(Z)$ with formal power series

$$f(Z)=\sum_{n=0}^{\infty}a_{n}Z^{n},$$ (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},\quad\quad f^{BL}(Z)=\frac{1}{Z}f\left(\frac{1}{Z}\right)=\sum_{n=0}^{\infty}\frac{a_{n}}{Z^{n+1}},$$ (20)

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},$$ (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}.$$ (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)),$$ (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.

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}}+\dots+\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,\ldots,m}|c_{j}|\cos(\phi_{j}-\theta).$$ (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}$$ (25)

which is decomposed as the sum of two functions

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

$$u^{-}(z)=\sum_{k>0}\hat{u}_{-k}e^{-ikz}.$$ (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^{-z}$ and $Z=e^{-iz}$ are made to $u^{+}$ and $u^{-}$ respectively to yield

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

$$f_{2}(Z)=u^{-}(i\ln Z),$$ (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.

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},\quad\quad\theta_{k}=\frac{2\pi k}{K}$$ (28)

for $j=0,1,2,\ldots,M-1$ and $k=0,1,2,\ldots,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 kn}{N}}=\sum_{n=0}^{N-1}a_{n}\frac{r_{j}^{n}}{n!}e^{-i\theta_{k}n},$$ (29)

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

$$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).$$ (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},\ldots,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},\ldots,Z^{*}_{s}\}$ for $f$, where $Z^{*}_{j}=1/y_{j}^{*}$, which are now the singularities of $f$ closest to 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}$.

(a)

(b)

Figure 14: 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.

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

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

$$\kappa=\frac{|\dot{\gamma}\times\ddot{\gamma}|}{|\dot{\gamma}|^{3}}.$$ (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 14, 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 the results found for the viscous case in (Caflisch et al., 2022). 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.

7 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. (2022) 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 have a region of very strong compression around the vortex centers and of elongation in the spiral arms that form. These form a peak in compression directly at the center and scale with decreasing vortex sheet thickness $\delta$. For the limit of $\delta\rightarrow 0$ the vortex strength forms 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 edge of the condensed vortex blobs. These are of opposite sign and describe the transition from the vortex blob formation to the spiral arms and additionally scale with 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. (2022). These quantities also show strong self-similarity in time and space. With the help of the curvature, vortex strength and turnover time, the temporal rescaling 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 rescaling was not found to be consistent with that of the BR-equation and similarities are clearly present. 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. (2022) in the viscous case, i.e. for Navier–Stokes. 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 vorticty 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 et al. (2023).

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 Kerr (1993); Hou & Li (2007) initial condition can be found in Yin et al. (2023). High resolution 3D computations considering different flow configurations, like Kerr (1993); Hou & Li (2007) and more recently the one by Moffatt & Kimura (2019a, b), will be published in forthcoming work.