Tobias Hertel
email: tobias.hertel@oca.eu

Nicolas Besse
email: nicolas.besse@oca.eu

Uriel Frisch

The Cauchy-Lagrange method for 3D-axisymmetric wall-bounded and potentially singular incompressible Euler flows

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.

Introduction

The setting is that of the basic three-dimensional incompressible Euler equations in a wall-bounded infinite cylinder D, with boundary ∂ D,

v t + v • ∇v = -∇p, ∇ • v = 0, in D, v • ν = 0, on ∂ D, v | t=0 = v 0 , ( 1
)
where ν is the unit outward normal vector to ∂ D and the velocity field v = v(a, t) is assumed to be axisymmetric and periodic along the cylinder. In the classical sense, the Euler equations are solvable, at least until a finite time of existence T > 0, whenever v 0 ∈ C1,α , α > 0. In this case, one defines the Lagrangian particle trajectories via the characteristic map X, unique solution to the characteristic equation ∂ t X (a, t) = v (X (a, t) , t) , X (a, 0) = a,

for a ∈ D, t ∈ [0, T ).

A historical remark is in order here. The above equations still constitute some of the hardest problems in physics and applied mathematics. Especially the Lagrangian view point, using the particle trajectories, has proved to be a useful tool to address those problems. After the postulation by Euler [START_REF] Euler | Principes généraux du mouvement des fluides[END_REF] in 1757, the hydro-dynamical part of the equations (without the rigid boundary conditions) were well known to Cauchy in 1815. He proposed rewriting the equations in Lagrangian coordinates, thereby introducing the now called Cauchy-invariants formulation (CIF) [START_REF] Cauchy | Sur l'état du fluide à une époque quelconque du mouvement[END_REF][START_REF] Frisch | Cauchy's almost forgotten Lagrangian formulation of the Euler equation for 3D incompressible flow[END_REF][START_REF] Besse | Geometric formulation of the Cauchy invariants for incompressible Euler flow in flat and curved spaces[END_REF]. In 1858, Helmholtz stated remarkable results concerning vorticity conservation laws. Even though these results were expressed in Lagrangian language, the proof was established using Eulerian variables. Shortly after this, Riemann noticed that these results could perhaps be derived directly using Lagrangian variables. Due to Riemann's suggestion, Goettingen university proposed a prize for giving a direct Lagrangian proof of Helmholtz's theorems. This prize was awarded to Hankel [START_REF] Villone | Hermann Hankel's "On the general theory of motion of fluids[END_REF], who proved Helmholtz's vorticity conservation laws in a purely Lagrangian framework using the Cauchy-invariants formulation, see also [START_REF] Frisch | A contemporary look at Hermann Hankel's 1861 pioneering work on Lagrangian fluid dynamics[END_REF]. However, it was only in the beginning of the nineteenth century, that a rigorous mathematical proof for the above existence and uniqueness result was proved by Gyunter [START_REF] Gyunter | Über ein Hauptproblem der Hydrodynamik[END_REF] and Lichtenstein [START_REF] Lichtenstein | Über einige Existenzprobleme der Hydrodynamic[END_REF] for the case of the whole space, and by Ebin and Marsden [START_REF] Ebin | Groups of diffeomorphisms and the motion of an incompressible fluid[END_REF] in 1970, for the case of bounded domains and manifolds with boundaries (see also [START_REF] Besse | Regularity of the Geodesic Flow of the Incompressible Euler Equations on a Manifold[END_REF] for the case of manifolds with boundaries and in the context of the Cauchy-Lagrange approach). At the present moment, one of the most important questions, namely the global existence and uniqueness of an initially smooth solution, is still open. Current mathematical and numerical work (e.g. [START_REF] Elgindi | Finite-time singularity formation for strong solutions to the axi-symmetric 3D Euler equations[END_REF][START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF][START_REF] Chen | Finite time blowup of 2D Boussinesq and 3D Euler equations with C 1,α velocity and boundary[END_REF]) suggests that the presence of boundaries plays a crucial role in answering this question.

Remarkably, the particle trajectories, also called characteristics, are analytic in time, i.e. the characteristic map X is the sum of an absolutely converging Taylor series locally in time [START_REF] Serfati | Étude mathématique de flammes infiniment minces en combustion[END_REF][START_REF] Serfati | Structures holomorphes à faible régularité spatiale en mécanique des fluides[END_REF][START_REF] Gamblin | Système d'Euler incompressible et régularité microlocale analytique[END_REF][START_REF] Shnirelman | On the analyticity of particle trajectories in the ideal incompressible fluid[END_REF][START_REF] Glass | Smoothness of the motion of a rigid body immersed in an incompressible perfect fluid[END_REF][START_REF] Frisch | A very smooth ride in a rough sea[END_REF][START_REF] Constantin | Analyticity of Lagrangian trajectories for well posed inviscid incompressible fluid models[END_REF][START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF][START_REF] Besse | Regularity of the Geodesic Flow of the Incompressible Euler Equations on a Manifold[END_REF]. Thus, we may conveniently write X (a, t; t 0 ) = ∞ s=0 X (s) (a) (t -t 0 ) s , t 0 ∈ [0, T ),

where t must lie within the disc of convergence with radius around t 0 . In the sequel of this article, we will drop the argument t 0 in eq. ( 3), because we use t 0 = 0 throughout the text. All the referenced proofs show that the temporal analyticity of X holds even under the lowest regularity assumptions on the initial flow, namely v 0 ∈ C 1,α , α > 0 (or a weak alternative). While most of the proofs in the cited references above make use of singular integral calculus and start directly from the Euler equations, the references [START_REF] Frisch | A very smooth ride in a rough sea[END_REF][START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF][START_REF] Besse | Regularity of the Geodesic Flow of the Incompressible Euler Equations on a Manifold[END_REF]] make use of the revived CIF. This reformulation of the incompressible Euler equations 1 is written entirely in terms of the characteristic map and gives rise to a constructive proof of the temporal analyticity of X. As a byproduct, the authors obtain recursion relations for the time-Taylor coefficients X (s) of eq. [START_REF] Beale | Remarks on the breakdown of smooth solutions for the 3-D Euler equations[END_REF]. In [START_REF] Podvigina | The Cauchy-Lagrangian method for numerical analysis of Euler flow[END_REF], it has been shown that this constructive nature of the CIF-approach leads to an effective algorithm for an incompressible Euler-flow in two dimensions with periodic boundary conditions. This numerical scheme, known as the Cauchy-Lagrange algorithm (CLA), belongs to the large class of semi-Lagrangian methods which are very popular in the fields of numerical weather prediction and numerical plasma physics, and also fluid mechanics and gas dynamics. The present work is inspired by the preceding papers and builds on the findings in [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF]. In this work, we describe the CLA for wall-bounded Euler flow and investigate numerically a periodic, and axisymmetric incompressible flow in the cylinder2 D.

A common numerical method to advance an incompressible Euler flow in time is the Runge-Kutta time-stepping scheme on an underlying Eulerian mesh. A major drawback of this method is the limitation of the time-step size depending on the mesh spacing and the flow itself, well known as the Courant-Friedrichs-Lewy (CFL) condition. This is problematic especially in turbulent flows, where a necessarily high resolution enforces narrow time-steps and, thus, slows down the overall algorithm. The time-step in the CLA is only limited by the convergence radius , which behaves roughly as cst./ ω 0 α where ω 0 is the vorticity field at t 0 and . α denotes the Hölder norm. Therefore, the algorithm is entirely independent of a CFL condition. This advantage makes the CLA a suitable candidate for simulations of potentially singular flows, where high resolutions are needed to accurately represent large developing gradients. A celebrated result connects the blow-up event to the time-integrated vorticity modulus [START_REF] Beale | Remarks on the breakdown of smooth solutions for the 3-D Euler equations[END_REF][START_REF] Ferrari | On the blow-up of solutions of the 3-D Euler equations in a bounded domain[END_REF]. More accurately, if a blow-up occurs at a finite time T * , then the timeintegrated L ∞ -norm of the vorticity ω until T * becomes infinite and in particular lim sup

t→T * ω L ∞ = ∞. (4)
In the spirit of this result, broad numerical studies have been executed in the last decades, e.g. [START_REF] Grauer | Numerical computation of 3D incompressible ideal fluids with swirl[END_REF][START_REF] Pumir | Development of singular solutions to the axisymmetric Euler equations[END_REF][START_REF] Kerr | Evidence for a singularity of the three-dimensional, incompressible Euler equations[END_REF], mainly searching for an initially smooth flow that eventually verifies [START_REF] Bennett | Lagrangian fluid dynamics[END_REF]. See also [START_REF] Gibbon | The three-dimensional Euler equations: Where do we stand?[END_REF] for a thorough inspection of the concerned literature. However, even with steadily growing computational power, it seemed that the blow-up problem could hardly be answered numerically, because earlier claims of finite time singularities were later attributed to numerical artifacts or under-resolution (see e.g. [START_REF] Shu | Small-scale structures in Boussinesq convection[END_REF][START_REF] Hou | Dynamic depletion of vortex stretching and non-blowup of the 3-D incompressible Euler equations[END_REF][START_REF] Bustamante | 3D Euler about a 2D symmetry plane[END_REF]). Nonetheless, efforts have been continued and now reference [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] provides ample numerical evidence for a blow-up of an axisymmetric flow in a cylinder. To verify their blow-up result, either a rigorous mathematical proof or more simulations of the same singular flow by other, different, numerical schemes have to be established. We believe that the here presented CLA can take on that challenge since almost no parallels to the algorithm in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] can be drawn, as discussed in Section 10.

The absence of a CFL condition makes the CLA a well-suited candidate to approach the strong gradients near the potential blow-up with larger time-steps on arbitrarily fine3 grids. Furthermore, special attention has been employed to solve the appearing Poisson problems directly for the derivatives of the solutions. This minimizes the errors in the computation of the time-Taylor coefficients X (s) and allows to choose larger time-steps in the disc of convergence. The CLA is implemented in Fortran for high performance and is found to be accurate, stable, and fast. Furthermore, a wide parameter palette is available, including different Poisson solving schemes, interpolation methods (with refinement options), truncation and time-stepping settings, limits for conserved quantities, de-aliasing options, etc. The present work is refined from [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF], where the implemented CLA and the connected theory is discussed in greater detail. The article is organized as follows. An overview of the different steps of the CLA is given in Section 2. The technical frame in form of the recursion relations for the time-Taylor coefficients X (s) and the used Helmholtz-Hodge decomposition (HHD) of X (s) in cylindrical coordinates are stated in Sections 3 and 4. Particular attention has been attributed to the appearing Poisson problems in the HHD. Their treatment can be found in Section 5 together with the used spectral Chebyshev-Fourier discretization. How the flow fields (velocity and vorticity) are re-interpolated onto the original Eulerian grid, after they have been computed along the particle trajectories, is described in Section 6. In Section 7, the choice of a possible time-step is discussed, and validation criteria for one iteration are presented in Section 8. Subsequently, the CLA is tested on stationary and swirl-free (non-stationary) solutions in Section 9, before it is applied to the potentially singular solution of [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] in Section 10.

Description of the Cauchy-Lagrange Algorithm

The main goal of the Cauchy-Lagrange algorithm (CLA) is the computation of the particle trajectory map X as defined in eq. ( 3). The Cauchy-invariants formulation of the three-dimensional incompressible Euler equations (1), together with the incompressibility condition, reads k=1,2,3 ∇ Ẋk (a, t) × ∇X k (a, t) = ω 0 (a), [START_REF] Besse | Adaptive multiresolution semi-Lagrangian discontinuous Galerkin methods for the Vlasov equations[END_REF] det ( DX(a, t)) = 1,

in the Cartesian coordinates a = (a 1 , a 2 , a 3 ) . The symbol ω 0 denotes the initial vorticity and D indicates the total derivative with respect to the spatial variable a, ∇ denotes the nabla operator and the dot-notation has been adopted to indicate temporal derivatives. Equations ( 5) and [START_REF] Besse | Regularity of the Geodesic Flow of the Incompressible Euler Equations on a Manifold[END_REF] give rise to recursion relations for the curl and divergence of the time-Taylor coefficients X (s) , such that the latter can be expressed in terms of lower order time-Taylor coefficients X (s ) , s < s. From [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF], we gather those recursion relations, but abbreviate here the appearing sums by calligraphic capital letters F and G to simplify the notation. Note that, in the given reference, the recursions are given in terms of the displacement map ξ = X -Id, but since the Taylor coefficients are time derivatives at t = 0, we have ξ (s) = X (s) , for s > 0 and X (0) = Id, the identity map. The first recursion step starts with X (1) = v 0 , and therefore, ∇ • X (1) = 0 and -∇ × X (1) = -ω 0 . For a time-Taylor coefficient with orders s > 1 we have ∇ • X (s) = F ( DX (s ) ) 0<s <s [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF] -∇ × X (s) = G ( DX (s ) ) 0<s <s ,

Subsequently, the time-Taylor coefficient X (s) is expressed by a Helmholtz-Hodge decomposition X (s) = ∇φ+∇×Ψ , where the potentials result from two Poisson problems, namely ∆ φ = ∇•X (s) and ∆ Ψ = -∇ × X (s) . Those Poisson problems are supplemented with appropriate boundary conditions that are obtained from the boundary defining map S : R 3 → R, for which S(X(∂ Ω, t)) = 0 holds true for all times t ≥ 0. In the case of a periodic cylinder, the explicit boundary conditions are specified in Section 4. After a predefined number of recursions S, we may proceed with the time-step insertion by directly computing

X (a, t) ≈ S s=0 X (s) (a) t s , (9)
provided a suitable time-step t has been found. The only limitation of the latter is the radius of convergence of the time Taylor series. This radius is typically of the form = C/ ω 0 C α , but the constant C = C(Ω, α) is difficult to obtain analytically. We have achieved good results in setting the time-step small enough such that the highest order terms do not contribute to the sum anymore for moderate S due to limitations in machine precision. For a detailed inspection of the time-stepping, we refer the reader to Section 7. Once a time-step has been set, we obtain the velocity field on the trajectories via the characteristic equation

u (X (a, t) , t) = X t (a, t) ≈ S s=1 sX (s) (a) t s-1 , ( 10
)
the vorticity is similarly obtained via the vorticity transport formula (i.e. [46, p. 20])

ω (X (a, t) , t) = DX(a, t) • ω 0 (a) ≈ S s=0 DX (s) (a) • ω 0 (a) t s . (11)
Numerically, the characteristics deform the underlying discrete grid and the velocity and vorticity fields are only known on that distorted grid, i.e. at the end points of the particle trajectories. Hence, an interpolation is needed to obtain the values of the advanced flow fields on the Eulerian grid. This three-dimensional scattered interpolation needs, in principal, only to be performed for the vorticity field, but it is better to re-interpolate the velocity field as well, because it appears in the 1st order term in the next time-Taylor series of the characteristics (in X (1) ). These interpolated flow fields are then taken as the new initial conditions and, therefore, we may set t 0 = 0 after each time-step, where t 0 appeared in eq. ( 3). In this paper, a three-dimensional axisymmetric flow in a cylindrical domain is modeled to understand the challenges and advantages of this novel algorithm for wall-bounded incompressible Euler flow. The axisymmetry of the problem simplifies both the Poisson solvers and the interpolation procedure to two dimensions. Best results in the interpolation step have been achieved using the cascade interpolation scheme [START_REF] Purser | An efficient interpolation procedure for high-order three-dimensional semi-Lagrangian models[END_REF].

Recursion relations for axisymmetric flow

In this section, we would like to introduce our notations for cylindrical domains and give the recursion relations for our particular setting of an axisymmetric and periodic flow in the cylindrical region

D(1, L) := a = (a 1 , a 2 , a 3 ) ∈ R | a 2 1 + a 2 2 ≤ 1 and a 3 ∈ [0, L] . (12)
Note that the real flow takes place inside of an infinite cylinder, but because of the periodicity we can restrict our computations to the region D(1, L). The Lagrangian particles, however, may leave this region and disperse into the whole cylinder, but they can always be remapped into D(1, L). See Section 6 for a more detailed discussion on that matter. An axisymmetric flow in

D(1, L) takes the form v(a, t) = v r (r, z, t)e r + v α (r, z, t)e α + v z (r, z, t)e z , (13)
where v and a = re r + ze z are represented in the cylindrical basis e r = (cos α, sin α, 0) , e α = (-sin α, cos α, 0) , e z = (0, 0, 1) . As usual we have, for a = (a 1 , a 2 , a 3 ) ,

r = a 2 1 + a 2 2 , α = A(a 1 , a 2 ), z = a 3 , with (r, α, z) ∈ [0, 1] × (-π, π] × [0, L] and A(a 1 , a 2 ) :=            0, if a 1 = a 2 = 0 arcsin( a2 r ), if a 1 = 0, a 2 = 0 arctan( a2 a1 ), if a 1 > 0 π -arcsin( a2 r ), if a 1 < 0 . (14)
The cylindrical vorticity corresponding to eq. ( 13) reads

ω(a, t) = ω r (r, z, t)e r + ω α (r, z, t)e α + ω z (r, z, t)e z , (15)
and is obtained by

ω = ∇ c × v = -v α z e r + (v r z -v z r )e α + 1 r (rv α ) r , e z , (16)
where

∇ c := e r ∂ r + 1 r e α ∂ α + e z ∂ z (17)
denotes the cylindrical nabla operator. The components of v and ω are independent of the angular variable α and periodic in z ∈ [0, L]. The α-dependence of the flow fields only lies in the basis vectors e r and e α . The same is true for the the characteristics X = X(a, t) and the displacement map ξ := X -Id, which are denoted by

X(a, t) = X r (r, z, t)e r + X α (r, z, t)e α + X z (r, z, t)e z , (18)
ξ(a, t) = R(r, z, t)e r + A(r, z, t)e α + Z(r, z, t)e z . (19)
We recall the Cauchy-invariants formula in Cartesian coordinates in terms of the displacement map ξ = (ξ 1 , ξ 2 , ξ 3 ), namely eqs. ( 16)-( 17) in [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF],

∇ × ξ + 1≤k≤3 ∇ ξk × ∇ξ k = ω 0 (20) ∇ • ξ + 1≤i<j≤3 (∂ i ξ i ∂ j ξ j -∂ i ξ j ∂ j ξ i ) + det ( Dξ) = 0. (21)
Here, the dot-notation is used to indicate temporal derivatives. An insertion of eq. ( 19) into eqs. [START_REF] Euler | Principes généraux du mouvement des fluides[END_REF] and ( 21) yields after a slight rearrangement

-∇ c × ξ = e r 1 r ∂ z ṘA -ȦR + e α Ṙ, R + Ȧ, A + Ż, Z + e z 1 r ∂ r ȦR -ṘA -ω 0 (22) ∇ c • ξ = 1 r A A, Z -A r + R R, Z -R r -Z z + R, Z , (23)
where f, g := f z g r -f r g z denotes the Poisson bracket for functions f = f (r, z) and g = g(r, z). An insertion of the formal time-Taylor series

Λ(r, z, t) = s≥1 Λ (s) (r, z) t s , Λ(r, z, t) = s≥1 sΛ (s) (r, z) t s-1 , Λ = R, A, Z, (24)
into the expressions of the rotation (eq. ( 22)) and divergence (eq. ( 23)) above yields component-wise

(-∇ c × ξ (s) ) r = -δ 1s ω r 0 + 1≤m<s m rs ∂ z R (m) A (s-m) -A (m) R (s-m) (25)
(-∇ c × ξ (s) ) α = -δ 1s ω α 0 + 1≤m<s m s R (m) , R (s-m) + A (m) , A (s-m) + Z (m) , Z (s-m) (26)
(-∇ c × ξ (s) ) z = -δ 1s ω z 0 + 1≤m<s m rs ∂ r A (m) R (s-m) -R (m) A (s-m) , (27)
and

∇ c • ξ (s) = R, Z (s) + 1≤m<s A (m) r A, Z (s-m) -A (s-m) r + R (m) r R, Z (s-m) -R (s-m) r -Z (s-m) z (28) with Λ, Z (n) := 1≤m<n Λ (m) z Z (n-m) r -Λ (m) r Z (s-m) z , n = 1, 2, . . . , such that Λ, Z = 1≤s Λ, Z (s)
t s for Λ = R, A. By convention, empty sums evaluate to zero.

Numerically it makes sense to store the time-Taylor coefficients of R, Z and A, Z . Although this increases storage demand, it speeds up the recurrence computations tremendously for high orders, as the summation terms in eq. ( 28) consist only of two factors instead of three.

In sight of the envisaged Helmholtz-Hodge decomposition in Section 4 we state the boundary value R(1, z, t) for t > 0. As a matter of fact, for a boundary point a, X(a, t) lies on the boundary for all t > 0. In the case of our periodic flow in the cylinder D(1, L), this boundary reduces to the mantle of the cylinder which is discribed by the kernel of the boundary map S(a 1 , a 2 , a 3 ) = a 2 1 + a 2 2 -1. An insertion of X(a, t), a ∈ ∂D(1, L), yields, in terms of the cylindrical components of ξ,

0 = S(X(a, t)) = X 2 1 (a, t) + X 2 2 (a, t) -1 = (R(1, z, t) + 1) 2 + A(1, z, t) 2 -1 = 2R(1, z, t) + R(1, z, t) 2 + A(1, z, t) 2 . ( 29
)
and, thus,

R(1, z, t) = - 1 2 R(1, z, t) 2 + A(1, z, t) 2 .
An insertion of the power series of the components (24) into the last term, with subsequent grouping of coefficient orders s ≥ 1, entails

R (s) (1, z) = - 1 2 s-1 n=1 R (n) R (s-n) + A (n) A (s-n) (1, z). ( 30
)
We have stated all the needed recursion formulas for the rotation, divergence and boundary terms of the coefficients ξ (s) (they equate the X (s) ), which are needed for the Helmholtz-Hodge decomposition that is discussed in the next section. Note that those recursion formulas depend only on (r, z) ∈ [0, 1]×[0, L] due to the imposed axisymmetry of the flow, which reduces the computational complexity tremendously when compared to general 3D flow.

Helmholtz-Hodge Decomposition and Computational Simplifications

In order to exploit the formerly found recursion relations of the divergence and rotation we state the Helmholtz-Hodge decomposition (HHD) of the trajectory's time-Taylor coefficients in cylindrical coordinates. In the following, let s ≥ 1 be fixed. ξ (s) can be expressed by a HHD in the cylindrical domain D(1, L) as

ξ (s) = ∇φ + ∇ × Ψ . (31)
The potentials φ and Ψ should read φ (s) and Ψ (s) , since those potentials depend on s. We omit the superscripts only to increase the readability. The scalar potential φ is uniquely described by the Poisson problem

∆ φ = ∇ • ξ (s) , in D(1, L), ∂ ν φ = ν • ξ (s) , on ∂ D(1, L), (32)
and the vector potential

Ψ by ∆ Ψ = -∇ × ξ (s) , ∇ • Ψ = 0, in D(1, L), Ψ × ν = 0, on ∂ D(1, L), (33)
where ν denotes the unit outward normal vector to the boundary, and in our case we simply have ν = e r . See [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF] for a detailed analysis of the HHD and the above statements for a bounded domain.

In principle, we treat a flow inside of an infinite cylinder, which is unbounded in vertical direction, but we do not need to impose a decay of the solution for z → ∞, because of the periodicity in [0, L].

In this particular case, the validity of the above statement is verified without difficulty by a slight adaptation of the arguments in [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF]. Note that the boundary conditions for Ψ are different from those in [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF]. Indeed, Ψ = 0 on the boundary of the domain is not sufficient for a general HHD in a bounded domain, but accounts only for a special situation. We refer the reader to [START_REF] Besse | Regularity of the Geodesic Flow of the Incompressible Euler Equations on a Manifold[END_REF], where the boundary value problem eq. ( 33) from above has been addressed in Remark 6. We reformulate the above problems in cylindrical coordinates a = (a 1 , a 2 , a 3 ) → (r cos α, r sin α, z), and keep the symbols for the vector potentials by an abuse of notation. Thus, we are looking for functions φ = φ(r, z) and

Ψ = Ψ (r, α, z), with Ψ (r, α, z) = Ψ r (r, z)e r + Ψ α (r, z)e α + Ψ z (r, z)e z , (34)
defined in the domain D(1, L). Under the assumption of axisymmetry we may fix momentarily α = 0 which implies e r = (1, 0, 0) . The boundary conditions for φ translate to

∂ r φ | r=1 = e r • ξ (s) | r=1 = R (s) (1, z),
and those for Ψ , namely Ψ | r=1 ×e r = 0, imply

Ψ z (1, z) = Ψ α (1, z) = 0 . (35)
Moreover, the divergence-free condition,

∇ • Ψ = 1 r Ψ r + ∂ r Ψ r + ∂ z Ψ z = 0
, is assumed to hold up to the boundary and delivers the boundary condition on Ψ r ,

Ψ r (1, z) + ∂ r Ψ r (1, z) = -∂ z Ψ z (1, z) = 0, (36)
where eq. ( 35) is used. Let us define the linear differential operators

L := 1 r ∂ r (r . ) = 1 r + ∂ r , (37)
L 1 := L(∂ r . ) = 1 r ∂ r (r∂ r . ) = ∂ 2 r + 1 r ∂ r , (38)
L 2 := ∂ r (L . ) = ∂ r 1 r ∂ r (r . ) = ∂ 2 r + 1 r ∂ r - 1 r 2 , (39)
G := L 1 + ∂ 2 z = ∂ 2 r + 1 r ∂ r + ∂ 2 z , (40)
H := L 2 + ∂ 2 z = ∂ 2 r + 1 r ∂ r - 1 r 2 + ∂ 2 z . (41)
Subsequently, the translation of eqs. [START_REF] Grauer | Numerical computation of 3D incompressible ideal fluids with swirl[END_REF] and [START_REF] Gyunter | Über ein Hauptproblem der Hydrodynamik[END_REF] in cylindrical coordinates entails four partial differential equations

G φ(r, z) = ∇ • ξ (s) (r, z) , ∂ r φ(0, z) = 0, (42)
∂ r φ(1, z) = R (s) (1, z), H Ψ r (r, z) = (-∇ × ξ (s) ) r (r, z), Ψ r (0, z) = 0, (43)
Ψ r (1, z) + ∂ r Ψ r (1, z) = 0, H Ψ α (r, z) = (-∇ × ξ (s) ) α (r, z), Ψ α (0, z) = 0, (44)
Ψ α (1, z) = 0, G Ψ z (r, z) = (-∇ × ξ (s) ) z (r, z), ∂ r Ψ z (0, z) = 0, (45)
Ψ z (1, z) = 0, for φ, Ψ r , Ψ α , Ψ z ∈ C 2 L ([0, 1] × [0, L])
, where the subscript ( . ) L indicates periodicity in [0, L]. The radial regularity at the boundary is understood to be a limit. The pole conditions stem from the continuity constraints, as the operators must stay bounded for r ∈ [0, 1]. After solving the above equations, one obtains the coefficients A (s) , R (s) and Z (s) in cylindrical coordinates via eq. ( 31), i. e.

R (s) = ∂ r φ -∂ z Ψ α (46)
A (s) = ∂ z Ψ r -∂ r Ψ z (47)
Z (s) = ∂ z φ + ∂ r Ψ α + 1 r Ψ α , (48)
where all functions only depend on (r, z) ∈ [0, 1] × [0, L].

Computational simplifications for the implementation

In consulting the recursion relations eqs. ( 25) to [START_REF] Gautschi | Orthogonal polynomials[END_REF], the overall algorithm demands the spatial derivatives of the time-Taylor coefficients R (s) , A (s) and Z (s) in radial and vertical direction, i.e. second derivatives on the potentials. Second-order derivatives imply large numerical errors, especially if produced via standard spectral differentiation methods, which potentially limits the possible timestep that can be inserted into the truncated time-Taylor series. Therefore, we aim to directly solve for the needed second-order derivatives, which is discussed in later sections. In the following, we state the derivatives of φ and the components of Ψ that need to be found. Fortunately, not all (mixed) second-order derivatives of the potentials need to be solved for, but only a few, which results from the simplifications that are discussed here. While eqs. [START_REF] Lamb | Hydrodynamics[END_REF], [START_REF] Lichtenstein | Über einige Existenzprobleme der Hydrodynamic[END_REF] and [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] are pure Neumann or Dirichlet type boundary value problems, is eq. ( 43) a mixed boundary value problem. However, in consulting eqs. [START_REF] Majda | Vorticity and Incompressible Flow[END_REF] to [START_REF] Mercier | Analyse numérique des méthodes spectrales[END_REF], the component Ψ r is only used for the construction of A (s) and derivatives thereof. In the following we show that it is not necessary to solve eqs. [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF] and [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] in the case of our cylindrical geometry. Instead of solving eqs. [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF] and [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] for Ψ r , Ψ z , we will directly solve for A (s) and its derivatives. As a matter of fact, we have

∂ z A (s) = -(∇ × ξ (s) ) r (49)
LA (s) = A (s) r + ∂ r A (s) = (∇ × ξ (s) ) z , (50)
stemming from the definition of the curl operator in cylindrical coordinates, namely

∇ c × ξ (s) = -A (s) z e r + (R (s) z -Z (s) r )e α + 1 r (rA (s) ) r e z , (51)
with ∇ c defined in eq. ( 17). The RHSs of eqs. ( 49) and ( 50) are given via the recursion relations eqs. ( 25) and [START_REF] Gamblin | Système d'Euler incompressible et régularité microlocale analytique[END_REF], which are computed and saved in the CLA. Integrating the operator L = 1 r ∂ r (r . ) in eq. ( 50) directly yields

A (s) (r, z) = 1 r r 0 t (∇ × ξ (s) ) z (t, z) dt (52)
A (s) 1 (r, z) := A (s) (r, z) r = 1 r 2 r 0 t (∇ × ξ (s) ) z (t, z) dt, (53)
where we have used that A (s) (0, z) = 0 for all z ∈ [0, L] and s ∈ N 0 . However, our implementation of the CLA solves directly for ∂ r A (s) by calculating ∂ r L -1 in a later defined way. Subsequently, A

1 is obtained in using eq. ( 50), i.e. A (s) s) , and

1 = (∇ × ξ (s) ) z -∂ r A (
A (s) = rA (s) 1 . Furthermore, provided we know ∂ r R (s) , ∂ z R (s) and R (s) 1 , we compute in a straight forward way ∂ r Z (s) = ∂ z R (s) + (-∇ × ξ (s) ) α (54)
∂ z Z (s) = -R (s) 1 -∂ r R (s) + ∇ • ξ (s) , (55)
where we have used again eq. ( 51) and the definition of the divergence, namely

∇ c • ξ (s) = 1 r (rR (s) ) r + Z (s) z . (56)
The RHSs of eqs. ( 54) and ( 55) are given by the recursion relations eqs. ( 26) and [START_REF] Gautschi | Orthogonal polynomials[END_REF]. Therefore, we only need to compute those expressions of the Helmholtz-Hodge potentials that appear in the following RHSs,

∂ r R (s) = ∂ 2 r φ -∂ r ∂ z Ψ α (57)
∂ z R (s) = ∂ z ∂ r φ -∂ 2 z Ψ α (58)
R (s) 1 = ∂ r φ r - ∂ z Ψ α r (59)
Z (s) = ∂ z φ + ∂ r Ψ α + 1 r Ψ α . ( 60
)
In fact, the coefficients R (s) , A (s) and Z (s) are not needed for the recurrence mechanism itself, but for the actual computation of the particle trajectories. This information is mainly needed for the interpolation step after the time-step insertion. The coefficients R (s) and A (s) are simply obtained by multiplying R (s)

1 and A (s)

1 by r respectively.

Discretization in space and Poisson solvers

A crucial part of the investigated algorithm is the accurate computation of the Calderon-Zygmund operators of order zero, namely the double derivatives of an inverse Laplacian, applied to the divergence and rotation of the trajectory's time-Taylor coefficients. The appearing errors from the second derivatives on top of the errors in the solutions themselves may recursively cascade into higher orders s. This amplification of errors can lead to a lower order scheme and a very limited time-step. Thus, the better the above-mentioned operators are computed, the better the time-step can be chosen and the more accurate the overall algorithm becomes with fewer time-steps. It is desirable to preserve the possible range of a time-step as much as possible by maximally reducing the accumulation of numerical errors.

In this first implementation of the CLA for a domain with boundaries, we refrain from finite difference methods for the Calderon-Zygmund operators, because of their slow convergence and limited accuracy in the derivatives. Instead, a pseudo-spectral approach was chosen and a representation of the trajectory's coefficients by shifted Chebyshev series in the radial dimension and Fourier series in the vertical dimension is applied:

X (s) (r, z) = ∞ |k|=0 ∞ l=0 X (s) kl T * l (r) e i kz , (61)
where T * l (r) = T l (2r -1) are the shifted Chebyshev polynomials and

k := 2πk L . (62)
A truncation to M + 1, N + 1 ∈ N modes yields a spectrally discrete approximation of the latter, i.e.

X (s) (r, z) ∼ = (I N M X (s) )(r, z) = M/2 |k|=0 N l=0 X (s) kl T * l (r) e i kz . (63)
See Section A for a basic overview of the shifted Chebyshev polynomials and series. The Chebyshev-Fourier basis has been chosen due to the underlying numerically fast discrete Fourier transforms, and also to initialize this novel method for wall-bounded incompressible Euler flow in using widely known and well studied spectral representations. The applied methods allow for fast numerical convergence and, thus, for an accurate representation of the solutions to the Helmholtz-Hodge decomposition and the flow itself in low to moderate resolutions. Although the choice for a spectral approach introduces difficulties for higher resolutions, as spectral derivatives become progressively more erroneous and unstable, we have found ways to overcome those pitfalls and find that the overall algorithm is stable, accurate, and fast. This is to a large extent due to sophisticated Poisson solvers, which solve directly for the second-order radial derivatives and divisions by the radial argument. One major advantage of shifted Chebyshev series is that they are obtained from function values known on the shifted Chebyshev extrema, namely ((cos (iπ/N ) + 1)/2) i=0,...,N . Those are clustered near (and contain) the boundary of the cylindrical domain as well as the pole r = 0. This clustering is important, as the numerical studies in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] indicate a singularity directly at the boundary. In summary, we choose to deal with unstable spectral derivatives in exchange for very good approximation properties and a high point density at the boundary and pole of the cylinder.

Decoupling of the Poisson problems.

Let in the following u = u(r, z) and w = w(r, z) be sufficiently smooth z-periodic functions [0, 1] × [0, L] → R. In sight of eqs. [START_REF] Lamb | Hydrodynamics[END_REF] and [START_REF] Lichtenstein | Über einige Existenzprobleme der Hydrodynamic[END_REF] and the fact that in our approach we will not solve eqs. [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF] and ( 45) explicitly, we will treat the following homogeneous model problems

Gu := ∂ 2 r + 1 r ∂ r + ∂ 2 z u(r, z) = f (r, z) ∂ r | r=0,1 u(r, z) = 0, (64)
and

Hw := ∂ 2 r + 1 r ∂ r -1 r 2 + ∂ 2 z w(r, z) = g(r, z), w(0, z) = w(1, z) = 0 . ( 65
)
As discussed in the discretization Section 5, we use the periodicity in vertical direction and anticipate an approximation of the solutions in truncating the Fourier expanded differential equations. In replacing the appearing functions by truncated Fourier series, i.e.

Λ M (r, z) = M/2 |k|=0 Λk (r) e i kz , Λ ∈ {u, f, w, g}, (66)
the equations decouple and for each k ∈ [0, M/2] we are left to solve the ordinary differential equations

(L 1 -k2 )û k (r) = û k (r) + 1 r û k (r) -k2 ûk (r) = fk (r), û k (0) = û k (1) = 0, (67)
and

(L 2 -k2 ) ŵk (r) = ŵ k (r) + 1 r ŵ k (r) -( 1 r 2 + k2 ) ŵk (r) = ĝk (r), ŵk (0) = ŵk (1) = 0, (68)
where L 1 and L 2 are defined in eqs. [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF] and [START_REF] Kerr | Evidence for a singularity of the three-dimensional, incompressible Euler equations[END_REF]. Since only real valued functions, whose complex modes verify ûk = û-k , are treated, we may restrict ourselves to positive modes. Note that eqs. (67) and (68) are commonly known as the modified Bessel ODEs of order 0 and 1 respectively. Their solutions, well known yet hard to compute numerically, are described later in this section.

Remark. To obtain the solution of the inhomogeneous problem eq. ( 42), with Neumann boundary data b = b(z) at r = 1, from a homogeneous problem eq. ( 64), one simply solves the latter with the modified right hand side f = f (r, z) = f (r, z) -G( r 2 2 b(z)) for a solution ũ = ũ(r, z). The inhomogeneous solution is then given by ũ(r, z)+ r 2 2 b(z). The smoothness requirement for b is neglected here, because the procedure is directly applied to the truncated problem, where all involved functions are C ∞ . The special case, where k = 0 can be treated explicitly, since the operators L, L 1 and L 2 are readily inverted. Note that L does not appear in the Poisson problems but is needed to solve directly for A (s) and its derivatives as described in Section 4. Suppose h = h(r), r ∈ [0, 1], is a generic function that is represented by a truncated shifted Chebyshev series

h(r) = N l=0 h l T * l (r) . ( 69
)
The inverses of the linear differential operators are

L -1 h(r) = 1 r r 0 t h(t) dt = N l=0 h l 1 r r 0 t T * l (t) dt , (70)
L -1 1 h(r) = r 0 1 t t 0 yh(y) dy dt = N l=0 h l r 0 1 t t 0 yT * l (y) dy dt , (71)
L -1 2 h(r) = r 0 t c(h) + t 0 h(x) dx dt r = N l=0 h l   r 0 t c(T * l ) + t 0 T * l (x) dx dt r   , (72)
where

c( . ) = -2 1 0 t t 0 ( . )(x) dx dt
ensures homogeneity at r = 1. Any bounded differential operator can then be applied to the inverses, and, thus, the derivatives of the solutions to the above ODEs can be obtained on a fixed grid via a matrix vector multiplication. For example, the solution of eq. ( 67), with k = 0, subject to the differential operator D r ∈ {∂ r , ∂ 2 r , ∂ r /r, 1/r} is approximated by

( D r û0 (x i )) i=0,...,N = ( D r L -1 1 T * l (x i )) i,l=0,...,N • ( f0l ) l=0,...,N , (73)
on the fixed grid (x i ) i=0,...,N . Here ( f0l ) l=0,...,N denote the coefficients of the shifted Chebyshev series that approximate f0 in eq. ( 67) with k = 0. In our CLA code, the matrices ( D r L -1 1 T * l (x i )) i,l=0,...,N are pre-computed for the needed derivatives and stored for continuous re-application. This method, to directly obtain the derivatives of the solution to an ODE, is very accurate and of usual complexity of O(N 2 ) for spectral solving methods with full matrices. If k > 0, then the operators L 1 -k2 and L 2 -k2 are known as the modified Bessel ODEs of order 0 and 1 respectively. Those special operators are also invertible and discussed in Section 5.3. But first we will apply to eqs. ( 67) and (68), with k > 0, well known spectral solvers and propose a slight modification, which allows us to use the here presented technique to solve for the radial derivatives directly.

Notation. The application of a bounded differential operator ∂ i r /r j to an inverse operator L -1 will be indicated by superscripts of the form L -1,(i)/r j , e.g. L -1,(2) := ∂ 2 r L -1 and L -1,(0)/r := 1 r L -1 .

Chebyshev-Galerkin solver based on a modified Galerkin method

Let us start with the homogeneous Dirichlet problem eq. (68). In [START_REF] Shen | Efficient spectral-Galerkin methods III: Polar and cylindrical geometries[END_REF], the author introduced a Chebyshev-Galerkin scheme for the homogeneous Helmholtz equation in a disc. As a matter of fact, the operators in use are almost identical to ours and, therefore, we may apply his findings to our case. A Galerkin scheme consists of a variational formulation of a given problem together with a suitable test function space. A finite-dimensional sub-space of the latter, spanned by a complete, sometimes orthogonal, set of basis functions, entails a discrete approximation of the solution to the problem. The Lax-Milgram theorem ensures existence and uniqueness if continuity and coercivity conditions on the variational form hold true. For a complete introduction on Galerkin schemes, we would like to refer the reader to [START_REF] Shen | Spectral methods: algorithms, analysis and applications[END_REF]. Here, we aim directly at a suitable formulation for Chebyshev polynomials.

Notation. The method described below is called Chebyshev-Galerkin (CG) method and abbreviated by CG Dir 1 , where the super-/subscripts indicate the Dirichlet boundary condition and the used set of basis functions (B 1 from eq. ( 80)) 4 . The Neumann case, which treats another operator, will be later referred to as CG N eu .

Starting from eq. ( 68), the radial interval is expanded via the coordinate transformation r = (t + 1)/2, t ∈ I := [-1, 1] to give

y (t) + 1 t + 1 y (t) - 1 (t + 1) 2 + k2 4 y(t) = 1 4 ĝk ((t + 1)/2), (74)
with y(t) := ŵk ((t + 1)/2), k ∈ [0, M/2] fixed, which we multiply by (t + 1) to obtain

(t + 1)y (t) + y (t) =(( . +1)y ) (t) - 1 t + 1 + k2 4 (t + 1) y(t) = h(t), (75)
where h(t) = 1 4 (t + 1)ĝ k ((t + 1)/2), k = (2πk)/L as in eq. ( 62), and y(0

) = y(1) = 0. A suitable test function space is Y 0 := H 1 0,ω (I)
, where the subscript ω indicates that the inner product in

H 1 0 (I) := {θ ∈ H 1 (I) = W 1,2 (I) | θ(-1) = θ(1) = 0 in the trace sense} (76)
is weighted by ω(t) = (1 -t 2 ) -1/2 , i.e. the Chebyshev weight from Section A. In the following let v ∈ Y 0 , we apply the inner product ( . , v) ω = I . vω dt to both sides of eq. ( 75), execute an integration by parts in the LHS and multiply by (-1) to obtain

I (t + 1)y (t)(vω) (t) dt + I 1 t + 1 y(t)v(t)ω(t) dt + k2 I 1 4 (t + 1)y(t)v(t)ω(t) dt = - I h(t)v(t)ω(t) dt . (77)
The variational formulation to our problem eq. ( 75) now reads as follows: find y ∈ Y 0 , such that eq. ( 77) holds for each v ∈ Y 0 . In defining

Y N 0 := {p ∈ P N | p(-1) = p(1) = 0} ⊂ Y 0 , (78)
with P N being the polynomial space of maximal degree N ∈ N, we may state the Galerkin scheme for our homogeneous Dirichlet problem eq. ( 75):

find y = y N ∈ Y N 0 s.t. eq. (77) holds ∀v = v N ∈ Y N 0 , (79)
in other words, we look for an approximation of the solution y ∈ Y 0 in the finite dimensional space Y N 0 . The advantage of the scheme in [START_REF] Shen | Efficient spectral-Galerkin methods III: Polar and cylindrical geometries[END_REF] (see also [START_REF] Shen | Spectral methods: algorithms, analysis and applications[END_REF]) is that it uses special basis functions that fulfill the boundary conditions exactly and allow additionally for simple structured solving matrices. Those basis functions in the Dirichlet case are given by

Φ l := T l -T l+2 , l = 0, . . . , N -2, B 1 := {Φ 0 , . . . , Φ N -2 } . ( 80
)
The set B 1 is a basis of Y N 0 and (79) must, therefore, only hold for any v = Φ 0 , . . . , Φ N -2 . Thus, we seek a solution in the form

y N (t) = N -2 l=0 a l Φ l (t), t ∈ I, (81)
which we insert into eq. ( 77) to find

N -2 l=0 a l I (t + 1)Φ l (t)(Φ i ω) (t) dt+ I 1 t + 1 Φ l (t)Φ i (t)ω(t) dt + k2 I (t + 1) 4 Φ l (t)Φ i (t)ω(t) dt = - I h(t)Φ i (t)ω(t) dt = - N j=0 ĝkj I 1 4 (t + 1)T j (t)Φ i (t)ω(t) dt , (82)
for each i = 0, . . . , N -2. In the last line we have inserted the Chebyshev approximation of g in order to approximate the integral on the RHS. In fact, since we know the values of ĝk on the shifted Chebyshev nodes, we have ĝk (r) = N j=0 ĝkj T * j (r) and, therefore, ĝk ((t + 1)/2) = N j=0 ĝkj T j (t). To complete the Galerkin scheme CG Dir 1 , we solve the following linear system of equations

(A + B -k2 C) • (a l ) l = -D • (ĝ kj ) j , (83)
where the entries of the solving matrices (Λ = (Λ i,l )) are given by

A i,l = 1 -1 (t + 1)Φ l (t)(Φ i ω) (t) dt, B i,l = 1 -1 1 t + 1 Φ l (t)Φ i (t)ω(t) dt, C i,l = 1 -1 1 4 (t + 1)Φ l (t)Φ i (t)ω(t) dt, D i,j = 1 -1 1 4 (t + 1)T j (t)Φ i (t)ω(t) dt, (84)
for i, l = 0, . . . , N -2 and j = 0, . . . , N . The above integrals integrate exactly and have to be computed only once. The entries of the first three matrices are stated in [START_REF] Shen | Efficient spectral-Galerkin methods III: Polar and cylindrical geometries[END_REF] and are not repeated here. The entries of D, which are not stated in [START_REF] Shen | Efficient spectral-Galerkin methods III: Polar and cylindrical geometries[END_REF], read

D i,j =                (δ i1 + 1)π/16, j = i -1, (δ i0 + 1)π/8, j = i, δ i0 π/16, j = i + 1, -π/8, j = i + 2, -π/16, j = i + 3,
for i = 0, . . . , N -2 and j = 0, . . . , N . The system (83) is solved for a solution in the form (81) and transformed to the solution ŵN k of eq. ( 68) in form of a shifted Chebyshev series, namely N j=0 ŵkj T * j , in using the variable transform t = 2r -1, r ∈ [0, 1] and the coefficient transform

ŵkj = a j -a j-2 , j = 0, . . . , N, where we set a -2 = a -1 = a N -1 = a N = 0.
Let us come to the Poisson problem with Neumann boundary conditions (67). The approach is identical to the Dirichlet case and, hence, abbreviated. The same procedure as before turns eq. ( 67) into

(( . + 1)z ) (t) - k2 4 (t + 1)z(t) = l(t), t ∈ I, z (-1) = z (1) = 0, (85)
with z(t) = ûk ((t+1)/2) and l(t) = 1 4 (t+1) fk ((t+1)/2), k ∈ [0, M/2] fixed. A variational formulation of the latter reads

I (t + 1)z (t)(vω) (t) dt + k2 I 1 4 (t + 1)z(t)v(t)ω(t) dt = - I h(t)v(t)ω(t) dt, (86)
for v ∈ Y := {θ ∈ H 1 ω (I) | ∂ θ(-1) = ∂ θ(1) = 0}.
We look for an approximation z N , N ∈ N, of the solution z in a finite dimensional polynomial subspace Y N := {p ∈ P N | p (-1) = p (1) = 0} ⊂ Y , that is spanned by the basis functions

Ψ l := T l -γ l T l+2 , l = 0, . . . , N -2, (87)
where

γ l = l 2 (2 + l) 2 .
Inserting the basis functions together with the representations z(t) = N -2 l=0 z l Ψ l (t) and fk ((t+1)/2) = N j=0 fkj T j (t) into eq. (86) entails the system of linear equations

(E + k2 F ) • (z l ) l≤N -2 = G • ( fkj ) j≤N , (88)
with matrices E, F ∈ R N -2×N -2 and G ∈ R N -2×N , whose entries are defined by

E i,l = 1 -1 (t + 1)Ψ l (t)(Ψ i ω) (t) dt, F i,l = 1 -1 1 4 (t + 1)Ψ l (t)Ψ i (t)ω(t) dt, G i,l = 1 -1 (t + 1) 4 T j (t)Ψ i (t)ω(t) dt .
Once eq. ( 88) is solved, we find the Chebyshev series representation of the approximated solution

ûN k = N j=0 ûkj T * j in using the coefficient transform ûkj = z j -γ j-2 z j-2 , j = 0, . . . , N, (89)
with γ j = 0 for j = -2, -1, N -1, N . It is rather tedious to find the exact values of the matrix entries, yet they can also be computed via Chebyshev-Gauss integration. This quadrature rule is accurate for polynomials with a maximal degree of 2n -1, where n is the number of nodes on which the integrands are evaluated (see eq. ( 190) in Section A). As a matter of fact, the above expressions are weighted integrals over polynomials (t + 1)

T i T j = 1 2 (t + 1)(T i+j + T |i-j| ).
Therefore, those polynomials are of maximal degree i + j + 1 and, for the application of the Chebyshev-Gauss quadrature, it suffices to evaluate the integrands on n = (i + j)/2 + 1 Chebyshev nodes (i.e. on Chebyshev zeros). The matrices have to be computed only once and to make sure that the values are correctly stored and to minimize rounding errors, the matrices in the present implementation are computed in quadruple precision and rounded to working precision. Also, the programming effort reduces to a minimum, because this method can be applied to all the preceding matrices in this section. Moreover, once computed for a high number of rows and columns, the matrices can be saved to disc in quadruple precision and read in for an extraction of any inferior or equal (Chebyshev-)dimension.

A slight modification of the Chebyshev-Galerkin solvers is undertaken in order to solve directly for the second vertical derivative of the periodic solution u = u(r, z), resp. w = w(r, z). By the help of this modification, one may, subsequently, use the methods of inverted operators from further above to solve directly for radial derivatives. Instead of solving for the modes ûk , the solution of eq. ( 67), we will solve directly for the terms -k2 ûk , which are the Fourier modes of the second vertical derivative of u. This is readily achieved in dividing the solving matrix by -k2 . More precisely, for the Neumann problem we will solve, instead of eq. ( 88),

Q • (ā l ) l = G • ( fkj ) j≤N , (90)
where Q = -(E + k2 F )/ k2 and āl = -k2 a l . Note that only the solving matrix is modified and not the RHS of eq. ( 67). Once we have solved eq. ( 90) for all k ∈ [0, M/2], the second vertical derivative of u may be back-transformed to physical space and subtracted from the RHS function to entail

L 1 u M (r, z) = f M (r, z) -∂ 2 z u M (r, z) =: f M (r, z) . (91)
The latter equation may be accurately solved for radial derivatives and divisions by means of the methods L -1,( 2)

1
and L -1,(1)/r 1 from before. Exactly the same procedure may be applied to the Dirichlet boundary problem eq. ( 68) with the associated linear system eq. ( 88). One solves with the test function u1(r) = e 1-cos(2πr) -1. The RHS function is analytically calculated as f1(r) = L1u1(r) and implemented exactly, the relative errors are computed in the l 2 -norm on the fixed grid. The legend references the graphs by the method that was used to obtain them. Here, the CG N eu method was used to solve L1u1 = f1 for the solution itself, to which standard spectral differentiation and division methods (see Section A) were applied. Those spectral differentiation and division methods, applied directly to the implemented function u1 (no solving is involved, u1 is directly projected into the finite dimensional Chebyshev-Fourier space), are referenced with spec.

-k-2 (A + B -k2 C) • (z l ) l≤N -2 = -D • (ĝ kj ) j , (92)
for zl = -k2 z l , and obtains the radial divisions and derivatives via

L 2 w M (r, z) = g M (r, z) -∂ 2 z w M (r, z) =: ḡM (r, z) . ( 93
)
Note that the term k = 2πk/L may become very large, especially for large grid sizes and small lengths of periodicity L > 0, and it may quickly dominate eqs. (67) and (68). Therefore, a backtransform of the second vertical derivative to physical space before subtraction from the RHS function is recommended. As a matter of fact, the condition number of the solving matrix is not altered by the division and, hence, the system is solved as accurately as for the solution itself.

Notation. The methods that divide the solving matrices by -k2 and, subsequently, apply the methods of inverted operators for the treatment of eqs. ( 91) and (93) will be referred to by CG Dir, (1) , CG Dir,(0)/r , CG N eu, (2) and CG N eu,(1)/r .

Figure 1 shows a comparison between the Chebyshev-Galerkin solver CG N eu from this section and the method of inverted operators for L -1 1 in the particular case k = 0. Of course, in this particular case, there is no division by -k2 in the solving matrices.

As visualized in the referenced figure, the solving method for L -1 1 produces the needed expressions close to machine precision. This implies for eq. ( 91) that one obtains the radial derivatives and division by r close to the accuracy of the formerly found second vertical derivative, which is in principle as accurate as the solution without the division by -k2 , because of the unaltered condition number.

Bessel-Chebyshev solver based on exact solutions to truncated problems

Initially, set the focus on the homogeneous Neumann Problem (64). Here, we will state the exact analytical solution for an approximated RHS and propose an algorithm that will give rise to a very precise approximation to the discrete solution and its derivatives. Moreover, as we have the analytical solutions at hand, we can get a formula for the second derivatives in a straight forward way and obtain a precise approximation of the latter. The resulting algorithm for this method is of complexity O(N 2 ) for each (relevant) Fourier mode, but its precision is striking, and efficient parallelization of the code assures the applicability of the method.

Let us depart directly from the Fourier truncated problem eq. (67), namely

(L 1 -k2 )û k (r) = û k (r) + 1 r û k (r) -k2 ûk (r) = fk (r), û k (0) = û k (1) = 0, where k = 2πk/L, k ∈ [0, M/2].
The case where k = 0, has been treated further above, so let us fix k > 0 in the following. If k > 0, then eq. ( 67) is well known as the modified Bessel ODE of order 0. Its fundamental solutions are the modified Bessel functions of order 0, namely I 0 ( k . ) and K 0 ( k . ).

Notation. In order to increase readability, for j = 0, 1, we will abbreviate

J [I j h]( k, r) := r 0 tI j ( kt)h(t)dt, J [K j h]( k, r) := r 0 tK j ( kt)h(t)dt, (94)
for a generic function h. If h = h(x) is a power function with exponent n, then we will write x n in place of h, i.e. J [I j x n ] and

J [K j x n ].
By an application of the method of variation of parameters (e.g. see [START_REF] King | Differential equations: linear, nonlinear, ordinary, partial[END_REF]), under consideration of the homogeneous boundary conditions, one obtains the complete analytical solution to problem (64) as ûk (r

) = I 0 ( kr) • J [K 0 fk ]( k, r) -K 0 ( kr) • J [I 0 fk ]( k, r) + c k I 0 ( kr) . ( 95
)
The solution can be verified by a direct insertion into the differential equation. The exact first order derivatives read

û k (r) = kI 1 ( kr) J [K 0 fk ]( k, r) + c k + kK 1 (kr)J [I 0 fk ]( k, r) (96)
û k (r) = 1 2 k -k(K 0 ( kr) + K 2 ( kr))J [I 0 fk ]( k, r) + k I 0 ( kr) + I 2 ( kr) J [K 0 fk ]( k, r) + c k +2r fk (r) I 0 ( kr)K 1 ( kr) + I 1 ( kr)K 0 ( kr) . (97)
Using eq. ( 96), one easily shows that

c k = -J [K 0 fk ]( k, 1) + K 1 ( k) I 1 ( k) • J [I 0 fk ]( k, 1)
ensures the Neumann homogeneity. Subsequently, we approximate the Fourier mode fk by a truncated shifted Chebyshev series,

fk (r) = N l=0 fkl T * l (r), (98)
and insert it into the precedingly found formulas. Since the expressions for the solution, and derivatives thereof, are linear in fk , we may pull out the sum and conveniently write

ûk (r) = N l=0 fkl N (0) ( k, l, r), (99)
where

N (0) ( k, l, r) := I 0 ( kr) • J [K 0 T * l ]( k, r) -K 0 ( kr) • J [I 0 T * l ]( k, r) + c * k I 0 ( kr), (100)
and

c * k = -J [K 0 T * l ]( k, 1) + K 1 ( k) I 1 ( k) • J [I 0 T * l ]( k, 1) .
Similarly, we find expressions for the derivatives, namely

û k (r) = N l=0 fkl N (1) ( k, l, r), (101)
û k (r) = N l=0 fkl N (2) ( k, l, r), (102)
with

N (1) ( k, l, r) := k I 1 ( kr) • J [K 0 T * l ]( k, r) -K 1 ( kr) • J [I 0 T * l ]( k, r) + c * k I 1 ( kr) , (103)
N (2) ( k, l, r) := 1 2 k -k(K 0 ( kr) + K 2 ( kr))J [I 0 T * l ]( k, r) + k I 0 ( kr) + I 2 ( kr) J [K 0 T * l ]( k, r) + c * k + 2rT * l (r) I 0 ( kr)K 1 ( kr) + I 1 ( kr)K 0 ( kr) . ( 104
)
To obtain the expression û k (r)/r, eq. ( 103) and eq. ( 101) are directly divided by r. We extend the definitions of the N (j) by the k = 0 terms in using the exact form of L -1 1 in eq. ( 71) (here in an integrated form):

N (0) (0, l, r) := ln(r) r 0 yT * l (y) dy - r 0 yT * l (y) ln(y) dy N (1) (0, l, r) := r 0 yT * l (y) dy r N (2) (0, l, r) := T * l (r) - r 0 yT * l (y) dy r 2 .
For the first derivative divided by r we define furthermore

N (1)/r := 1 r N (1) . (105)
Hence, given a Chebyshev-Fourier approximation of the RHS of (64) in the domain [0, 1] × [0, L], with modes fkl , the derivatives of the solution u M N = u M N (r, z) to the spectrally truncated Neumann problem take the form

∂ j r u M N (r, z) = M |k|=0 N l=0 fkl N (j) ( k, l, r) e i kz , j = 0, 1, 2 . (106)
Equally, the first derivative divided by r is given by

1 r ∂ r u M N (r, z) = M |k|=0 N l=0 fkl N (1)/r ( k, l, r) e i kz . (107)
Notation. These particular methods of matrix-vector multiplications to obtain ∂ 2 r u M N k (r, z) and ∂ r u M N k (r, z)/r will be referred to as the Bessel-Chebyshev (BC) method for the Neumann problem and denoted by N (2) and N (1)/r respectively. Other expressions do not need to be computed in the current implementation of the CLA.

Let us come to the Dirichlet problem (65), for which we proceed in the same manner as for the Neumann problem before. We recall the truncated problem eq. ( 68), for all k ∈ [0, M/2]:

(L 2 - k2 ) ŵk (r) = ŵ k (r) + 1 r ŵ k (r) -( 1 r 2 + k2 ) ŵk (r) = ĝk (r), ŵk (0) = ŵk (1) = 0
, where we further assume that ĝk (r) is approximated by a truncated Chebyshev series

N l=0 ĝkl T * l (r). • For k = 0, the operator L 2 = ∂ r 1 r ∂ r (r .
) integrates exactly, as already seen in eq. ( 72), to give

ŵ0 (r) = r 0 t c + t 0 ĝ0 (x) dx dt r = N l=0 ĝ0l r 0 t c l + t 0 T * l (x) dx dt r =: N l=0 ĝkl D (0) (0, l, r), ∂ r ŵ0 (r) = N l=0 ĝ0l   - r 0 t c l + t 0 T * l (x) dx dt r 2 + c l + r 0 T * l (x) dx   (108) =: N l=0 ĝkl D (1) (0, l, r), ∂ 2 r ŵ0 (r) = N l=0 ĝ0l 1 r   2 r 0 t c l + t 0 T * l (x) dx dt r 2 -c l - r 0 T * l (x) dx   + T * l (r) =: N l=0 ĝkl D (2) (0, l, r), (109)
where

c = -2 1 0 t t 0 ĝ0 (x) dx dt and c l = -2 1 0 t t 0 T * l (x) dx dt .
• For k > 0, eq. ( 68) is the modified Bessel ODE of first kind, whose fundamental solutions are I 1 and K 1 . The particular solutions to the homogeneous Dirichlet problem eq. ( 68) are, by the method of variation of parameters,

ŵk (r) = I 1 ( kr) J [K 1 ĝk ]( k, r) + d k -K 1 ( kr)J [I 1 ĝk ]( k, r), (110)
with

d k = -J [K 1 ĝk ]( k, 1) + K 1 ( k) I 1 ( k) • J [I 1 ĝk ]( k, 1)
to enforce homogeneity at r = 1. An insertion of the shifted Chebyshev representation of ĝk (r) and subsequent simplifications, which exploit the linearity of the involved integrals, give

ŵk (r) = N l=0 ĝkl I 1 ( kr) J [K 1 T * l ]( k, r) + d * k -K 1 ( kr)J [I 1 T * l ]( k, r) =: N l=0 ĝkl D (0) ( k, l, r), (111)
with

d * k = -J [K 1 T * l ]( k, 1) + K 1 ( k) I 1 ( k) • J [I 1 T * l ]( k, 1) .
The same procedure applied to the derivatives of the analytic solution produces

∂ r ŵk (r) = N l=0 ĝkl k 2 (I 0 ( kr) + I 2 ( kr)) J [K 1 T * l ]( k, r) + d * k + (K 0 ( kr) + K 2 ( kr))J [I 1 T * l ]( k, r) =: N l=0 ĝkl D (1) ( k, l, r) ∂ 2 r ŵk (r) = N l=0 ĝkl k 4 -k(3K 1 ( kr) + K 3 ( kr))J [I 1 T * l ]( k, r) + k(3I 1 ( kr) + I 3 ( kr)) J [K 1 T * l ]( k, r) + d * k + 2rT * l (r) I 1 ( kr)K 0 ( kr) + I 1 ( kr)K 2 ( kr) + I 0 ( kr)K 1 ( kr) + I 2 ( kr)K 1 ( kr) =: N l=0 ĝkl D (2) ( k, l, r) .
To obtain the expression ŵk (r)/r, eq. ( 111) is directly divided by r, such that we may define

D (0)/r := 1 r D (0) . ( 112
)
We have found the representations of the exact derivatives of the solution w M N = w M N (r, z) to the spectrally discretized version of eq. ( 65) in the form

∂ j r w M N (r, z) = M |k|=0 N l=0 ĝkl D (j) ( k, l, r) e i kz , j = 0, 1, 2, (113)
where D (0)/r replaces D (j) to obtain w M N (r, z)/r.

Notation. We will refer to the particular methods of matrix-vector multiplication to obtain ∂ r w M N (r, z) and w M N (r, z)/r as D (1) and D (0)/r respectively. We will refer to them as Bessel-Chebyshev (BC) methods for the Dirichlet problem. Other expressions do not need to be computed in the current implementation of the CLA.

Computation of the solving matrices Above we have stated the exact solutions and their derivatives, eqs. ( 106) and ( 113), to the truncated problems eqs. ( 67) and ( 68) respectively. In order to calculate them numerically, the exact values of the expressions N ( . ) and D ( . ) are needed, where ( . ) can be any of (0), ( 1), ( 2), (1)/r and (0)/r as defined above. Given a discrete set of points {r 0 , . . . , r N } ⊂ [0, 1], we aim at the computation of the 3-dimensional cubic matrices N ( . ) and D ( . ) with entries

N ( . ) kli := N ( . ) ( 2π L k, l, r i ), D ( . ) kli := D ( . ) ( 2π L k, l, r i ), (114)
for k = 0, . . . , M/2 and l, i = 0, . . . , N , where M + 1 and N + 1 are the Fourier and Chebyshev dimensions respectively. However, not all matrices have to be computed, if only a particular derivative is needed. Also, inferior orders of differentiation may be obtained by integration of a higher derivative if more orders are needed. The latter only works provided a supplementary function value in the interval is known. For instance, one could integrate the second derivative of the Neumann solution to obtain the first derivative, in knowing that the last must vanish at r = 0 to keep the operator L 1 bounded. Once the matrices N ( . ) , D ( . ) are computed, one may save them to the disc and reuse them any number of times. We will see in the following that the numerical computation of the above entries is by no means straightforward or low in complexity and a non-negligible effort has been employed to obtain them. Therefore, the needed matrices have been computed for a moderately high number of Chebyshev-Fourier modes on the nested grid {r i } 0≤i≤N of the shifted Chebyshev-extrema. This enables us to use the matrix for M/2 = N = 512 modes (513 3 entries) to solve on all (2 n + 1) × 2 m with n ≤ 9, m ≤ 10, physical Chebyshev-Fourier sub-grids of 513 × 1024 points. Note that the odd numbers 2 n + 1 have been chosen deliberately because the underlying fast Chebyshev transform on the Chebyshev-extrema grid, namely the DCT-I in FFTW3, shows the best performance for these numbers (see [START_REF] Frigo | The design and implementation of FFTW3[END_REF][START_REF] Frigo | FFTW3[END_REF] for details on the FFTW3).

The most involved part of the matrix computations is the evaluation of the integrals in (94). The modified Bessel functions expose such extreme behavior at r = 0 and r = 1 that any standard quadrature formula even in quadruple precision is insufficient to find accurate values. The functions K j , j = 0, 1, . . . possess (logarithmic) singularities at the pole r = 0 and I j , j = 0, 1, . . . grow exponentially for large arguments ( kr). Also, these integrals need to be evaluated to a very large number of digits, because the summands in the exact analytical solutions and their derivatives show violent cancelations. To be sure to calculate all entries correctly to more than 36 digits, a precision of roughly

4 5 2π
L k digits is needed, being especially costly for small heights L or large k, but the cancellations do dependent on the actual values of 2π L k, l and r so that the working precision may be adapted. The amount of needed computational power and time seems to be in no proportion to the seemingly standard goal of inverting a Calderon-Zygmund operator of order zero. For very accurate computations of incompressible Euler flow, with a possible goal of detecting finite-time singularities, however, this one-time effort might be appropriate. Because of steep singularities at r = 0 and exponential growth at r = 1, the integrals in (94) are not amenable to treatment by standard quadrature schemes and we will thus proceed in mathematical theory, which will allow for a direct evaluation of the latter integrals in terms of Hypergeometrical and Meijer-G functions.

Lemma 1. The following holds for k = 2π

L k, n ∈ N 0 and r ∈ [0, 1],

J [I 0 x n ]( k, r) = r 0 I 0 ( kt)t n+1 dt = r n+2 1 F 2 n 2 + 1; 1, n 2 + 2; k2 r 2 4 n + 2 , (115)
J [K 0 x n ]( k, r) = r 0 K 0 ( kt)t n+1 dt = r n+2 (n + 2) 2 (n + 2)K 0 ( kr) 1 F 2 1; n 2 + 1, n 2 + 2; k2 r 2 4 + krK 1 ( kr) 1 F 2 1; n 2 + 2, n 2 + 2; k2 r 2 4 , (116)
J [I 1 x n ]( k, r) = r 0 I 1 ( kt)t n+1 dt = t n+1 k I 0 ( kt) t=r t=0 -(n + 1)J [I 0 x n ]( k, r) (117) = 1 4 kr n+3 Γ n+3 2 Γ ( n+5 2 ) 1 F 2 n + 3 2 ; 2, n + 5 2 ; k2 r 2 4 , (118)
J [K 1 x n ]( k, r) = r 0 K 1 ( kt)t j+1 dt = 1 4 r n+2 G 2,1 1,3 kr 2 , 1 2 ; -n 2 -1 2 , 1 2 , -n 2 -1 (119) (n>0)
= -

t n+1 k K 0 ( kt) t=r t=0 + (n + 1)J [K 0 x n ]( k, r), (120)
J [K 1 x 0 ]( k, r) = πr(L 0 (kr)K 1 (kr) + L 1 (kr)K 0 (kr)) 2k , (121)
where 1 F 2 is a generalized hypergeometric function, G the Meijer-G function and L the Struve-L function, see [START_REF] Olver | NIST handbook of mathematical functions hardback and CD-ROM[END_REF] for a comprehensive introduction to these special functions.

Proof. The identities eqs. ( 117) and (120) are obtained by partial integration. Equation ( 121) is obtained from eq. ( 6.561-4) on page 676 in [START_REF] Jeffrey | Table of integrals, series, and products[END_REF] by setting a = kr and a change of variables. Equations (115) and (118) are adapted from eq. ( 84) on page 198 in [START_REF] Bateman | Tables of integral transforms[END_REF] again via a change of variables. Equations ( 116) and (119) originate from eq. ( 6.592-2) on page 690 in [START_REF] Jeffrey | Table of integrals, series, and products[END_REF] and references therein, see also [START_REF] Oberhettinger | Tables of Bessel Transforms[END_REF] page 365. These identities have additionally been verified by Mathematica's symbolic integration and via its arbitrary precision integration and evaluation for arbitrarily chosen values of k, r and n.

We may now use the explicit polynomial expression (183) of a shifted Chebyshev polynomial T * l (r) = T l (2r -1) = l j=0 a j,l r j , where a j,l = 1, l = 0, j = 0,

4 j l(-1) l-j (j+l-1)! (2j)!(l-j)!
, l > 0, j = 0, . . . , l, and find the explicit formulas for the integrals in (94) by inserting I 0 , I 1 , K 0 and K 1 into

J [ . T * l ]( k, r) = l j=0 a j,l J [ . x j ]( k, r) . (122)
The cubic matrices N ( . ) and D ( . ) from eq. ( 114) can then be obtained as a sum of special functions which are evaluated to a high enough precision. The working precision for these computations has to be large enough to provide for the range of the cancellations. The evaluation of the involved special functions is very costly and can only be achieved by modern arbitrary precision packages, which have the needed functions implemented. An efficient parallelization is also mandatory for moderate to high resolutions. To reduce the computational amount, eqs. ( 117) and (120) are implemented and re-use the entries of J [I 0 x n ]( k, r) and J [K 0 x n ]( k, r), n ≥ 1. Currently, eq. ( 119) is implemented for J [K 0 x n ]( k, r), n = 0, but eq. ( 121) is faster and the involved Struve-L functions are easier to find as components of arbitrary precision libraries. The Library that is used to compute the values of the solving matrices N ( . ) and D ( . ) is ARBLIB [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. This C library uses ball arithmetic and has a large number of special functions implemented. Figure 2 shows an example test and comparison of the methods CG N eu,(1)/r and N (1)/r , as well as CG Dir,(1) 1 and D (1) . The analytic test function is chosen to have several Chebyshev-Fourier modes whose modulus are above machine epsilon, it is given in the caption of the figure. All the presented solving methods prove to be stable and accurate for L = 2π. In this case the resulting relative errors almost overlap. If L = 1/6, then the accuracy levels shift higher but the methods stay stable and comparatively accurate. The blue-dotted curves in Figure 2, labeled by 'spec.' are the result of standard spectral projection of the solutions themselves with no solving involved, i.e. the relations in eqs. (191) and (192) are applied to the Chebyshev modes directly. In this test scenario, the BC methods N ( . ) and D ( . ) give the most accurate results, their relative errors stay close to machine epsilon even for L = 1/6, where k becomes large. In summary, both, the CG and BC methods are far superior to the standard spectral methods for the derivative and division in the radial argument. For further tests of the presented Poisson solving methods, we would like to refer the reader to [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF].

Remark. By the virtue of Lemma 1, the presented methods for the computation of the analytical solutions to truncated problems may also be applied to other polynomial basis functions. Instead of Chebyshev polynomials, one may use Lagrange, Laguerre, Jacobi, or any other useful polynomials in the approximation of the RHS function. Also, the choice of evaluation points may be deliberately chosen and can differ from the grid points that are related to the basis functions. The CLA code uses Chebyshev polynomials, mainly because of their connection to Fourier transforms, which are computed by the FFT algorithm.

Flow fields along the trajectories and spatial interpolation

Once sufficiently many coefficients X (s) have been computed, a suitable time-step can be inserted into the truncated time-Taylor series in order to approximate the new particle positions via

X(a, t) ∼ = S s=0 X (s) (a) t s = S s=0 R (s) (r, z)e r + A (s) (r, z)e α + Z (s) (r, z)e z t s , (123)
where S ∈ N is the truncation order. Clearly, the cylindrical components of X are approximated by their respective truncated time-Taylor series, e.g.

X r (r, z, t) ∼ = S s=0 R (s) (r, z) t s , R (0) = r, (124)
and X α , X z alike. Due to the assumed axisymmetry of the underlying flow, it is sufficient to compute the time-Taylor coefficients for particles that lie initially on a slice of the cylinder D(1, L). The discrete grid points on that slice, seen as distinct particles, disperse into the cylinder with time. Some particles may also leave the region D(1, L) since the actual flow takes place in an infinite cylinder along the vertical axis. The periodicity property allows us to re-map those leaving particles back inside the modeled region. It is important to note, that, again, due to the axisymmetric, the angular displacements of the particles from a chosen slice can be neglected and only the radial and vertical movements of the particles are relevant. Therefore, a scattered two-dimensional interpolation has to be executed to interpolate the flow fields, known only at the endpoints of the particle trajectories, onto the regular Chebyshev-Fourier grid.

The velocity field along the trajectories is given by the characteristic equation ∂ t X(a, t) = v(X(a, t), t), which reads in cylindrical coordinates

∂ t ( Re R + Ze z ) = v( R, Θ, Z, t) = v r ( R, Z, t 1 )e R + v α ( R, Z, t 1 )e Θ + v z ( R, Z, t 1 )e z . (125)
Here, we have used the representation of the characteristic map [START_REF] Ebin | Groups of diffeomorphisms and the motion of an incompressible fluid[END_REF] in its local basis

X(r, α, z, t) = R(r, z, t)e R + Z(r, z, t)e z ( 126
)
where e R = (cos Θ, sin Θ, 0) and e Θ = (-sin Θ, cos Θ, 0) with

R := (X r ) 2 + (X α ) 2 , Θ := A (X α , X r ) + α, Z := X z , (127)
and A as in eq. ( 14). With these definitions we obtain the correspondences

v r ( R, Z, t) = Ẋr cos (A (X r , X α )) + Ẋα sin (A (X r , X α )) , (128)
v α ( R, Z, t) = -Ẋr sin (A (X r , X α )) + Ẋα cos (A (X r , X α )) , (129)
v z ( R, Z, t) = Ẋz . ( 130
)
The vorticity ω(r, α, z, t) = ω r (r, z, t)e r + ω α (r, z, t)e α + ω z (r, z, t)e z is interpolated in the same way once it is obtained via the vorticity transport formula ω (X (a, t) , t) = DX(a, t) • ω 0 (a), which entails

ω r ( R, Z, t) = cos(A (X r , X α )) - X α ω α 0 r + X r r ω r 0 + X r z ω z 0 + sin(A (X r , X α )) X α r ω r 0 + X α z ω z 0 + X r ω α 0 r , (131)
ω α ( R, Z, t) = -sin(A (X r , X α )) - X α ω α 0 r + X r r ω r 0 + X r z ω z 0 + cos(A (X r , X α )) X α r ω r 0 + X α z ω z 0 + X r ω α 0 r , (132)
ω z ( R, Z, t) = ω r 0 X z r + ω z 0 X z z . (133)
The displacement inside the cylinder and the relevant radial and vertical movements inside of the slice in the cylinder are illustrated in Figure 3 (a). After a time-step insertion and the computation of the components X r , X α and X z , we seek to interpolate the components of the flow fields from the scattered Lagrangian points ( R(r i , z j , t 1 ), Z(r i , z j , t 1 )) i,j back onto the fixed grid (r i , z j ) i,j . Clearly, the difficulty here lies in the fact that one has to interpolate from an irregular grid onto a regular grid. Furthermore, particles may leave the convex hull of the Chebyshev-Fourier grid points in vertical direction or show a contraction, such that there are grid points that do not lie in the convex hull of the scattered particles. In the latter case one could perform an extrapolation, but it is better to use the periodicity of the flow and simply include more particles from above and below the first and last horizontal grid lines as visualized in Figure 3 (b) .

Several methods for the needed two-dimensional scattered interpolations have been tested, such as radial basis functions, fitting of bi-cubic splines, as well as Delaunay tesselation with subsequent linear barycentric interpolation. However, we found that the best results by far are achieved with the cascade interpolation [START_REF] Purser | An efficient interpolation procedure for high-order three-dimensional semi-Lagrangian models[END_REF]. This method was developed explicitly for the setting of Lagrangian displacements in a fluid flow and reduces the two-dimensional problem to a sequence of one-dimensional problems. Those one-dimensional problems may then be tackled by any other (scattered) interpolation method, e.g. Lagrange or B-spline interpolations, which are both implemented. We refer to [START_REF] Purser | An efficient interpolation procedure for high-order three-dimensional semi-Lagrangian models[END_REF][START_REF] Podvigina | The Cauchy-Lagrangian method for numerical analysis of Euler flow[END_REF][START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF] for details on the cascade interpolation.

For brevity, we recall the Lagrange polynomial interpolation here and refer the reader to [START_REF] Boor | A practical guide to splines[END_REF] for a an introduction to B-splines and the involved interpolation schemes. For a number n ∈ N of strictly increasing sites x 0 < x 1 < . . . < x n-1 , the Lagrange basis polynomials of order n (and degree n -1) are defined by

I l (x) := 0≤i≤n-1 i =l x -x i x l -x i , l = 0, . . . , n -1 . (134)
A function y, with function values y(x 0 ), . . . , y(x n-1 ), may be interpolated by a Lagrange polynomial as

y(x) ∼ = n-1 l=0 y(x l )I l (x) . (135)
The fact that this method does not involve a linear system, which needs to be solved, makes it very fast, and numerous implementations can be found in public libraries 5 . In the implemented Lagrange routine, the sites are chosen such that the interpolation point is in their center, if possible, i.e. a centered Lagrange polynomial is computed with the same amount of sites to one side as to the other side of the interpolation point. In vertical direction, this can always be fulfilled, because of the periodicity of the interpolated functions. In radial direction, especially close to the boundary or pole, a periodicity, or parity with respect to the interval end points, is not known and, thus, no valid extension outside of the interval [0, 1] can be constructed. In this case, the interpolation polynomial is still constructed as above with sites that lie closest to the boundary or pole. In this scenario, the Chebyshev extrema (also their zeros) are well known to reduce the Runge phenomenon [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF][START_REF] Boor | A practical guide to splines[END_REF], such that the Lagrange interpolation is stable even if an un-centered site distribution is used.

In the B-spline interpolation, one has to solve for the coefficients α i in

n-1 i=0 α i N i,p-1 (x j ) = y(x j ), (136)
where N i,p-1 are the B-splines of order p (and degree p -1) based on n + p non-decreasing knots u j , j = 0, . . . , n + p -1. The x j are the strictly increasing interpolation points from before. A B-spline interpolation module, which is based on routines from [START_REF] Boor | A practical guide to splines[END_REF], is publicly available on Github6 . Parts of this module were extracted and included in the CLA code after an adaptation to our needs. In these routines, the n + p + 1 knots u i of the B-splines from above are distributed in the interpolation interval, such that 2p knots are placed on the interval endpoints, p on each end point, and the rest is placed on the data points starting from x p/2 , if p is even, and in between the data sites, if p is odd, starting after x (p-1)/2 . See [START_REF] Boor | A practical guide to splines[END_REF] for more information on knot placing.

Time-Stepping

A central question to the CLA is how to obtain the radius of convergence of eq. (3) that limits the possible time-step in an accurate way, or, if that's not possible, how to choose an appropriate timestep. Unfortunately, there is currently no explicit formula for the convergence radius available. Thus, the choice of the time-step becomes a challenging aspect of this algorithm. In principle, it would be possible to obtain an explicit formula for the convergence radius from the constructive proof in [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF]. However, the estimates that are used in [START_REF] Besse | A constructive approach to regularity of Lagrangian trajectories for incompressible Euler flow in a bounded domain[END_REF] involve constants that originate from Schauder estimates of the Newton potentials. Those constants are in general hard to compute and slight variations in their determined value might have a strong influence on the value of the later computed radius of convergence. As a result, in the space C 1,α , α > 0, the convergence radius α > 0 for (3) can be bound by

α ≥ C ω 0 α , (137)
where C = C(Ω, α) is a positive constant that depends only on the geometry of the domain Ω and . α denotes the Hölder norm. Classically, there are other ways to approximate the radius of convergence, e.g. via the Cauchy-Hadamard formula 1 = lim sup s→∞

X (s) 1 s , (138)
or via the ratio rule

= lim s→∞ X (s-1) X (s) , (139)
if the limit exists. The norms in the above expressions can be those of C 1,α , but also of H s , s > 5/2 or simply of an L p space, in short, from any space where the convergence needs to be established. The Cauchy-Hadamard formula eq. ( 138) is less useful in numerical computations, because of the strong growth of the coefficients X (s) and the difficulty that arises when computing the root of large orders s. A straight forward attempt to use the ratio rule can be made in computing a large number S of time-Taylor coefficients and use extrapolation techniques to approximate the limits in eq. ( 138) or ( 139). This method is emphasized here for a stationary flow from [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF]. The regarded stationary vorticity in the periodic cylinder D(1, 2π) reads

ω r (r, z) = c 2 1 + 1 sin(z)J 1 (rc 1 ) , ω α (r, z) = c 2 1 + 1 cos(z)J 1 (rc 1 ) , ω z (r, z) = c 1 c 2 1 + 1 cos(z)J 0 (rc 1 ) ,
where r ∈ [0, 1], z ∈ [0, 2π] and c 1 is the first positive root of the Bessel function J 1 . The coefficients X (s) , 0 ≤ s ≤ 100 for this flow have been computed via the recursion relations stated in Section 3. Figure 4 shows an application of the ratio rule, which seems to suggest a bound of the convergence radius in L ∞ around 0.4. At this point, we would like to remind the reader that a time-step below an approximated radius of convergence stays valid regardless of the underlying resolution, because this radius is independent of resolution. These computations are, however, not very practical, because many coefficients need to be computed accurately, which slows down the CLA.

Here we would like to present a different approach, which is based on a self-adapting time-stepping scheme. In this method, the convergence radius is not computed in the running program, but only once in the beginning to give a valid initial time-step. More precisely, let the truncated time-Taylor series of X (s) , evaluated at t 1 , be abbreviated by

b s (t 1 ) = s s=0 X (s) t s 1 ( 140
)
and define

F (s, t 1 ) := b s-1 (t 1 ) -b s (t 1 ) ∞ b s-1 (t 1 ) ∞ ( 141
)
where b s denotes the temporal derivative of b s . Then after each time-step insertion the following property is checked:

ε m > F (S -1, t 1 ) > F (S, t 1 ), (142)
where ε m is machine epsilon, i.e. ∼ 10 -16 in double precision computations. In case this property is not fulfilled by a certain time-step t 1 , then this time-step may be adapted accordingly. Condition ( 142) is meant to assure that additional terms of order higher than the truncation order S would not alter the result of the truncated time-Taylor sum in the working precision. Surely this approach is only valid, if the coefficients show a point-wise monotone growth, which must be checked before hand, i.e. by computing a larger number of coefficients in certain iteration intervals. In fact, this method also accounts for alternating signs in the entries of X (s) (a), which is internally given by a two dimensional array. The derivative b s is used in definition eq. ( 141), because it represents the approximation of the velocity field by the characteristic equation, and it implies that the same property (eq. ( 142)) holds with b s instead of b s in the definition of F (s, t 1 ).

The method adjusts a suitable and initially fixed time-step at the event of an insufficient number of coefficients. But it can also be used to implement a fully adaptive time-stepping scheme. In this scenario, with a fixed truncation order S, a time-step is set to fulfil

F (S, t 1 ) ∼ = ε m . ( 143
)
This can be achieved by slowly reducing an initial value of t 1 until eq. ( 143) is met. In this scenario, higher numbers of S imply a larger time-step due to the growth behavior of the coefficients: It holds X (s) ∞ < cst./ s for s large enough, where is the minimal radius of point-wise convergence. The time-step is chosen such that the time-Taylor sum and its first temporal derivative, truncated to S coefficients, is not altered by the last terms. This implies a bound on the last term in the Taylor sum of the form |X (S) |t S < ε, for some ε > 0. Thus, the bound from before allows choosing t closer to for larger S.

Conserved quantities and testing of the implementation

This section is intended to recall the de facto standard methods for validating numerical simulations of the (axisymmetric) incompressible Euler equations and adds original methods for the CLA. To date, there is no reliable method to verify the validity of a non-stationary incompressible Euler flow. There are, however, many propositions to approach this problem. A non-exhaustive list of numerical validation techniques can be found in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF]. The most prominent conserved quantities to watch in an incompressible Euler simulation are the total kinetic energy7 (E) the total helicity (H), as well as the total angular momentum (I 1 ) and total angular helicity (H 1 ). They are given respectively by the following expressions .

E = 1 2 1 0 L 0 ((v r ) 2 + (v α ) 2 + (v z ) 2 )r drdz (144) H = 1 0 L 0 (v r ω r + v α ω α + v z ω z )r drdz (145)
I 1 = 1 0 L 0 v α r 2 drdz (146)
H 1 = 1 0 L 0 v α ω α r drdz . ( 147
)
If the velocity field v solves the incompressible Euler equations and ω is the corresponding vorticity field, then the above quantities are conserved for all times. The converse is clearly not true. The check for this conservation can, therefore, only serve as a control mechanism to verify that the computed solution is not too far from the actual solution. Furthermore, the conservation of such quantities should not be trusted too much. Indeed, tests with stationary solutions show that the relative energy loss does not correlate well to the actual relative errors in the flow fields, which could give a false impression on the validity of the computation.

In the literature, the majority of Euler flow simulations are carried out on an Eulerian grid. The Lagrangian viewpoint allows in return to carry out additional quality measurements regarding an implementation. Note, however, that since we discuss a semi-Lagrangian algorithm, we only have access to a set of characteristics departing from a fixed grid within one time-step. This set of characteristics is different from one time-step to another. The following validation criteria are thus only applicable within one time-step and not with respect to the initial state (as it would be in a purely Lagrangian scheme without interpolation). One of those criteria is the movement of the boundary particles. Theoretically, they stay on the boundary, but numerically, an error in the distance to the boundary after every iteration can be measured. In our cylindrical domain, this is easily achieved in computing max

z∈[0,L] |1 -R(1, z, t 1 ) 2 + A(1, z, t 1 ) 2 |, ( 148
)
where t 1 is a time-step. In the aftermath, this error is usually very close to the error in the energy conservation, even though the current implementation of the CLA resets by default the boundary value to 1 after every iteration.

Unique to the implemented Cauchy-Lagrange formalism is also the possibility to verify directly the Cauchy invariants equation (CIE), eq. ( 5), together with the incompressibility condition eq. ( 6) after every time-step insertion. Instead of checking only some necessary conditions for a solution, this method verifies directly the CIE and thus, the Euler equations. On the contrary, in a purely Eulerian setting, the verification of the Euler equations appears to be tricky and impractical. Note that the initial vorticity in eq. ( 5) is replaced with the current vorticity after the time-step on the grid, hence, the CIE is not verified with the initial vorticity at t = 0, but with that at time t n-1 , as mentioned before.

Another conserved quantity is the circulation about a closed material curve C that moves with the fluid ( [START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Bennett | Lagrangian fluid dynamics[END_REF]). Kelvin's circulation theorem states that

∂ t Γ (t) = 0, (149)
with

Γ (t) = X(C,t) v(x, t) • dσ(x) = 1 0 v(X(C(s), t), t) • ∂ s X(C(s), t) ds, ( 150
) if C is parametrized by ∈ [0, 1], such that C(0) = C(1)
. Differently formulated, the quantity

C v(a, t) • dσ(a) (151)
is constant under the material derivative

D t = ∂ t + (v • ∇).
Numerically, this conservation is difficult to verify for arbitrary closed curves, but is in principle possible with spectral representations and the ability to evaluate them anywhere in the covered domain. An Evaluation of a non-trivial closed curve in the cylindrical domain D(1, L) would need to be done very thoroughly in order to let the integral converge accurately to its conserved value. A considerable numerical complexity and large errors seem to be implied, which practically rules out this circulation method for general curves in ambitious Euler flow simulations, i.e. with the purpose of finding finite time singularities. In return, for simple circular curves around the axis, this task is readily feasible. Such a curve is of the form

re r + ze z , for fixed r ∈ [0, 1], z ∈ [0, L], (152)
and is parametrized by α ∈ [0, 2π]. With the notation of eq. ( 126), Γ in cylindrical coordinates takes the form

2π 0 v( R, Θ, Z, t) • ∂ α ( Re R + Ze z ) dα = 2π 0 v( R, Θ, Z, t) • Re Θ dα (153) = 2π 0 v r ( R, Z, t)e R + v α ( R, Z, t)e Θ + v z ( R, Z, t)e z • Re Θ dα = 2π R v α ( R, Z, t) . ( 154
)
Kelvins circulation theorem tells us that Equation (154) stays constant, which can be easily verified in every time iteration via

rv α (r, z, t i ) -R(r, z, t i+1 )v α ( R(r, z, t i+1 ), Z(r, z, t i+1 ), t i+1 ) l 2 (D # ) rv α (r, z, t i ) l 2 (D # ) , ( 155
)
where t i is the time after i iterations and D # the underlying grid. This verification method is quite simple, low in numerical complexity and implemented in the CLA code. Note that in eq. ( 155) one cannot use the initial velocity v α 0 = v α 0 (r, z, 0), because, in the CLA implementation, R(r, z, t i+1 ) is the spatial displacement with respect to the underlying grid within one time-step and not within the whole integration period [0, t i+1 ]. In order to compare to the initial velocity v α 0 , one would need to track the particle positions ( R(r, z, t), Z(r, z, t)) from t = 0 to the time where the conservation of eq. ( 155) is checked.

Testing of the CLA

Tests on Stationary Solutions

The stationary velocity field, which is the principal test flow for the CLA and taken from [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF], reads component-wise

v r (r, z) = κJ 1 (c j r) sin(κz), v α (r, z) = (c j + κ 2 )J 1 (c j r) cos(κz), (156)
v z (r, z) = c j J 0 (c j r) cos(κz) .

The J k denote the k-th order Bessel functions of first kind, c j is the j-th positive root of J 1 , and κ can be any natural number. In taking the cylindrical curl of the velocity field, one obtains the corresponding vorticity components, namely

ω r (r, z) = κ c 2 j + κ 2 J 1 (rc j ) sin(κz), ω α (r, z) = c 2 j + κ 2 J 1 (rc j ) cos(κz), (157)
ω z (r, z) = c j c 2 j + κ 2 J 0 (rc j ) cos(κz) .

See Figure 6 for an illustration of this initial vorticity field with two different values of κ and roots c j .

One of the advantages of knowing such stationary solutions, is the fact that we are able to compute analytically the time-Taylor coefficients to some order. The details on this method are stated in Section B.

In the following, if not stated otherwise, then the relative errors in vector-fields are computed on the underlying Chebyshev-Fourier grid D # in the l 2 -norm. More precisely, for a two-dimensional array g, which is approximated by g, we compute

g -g l 2 (D # ) g l 2 (D # ) . ( 158
)
The l 2 -norm returns a scalar with the square root of the sum of all elements squared, and if g equals zero, then only the difference g-g is measured with the given norm (absolute error). The relative error in vector fields with higher dimensions is defined as the maximum over the l 2 -relative errors in all the (cylindrical) components. Figure 7 shows the relative errors in the first nine time-Taylor coefficients of X, belonging to the stationary flow from above, and computed via the recursion mechanism stated in Sections 3 and 4. The figure shows that the coefficients are very accurately computed for sufficiently large grid sizes. The effect of the resolution becomes further apparent in Figure 8, where the relative errors in the analytically known velocity field is visualized. Even low resolutions deliver good results due to the approximation qualities of the Chebyshev-Fourier basis.

From our investigations in Section 7 (and in particular Figure 4) we assume a valid time-step for the given stationary flow below 0.4. It was further analyzed in Figure 5 that a time-step should lie even below 0.1 to assure accuracy in a double precision computation. For a larger time-step there would be many more coefficients needed, but the Real64 number format in the arrays is ultimately not able to contain large order coefficients X (s) , as their norms grow steeply with the order s. This 158)) in the computed coefficients X (1) , . . . , X (9) , which belong to the flow in eq. ( 156) with κ = 1 and c1. The Poisson problems were solved in double precision via the Chebyshev-Galerkin solver (left) and in quadruple precision via the Bessel-Chebyshev solver (right). The true coefficients were obtained via symbolic differentiation inside a computer algebra system by the method described in Section B. limitation can be observed in Figure 5 as a sudden drop of values in the ratio formula. See in particular the relation between the critical orders and the employed working precision in Figure 5. Therefore, we have tested the time-steps 0.1, 0.05 and 0.02 for a truncation order S = 25. The time-step is maintained as long as eq. ( 142) is satisfied. If the errors in the computation grow such that eq. ( 142) is not verified anymore, then the time-step is halved. The results are illustrated in the left panel of Figure 9. We clearly see that larger time-steps yield better results, which is a common property of semi-Lagrangian schemes and stems from the fact that fewer interpolations are carried out to arrive to a certain point of time. In the right panel of Figure 9 we show the dependence of the relative errors in the velocity field to the order of B-spline functions inside the cascade interpolation. Lower to moderate orders produce the best relative errors over the majority of the integration time span, also if Lagrange interpolation is used instead of B-spline interpolation. This can be explained by the amplification of errors, namely that high order B-splines, or polynomials in the Lagrange interpolation, amplify rounding and approximation errors. Stationary solutions are good candidates for testing the CLA, because one can directly compare the true relative error in the flow fields to other validation criteria, for instance, those that were discussed in Section 8. If a good agreement between the relative error in a flow field and a relative change in a conserved quantity is found, then this correlation might help to assert the quality of a simulation once the flow fields are not known analytically, i.e. for non-stationary flows. In Figure 10, we observe that the relative changes in most of the classical conserved quantities provide poor correlations to the actual relative errors in the flow fields. Only the total angular momentum (or first-order Casimir) aligns roughly with the error in the velocity field. A better alignment is found in the preservation of Kelvin circulation along perfect circles in planes perpendicular to the vertical axis, see Figure 11.

Other validation criteria that are particular to the CLA, such as the movement of the boundary particles, the Cauchy-invariants formula (CIF), or the incompressibility condition, do neither show good indications of the validity of a simulation. The limiting factor here is that those expressions are checked only within one time-step and not in the whole integration period. The initial vorticity is replaced by that of the proceeding time-step before a new recursion cycle starts. To verify the CIF with the initial vorticity, one would also need to know the characteristics since the starting time and the derivatives thereof. Such verification is not yet implemented in the current CLA code.

Tests on a non-stationary flow without swirl

Another possibility to further test the CLA is given by swirl-free flows, whose angular component stays exactly zero at all times. As the here discussed swirl-free flows are non-stationary, one cannot directly measure a relative error in the flow fields, but fortunately they are better understood compared to general three-dimensional flows. For instance, it is shown in [START_REF] Shirota | Note on global existence for axially symmetric solutions of the Euler system[END_REF] that the axisymmetric Euler equations in a (wall-bounded) cylindrical domain Ω are well-posed globally in time for swirl-free initial data in H s (Ω) for s ≥ 3. It is shown furthermore that s > 5/2 is sufficient in case Ω = R 3 . See also [START_REF] Ukhovskii | Axially symmetric flows of ideal and viscous fluids filling the whole space[END_REF][START_REF] Raymond | Remarks on axisymmetric solutions of the incompressible Euler system[END_REF][START_REF] Shirota | Note on global existence for axially symmetric solutions of the Euler system[END_REF][START_REF] Majda | Vorticity and Incompressible Flow[END_REF] for connected results in R 3 with s > 7/2. Moreover, we can verify numerically, next to the fact that the angular component stays zero at all times, that the pseudo-vorticity ω α /r stays constant along the trajectories and verifies a maximum principle (see e.g. [START_REF] Majda | Vorticity and Incompressible Flow[END_REF]). More precisely, in the notation of eq. ( 126), we have ω α ( R(r, z, t), Z(r, z, t), t)

R(r, z, t) = ω α (r, z, 0) r . (159)
In order to verify these facts, we simulate a swirl-and divergence-free flow with initial components

v r 0 (r, z) = 2π L (1 -r) n r a+2 cos (2πz/L) , v α 0 (r, z) = 0, v z 0 (r, z) = -(1 -r) n-1 r a+1 sin (2πz/L) (a -r(a + n + 3) + 3),
where a ≥ 1 and n ∈ N. As ω = ∇ × v, and in using eq. ( 16), the radial and vertical component of the vorticity field vanish identically and the initial angular vorticity component reads

ω α 0 (r, z) = 1 L 2 (1 -r) n-2 r a sin (2πz/L) L 2 ((a + 1)(a + 3) + r 2 (a + n + 1)(a + n + 3) -r(2a 2 + 2a(n + 4) + 5n + 6)) -4π 2 (r -1) 2 r 2 .
As a matter of fact, there holds

v α (r, z, t) = ω r (r, z, t) = ω z (r, z, t) = 0 (160)
for all times t. These identities are exactly preserved by our implementation of the CLA, which is visualized exemplary in Figure 12, right panel. Therefore, all the conserved quantities which only depend on those components, such as Kelvin circulation etc., stay exactly zero in the computation. Also, the helicity evaluates exactly to zero, since for its computation v r ω r + v α ω α + v z ω z = 0 is integrated over the (r, z)-plane. As for the stationary solution in Section 9.1, the results for the swirlfree flow depend on resolution, the interpolation order, and the chosen time-step. Higher resolutions, lower to moderate interpolation orders, and larger time-steps increase the integration period until a critical energy loss occurs. The time-step t = 0.1 was maintained over the whole simulation and the run was stopped, when energy conservation exceeded 10 -6 . The total kinetic energy, in Figure 12 left panel, seems to be well preserved until t = 1.5. However, we know from the tests with the stationary solution (see Figure 10 in particular), that the actual error in the flow fields can behave differently to the error in the energy conservation. The min/max principle, in Figure 12 right panel, is well demonstrated too, even though a final shift in the min/max values on the fixed grid of around 10 -5 is observed. This shift, however, might simply stem from the fact that the real extrema do not lie on the grid and move towards, or away from the closest grid points. Higher resolutions yield better results also with respect to this point, although very large resolutions imply other inaccuracies, such as accumulating rounding errors.

Application to a potentially singular flow

In [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], the authors compute an axisymmetric and incompressible Euler Flow in the lower quarter of the cylindric region D(1, 1/6), and claim a finite time singularity of the maximal vorticity modulus for the initial swirl-only velocity field with v r 0 = v z 0 = 0 and

v α 0 (r, z) = 100 r e -30(1-r 2 ) 4 sin(2πz/L), L = 1/6 . (161)
A 3×10 8 -fold increase of the quantity ω L ∞ is reported in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], where the maximum of |ω| is attained on the boundary ∂ D(1, 1/6) (at r = 1) and more precisely, on the symmetry axis at z = 0. Their investigations depart from a totally different formulation of the incompressible Euler equations than that in the CLA. In using a transform of the stream-vorticity formulation of the incompressible Euler equations, one obtains

v 1,t + v r v 1,r + v z v 1,z = 2v 1 ψ 1,z , (162)
ω 1,t + v r ω 1,r + v z ω 1,z = (v 2 1 ) z , (163)
-∂ 2 r + (3/r)∂ r + ∂ z ψ 1 = ω 1 , (164)
with

v 1 = v α /r, ω 1 = ω α /r, ψ 1 = ψ α /r, (165)
where ψ α is the angular component of the vector stream function, satisfying -∆ ψ = ω = ∇ × v. The sub-scripts, different from 1, indicate partial derivatives as usual. See [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF][START_REF] Majda | Vorticity and Incompressible Flow[END_REF] for further details on the construction of the above formulation. The radial and vertical components of the velocity field are then given by v r = rψ 1,z , v z = 2ψ 1 + rψ 1,r .

The Euler equations in this form are numerically convenient, because they avoid divisions by the radial argument r. Note that the CLA obtains the terms A 1 and R 1 in eqs. ( 53) and ( 59) also without divisions by r in its effective implementation.

In [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], eq. ( 164) is solved via a 6th order B-Spline based Galerkin solver, and a 4th order Runge-Kutta time-stepping scheme is subsequently applied to eqs. ( 162) and (163); the appearing derivatives are obtained by a 6th order centered difference formula. However, the most important part of their method is an adaptive moving mesh, which keeps shifting grid points into the region with peak vorticity values. In [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], this adaptive mesh allows to reach an effective maximal resolution of around (10 12 ) 2 per unit square near the location of the potential singularity. Another important facet of their implementation is the use of the preserved symmetry properties of the initial velocity field, which allows to only model the flow in the quarter cylinder D(1, L/4) and start with a much higher point density. The CLA has not yet been adapted to exploit those flow symmetries, it always computes the flow in the full cylindrical domain D(1, L), here with L = 1/6. This can be seen as a draw-back of our current implementation because a much higher grid size has to be employed to obtain similar initial local resolutions. However, even with this drawback, the computational results are convincing, which highlights the potential of the CLA as a whole.

As in Section 9, we will investigate different parameter settings to show their influence on the outcome of a simulation. In all the presented runs, the computations were aborted if the time-step dropped below 5 × 10 -8 . The time-step is halved, whenever the precision criterion in eq. ( 142) is not fulfilled at the truncation order S, i.e. if F (S, t) > ε m for a time-step t. This halving takes place until criterion (142) is matched or the code aborted due to insufficient energy or helicity conservation 8 . The results here are mainly measured by means of the conservation of helicity, energy, and Kelvin circulation. Energy and helicity have not proven to be particularly good indicators for correctness, at least for the stationary solution examined in Section 9.1, but they are very common in computational fluid dynamics, and in using them, the results become more comparable to other numerical methods. However, if the limit of 10 -5 in relative energy or absolute9 helicity change is reached, then the solution is likely to be inaccurate. Fortunately, some maximal vorticity values ω ∞ at specific time instances are stated in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] and are, therefore, directly comparable. The list of those maximal values is given in Table 1 and is appended by an additional value taken from [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF]. The author of [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF] has repeated [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] and [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF]. They will serve as orientation points for our simulations and are marked as black crossed in the maximal vorticity plots. At least three algorithms [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF][START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF][START_REF] Kolluru | A pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation: tygers and thermalization[END_REF] have good agreement on the value at t = 3.1 • 10 - the simulation of the same equations as in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], i.e. system (162)-( 164), with a spectral Chebyshev-Fourier solver on a fixed grid with 769 radial and 2048 vertical points and a Runge-Kutta time-stepping scheme (for eqs. ( 162) and ( 163)). It is important to note that [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF] also models the equations on the quarter cylinder D(1, L/4). It is further stated in [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF] that one additional maximal vorticity value was communicated by one of the authors of [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF]. This additional value at t = 0.0031, given in the third column of Table 1, is verified in [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF] with a relative error of 1.9 × 10 -8 . Also very recently the authors of [START_REF] Kolluru | A pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation: tygers and thermalization[END_REF] develop a Chebyshev-Fourier pseudospectral method with a fourth-order Runge-Kutta timestepping scheme to solve equations ( 162)-( 164), and to study the Luo-Hou initial condition [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF]. The results of [START_REF] Kolluru | A pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation: tygers and thermalization[END_REF] confirm the simulations of [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF][START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] at least until the time t = 0.0031. In addition, using an analyticity-strip method and an extrapolation method based on a linear least-square regression, the authors of [START_REF] Kolluru | A pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation: tygers and thermalization[END_REF] obtain an estimate of the time at which the (potential) singularity occurs, and this Here and in the following plots, the black crosses in the vorticity plot are the exact coordinates from Table 1 at times t = 0.003, t = 0.0031 and t = 0.0034. Also the helicity panel displays absolute errors, because the initial helicity vanishes for this flow.

time value is consistent with that of [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF]. Therefore, there are at least three independent simulations that agree on the leading digits of a maximal vorticity value at a time instance not too far away from the reported temporal singularity. In the following plots of the maximal vorticity, the second, third, and fourth column values in Table 1 are indicated by black crosses to show qualitatively how far the CLA can approach them in those test runs with various parameter settings. These approximations are given quantitatively in Table 2.

Figure 13 shows the behavior of the error in total kinetic energy, helicity, and Kelvin circulation, as described in Section 8, as well as the vorticity growth. As seen in the figure, the result strongly depends on the chosen resolution, when all other run-time parameters stay fixed. The curves of the Kelvin circulation error, which will henceforth only be referred to as the Kelvin conservation, allow to easily distinguish the low and moderate grid sizes. The helicity conservation also shows noticeable differences, while the energy conservation curves are less clearly attributable to the different resolutions. Also, all curves show very similar behavior with respect to their measured quantity. The appearing sudden drops in energy and helicity conservation are usually present, but their cause is not yet known to the authors of this work. The maximal vorticity values that were computed on the largest grid agree qualitatively with the first values in Table 1, but their exact fulfillment will be studied later together with runs in higher resolutions. Note that even for the largest grid in Figure 13, only 256 vertical grid points lie in the modeled region D(1, L/4) of [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] and [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF]. For all the curves in Figure 13, the Bessel-Chebyshev solver (BC) was used, but the Chebyshev-Galerkin (CG) method produces very similar results. Regarding the second and fourth panel of Figure 13, the maximal vorticity seems to fall off from its strong growth after a sudden change in the Kelvin conservation. This could indicate that the breaks in the Kelvin conservation and the fall-off in the vorticity growth are connected and yield inaccurate computations from these points on.

Concerning the time-steps, the method that was used to obtain the approximate convergence radius for the stationary solution in Figure 4 fails for the given initial velocity field in eq. ( 161). The initial time-Taylor coefficients X (s) cannot be accurately represented for s 20 even in quadruple precision, due to a strong growth of their l ∞ -norms. The early development of the quotient X (s-1) /X (s) seems to point to 10 -4 , see [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF]Fig. 5.4]. However, time-steps between 5 × 10 -6 and 5 × 10 -5 have proven to provide good results. The time-step dependence of the outcome is plotted in Figure 14. Remarkably, the absolute helicity change is lowest for the largest time-step, while the conservation of kinetic energy and Kelvin circulation seems better with the smallest time-step. Also, if we consider the vorticity growth, then large time-steps and their corresponding vorticity growth falls off later than the growth with smaller time-steps. When searching for sudden breaks in the Kelvin conservation, then it appears that those happen later for larger time-steps, and the same is true for the energy conservation curves.

The breaks from a rather regular behavior of the conservation curves seem to indicate a loss of accuracy and, thus, we interpret a later appearance of those breaks as a sign for better stability. Moreover, the runs with larger time-steps take a fraction of the total execution time compared to smaller time-steps, while still providing adequate, if not better, results. The latter property is quite common among semi-Lagrangian methods and helps us here to achieve very large resolutions (see further below). The initial time-steps are not maintained throughout the computation, they are adapted to allow for precision and conservation criteria as mentioned above. However, these adaptations do not happen very often and one time-step can usually be maintained over many iterations. The adaptive time-stepping procedure discussed in the end of Section 7 is examined in Figure 15 for the truncation orders S = 8, 10, 12, 15. Again, higher truncation orders, thus larger time-steps, show a break out in the Kelvin conservation slightly later and seem to maintain the strong vorticity growth slightly longer than their lower order competitors. The helicity conservation curves are very close to each other, but the high orders give lower absolute changes than the small orders. The energy and Kelvin conservation curves behave seemingly opposite to the helicity behavior. There, the relative errors are lower for lower truncation orders, but those values only depend on the velocity field and are measured by a quotient, not a difference.

Extensive tests have been carried out, involving higher working precision, different interpolation orders, the use of de-aliasing techniques, as well as refinement strategies in the interpolation step. We would like to refer to [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF]Chapter 7.2] for illustrations and discussions of those results. In summary, a higher working precision does not necessarily lead to better results 10 . The reason for that simply lies in the fast evolution of the turbulent behavior of the flow. Ever smaller Fourier-scales fill up rapidly and the flow quickly becomes under-resolved, as seen in Figure 16. Such an under-resolution cannot be dealt with by increasing the working precision, but only by increasing the resolution. Regarding the interpolation orders, lower to moderate orders give better performance, as for the stationary flow, due to lesser error amplification. De-aliasing methods 11 seem not to improve the results and refinement methods in the interpolation fall short also because of the under-resolution.

Clearly, the strongest influence is given by the resolution and will be further investigated. Figure 17 shows a comparison between a few higher resolutions. The largest grid size, which has been plotted, is 1025 × 10240. The grid spacing at the boundary of this flow is 2.35 • 10 -6 in radial direction (due to the clustering of the Chebyshev extrema) and 1.63 • 10 -5 in vertical direction. Our mesh spacings are of course still very far from the mesh spacings that were achieved in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] close to the reported singularity, but they are comparable to those in [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF]. For higher resolutions, the curves belonging to conserved quantities stay longer in a regular shape and break their behavior later, when compared to lower resolution runs, which is also visible in Figure 13. The later those breaks appear, the better is the approach of the maximal vorticity to the value at t = 0.0034 given in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] and recalled in Table 1. This qualitative observation can be examined more closely, together with other values from Table 1, by comparing the exact values in the following Table 2 Even larger grid sizes and runs with a higher working precision at those high resolutions are possible. However, larger grid sizes or a higher working-precision certainly increase the memory consumption of the CLA. The current implementation is particularly memory demanding, because all the components of the coefficients, as well as their spatial derivatives and divisions by r, up to the truncation order S, are stored for direct insertion into the recurrence relations. In return, this allows us to avoid numerical errors from standard (spectral) differentiation schemes. Surely, for the particular flow that is investigated here and in the given references, namely a blow-up candidate for the incompressible Euler equations, it is much better and more efficient to exploit the underlying symmetries, as in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] and [START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF], which allows modeling the flow in only a quarter of the periodicity domain. In the case of the CLA, such symmetry considerations would cut the memory demands tremendously and would thus allow to achieve highly resolved Euler simulations in a fraction of the time of CFL bound algorithms for the same problem.

As the presented work gives a first realization of the CLA in a wall-bounded domain, it was intended to provide a code that can simulate any axisymmetric, periodic flow, which does not necessarily possess special symmetries. Figure 18 displays the profiles of the flow field components at r = 1.

The vertical derivatives of the angular and vertical velocity components seem to become singular. The radial velocity component should stay zero at r = 1 for all times, due to the no-slip boundary condition, but errors build up especially at the symmetry axes of the flow. Better preservation of the no-slip boundary condition may indicate a higher validity of the simulation and is tracked in our implementation of the CLA. Here, the evolution of the components is only displayed until t = 3.0 × 10 -3 , the growth of the local extrema is too strong to be conveniently visualized. For surface plots of the flow field components, we would like to refer the reader to [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF].

Conclusions

In this work, we have implemented and examined the Cauchy-Lagrange algorithm (CLA) for an axisymmetric flow in a cylindrical domain with boundaries. We have systematically investigated several methods for the resulting Poisson problems, the choice of suitable time-steps, interpolation of the flow fields from the distorted Lagrangian grid back onto the regular Eulerian grid, as well as validation criteria for the computed flow. Special attention has been given to the Calderon-Zygmund operators of order zero, namely the second-order derivatives of an inverse Laplacian. The accuracy in those derivatives of the solutions to the Poisson problems, which are boundary value problems, have strong implications on the overall quality of a simulation with the CLA. In this respect, we have avoided error-prone standard differentiation methods in spectral space and developed effective techniques to directly solve for the needed derivatives. Also, common to other semi-Lagrangian methods, the interpolations of the flow fields introduce larger errors and, thus, larger time-steps yield better results, as fewer interpolations have to be carried out. In the CLA, the time-steps are only bound by the convergence radius of the time-Taylor series of the trajectories, and, therefore, they are independent of any CFL condition. This means that, for a given initial velocity field, we can use a certain time-step, which is a fraction of the convergence radius, regardless of the underlying mesh. However, due to the finite working precision in a computation, 13 the time-steps cannot be chosen too close to the radius of convergence, as demonstrated in Figure 5. In our tests, time-steps with the size of around one fourth of the convergence radius in a double precision computation produced stable and accurate results. Larger fractions of the convergence radius for a time-step may be achieved, using a higher working precision. Our current implementation of the CLA allows to choose between double, extended double, or quadruple precision, whenever the Bessel-Chebyshev (BC) solver is selected. Unfortunately, the BC method is currently restricted to 513 × 1024 Chebyshev-Fourier points with reasons explained in Section 5.3. The Chebyshev-Galerkin (CG) method only works in double precision, but it can also be adapted to work in multiple precision modes. 14Regarding the simulation of the potentially singular flow treated in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], we have seen that a high resolution is crucial to catch the strong small-scale growth of the vorticity field. Using high resolutions, we were able to achieve good agreement with the maximal vorticity values given in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF] at times, which are not too far from the potential singularity. Currently, our CLA implementation only works on a fixed product grid, but a (locally) adaptive mesh, as used in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF], seems to be mandatory to reproduce their extreme maximal vorticity values. Furthermore, the CLA results presented here compute the flow in the whole cylinder D(1, L), instead of in the quarter cylinder D(1, L/4) (as in [START_REF] Luo | Toward the Finite-Time Blowup of the 3D Axisymmetric Euler Equations: A Numerical Investigation[END_REF][START_REF] Barkley | A fluid mechanic's analysis of the teacup singularity[END_REF]), which can be seen as a draw-back in computing this particular flow. Future work may consist of restricting the CLA to the quarter cylinder in combination with a locally adaptive mesh, where grid points are shifted into the region of high vorticity values. For this, however, the Chebyshev-Fourier series representation of the flow-fields must be dropped, which could result in a slower algorithm due to the non-applicability of fast Fourier transforms. A solution, like the one performed in [START_REF] Besse | Adaptive multiresolution semi-Lagrangian discontinuous Galerkin methods for the Vlasov equations[END_REF], could be to replace the Chebyshev-Fourier representation by a multi-wavelets representation to obtain a local adaptive multi-scales representation of the solution on a non-uniform dyadic grid. Indeed, in the adaptive multiresolution semi-Lagrangian method designed in [START_REF] Besse | Adaptive multiresolution semi-Lagrangian discontinuous Galerkin methods for the Vlasov equations[END_REF] the authors use interpolating wavelets (interpolets) which are built on polynomial spaces (Lagrange, Legendre, Chebyshev, ...) and enjoy a fast wavelet transform (FWT) like the FFT for the Chebyshev-Fourier representation.

An easier alternative approach to spatial adaptivity is to implement a global refinement method via uniform dyadic grids. This can readily be achieved in our implementation, provided we stay in the full cylinder. More precisely, in exploiting the nesting property of our Chebyshev-Fourier grid of sizes (2 N + 1) × 2 M , one may double the grid size, whenever the spectral dimension is insufficient to accurately represent the flow fields, i.e. when the smallest modes become non-zero. This refinement mechanism is described in detail in [START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF]Section 4.2], but it is yet only implemented in the CLA to improve the interpolation of the flow fields back onto the Eulerian grid. 15 In using a global refinement method with the mentioned dyadic grids, we can approach a potentially singular time of a flow much faster and with less computer resources, because we can start a simulation on a coarse grid and refine it when needed. Additional perspectives for future work can be the extension of the presented CLA to fully three-dimensional non-axisymmetric flow and to flows in other domains with boundaries. In this spirit, it is crucial to simulate other potentially singular flows in bounded domains and check numerically the location of blowing-up gradients in the velocity field. Moreover, in the fully three-dimensional case without symmetries and for potential singularities living on a low-dimensional subspace (typically a point or a one-dimensional curve) an adaptive multiresolution Cauchy-Lagrange method, based on a multi-wavelets representation, should be essential. The presence of boundaries may, after all, be a driving factor for the appearance of finite-time singularities.

Appendices

A. Chebyshev Polynomials and Series Representation

Introductions to Chebyshev polynomials can be found plentiful in the literature (e.g. [START_REF] Mason | Chebyshev polynomials[END_REF], [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF], [START_REF] Gautschi | Orthogonal polynomials[END_REF]). Here, we will only recall from the given references the definitions and properties that are used in our computations.

Definition. The Chebyshev polynomial of the first kind of order n is defined as

The Chebyshev extrema, i.e. the local extrema of T n (x) in the interval I, are given by

x i = cos iπ n , i ∈ {0, . . . n}, (170)
such that T n (x i ) = (-1) i , the zeros are given by xi = cos 2i + 1 2n π , i ∈ {0, . . . , n -1} .

Without difficulty, one verifies the nesting property of the Chebyshev extrema, that is

cos iπ n 0≤i≤n ⊂ cos iπ 2n 0≤i≤2n , ∀n ∈ N, (172)
or, in other words, the extrema of T n are also extrema of T 2n . Furthermore, there holds the orthogonality relation

1 -1 T l (t)T k (t)ω(t) dt = 1 -1 T l (t)T k (t) √ 1 -t 2 dt = c l π 2 δ lk , (173)
where

c l := δ l0 + 1 = 2, if l = 0, 1, if l > 0, (174)
and δ ( . ) denotes the Kronecker symbol. A discrete version of the latter relation, based on the Chebyshev extrema (x i ) 0≤i≤N ∈N from above, reads

N i=0 1 c i c N -i T l (x i )T k (x i ) = c l c N -l 2 N δ lk , (175)
for l, k = 0, . . . , N . In the literature the above summation

with the special case

xT n = 1 2 (T n+1 (x) + T |n-1| (x)). ( 177
)
The derivatives of the Chebyshev polynomials fulfill the relations

T 0 = T 1 , 4T 1 = T 2 , 2T n = T n+1 n + 1 - T n-1 n -1 , n ≥ 2, (178)
as well as

T n = n-2 k=1 k+n even 1 c k n(n 2 -k 2 )T k . (179)
The

From eq. (178), in using eq. ( 180), we directly obtain the derivative relations , n ≥ 3,

2T * 0 = T * 1 , 8T * 1 = T * 2 , 4T * n = T * n+1 n + 1 - T * n-1 n -1 , n ≥ 2, (184)
where the constant of integration has again been omitted. In the following we assume to be given function values u * i = u(x * i ), i = 0, 1, . . . , N , of a continuous function u : [0, 1] → R on the shifted Chebyshev extrema

x * i :=

x i + 1 2 ∈ [0, 1], (187)
with x i as in eq. ( 170). Note that the shifted Chebyshev extrema fulfill the same nesting property (172) as the non-shifted Chebyshev extrema. An approximation u N of u is given by the shifted Chebyshev series

u N (r) = N l=0 u l T * l (r), (188)
such that u N (x * i ) = u * i . The series coefficients u l are given by

u l = 1 c l c N -l 1 N u * 0 + (-1) l u * N + 2 N -1 i=1 u * i cos ilπ N . (189)
The bracketed term is commonly known as the type-I discrete cosine transform, or DCT-I, of the vector (u * i ) 0≤i≤N ; it is implemented in numerous open source software 16 . The transform from function values based on the Chebyshev zeros to the Chebyshev coefficients u l is defined by a DCT-II, whose backtransform is a DCT-III. However, for our numerical needs, namely the inclusion of interval end-points (boundary points of the cylinder) and the nesting property (172) that is not fulfilled by the Chebyshev zeros, we use the shifted Chebyshev extrema throughout this work. However, the zeros may be used for integration purposes in form of the Gauss-integration scheme, which is reviewed here briefly as well. See [START_REF] Canuto | Spectral methods: fundamentals in single domains[END_REF][START_REF] Mercier | Analyse numérique des méthodes spectrales[END_REF] for the following statements and proofs thereof. Gauss integration for Chebyshev Polynomials. Let w 0 , . . . , w n-1 be the solution of the linear system

If a function g = g(r) is given by a shifted Chebyshev series with coefficients {g (0) l } l=0,...,N , then the first-order derivative can be expressed by a shifted Chebyshev series as well, with coefficients {g

(1)

N +1 = g (1)
N = 0 and c l as in (174). Higher order derivatives can be obtained by a re-application of this relation.

Figure 1 :

 1 Figure 1: Convergence and accuracy verification for the methods L -1,(2) 1 and L -1,(1)/r 1

Figure 2 :

 2 Figure 2: Derivative of the approximated solution to eq. (64) divided by r and derivative of the approximated solution to eq. (65) for L = 2π and L = 1/6. The analytic solutions are u(r, z) = w(r, z) = u1(r) e sin(2πz/L) , with u1 as in Figure 1. The relative errors are computed with the l 2 -norm on the fixed Chebyshev-Fourier grid.

( a )

 a Particles move onto the slice that contains the Chebyshev-Fourier grid. (b) Particle displacement on the (r, z)plane on a 5 × 4 grid. The upper row of particles have been added by periodicity.

Figure 3 :

 3 Figure 3: The modelled slice of the cylindrical domain D(1, 2) with the underlying Chebyshev-Fourier grid is shown. Only the radial and vertical movements of the particles are relevant. The top particle positions in (b) are obtained by periodicity of the flow from those that correspond to the first horizontal grid line.

Figure 4 :

 4 Figure 4: Approximation of the limit in eq. (139) in the l ∞ -norm on the underlying gird. The curves overlap to at least s = 30, from which a possible limit at around 0.4 can be inferred. The sudden drop in values indicates that the working precision is insufficient to accurately represent the time-Taylor coefficients.

Figure 5 :

 5 Figure 5: Plot of the function F defined in (142) for different time-steps. The recursion mechanism may be stopped, when the values of F drop below machine epsilon εm, here in double precision and indicated by a dotted line in the plot. All higher order terms X (s) after the criterion has been met give the same summation result in the time-Taylor series in the working precision. A break in the linear behavior of F indicates the depletion of the number range in double precision in the coefficients X (s) .

ωFigure 6 :

 6 Figure 6: Stationary flow components. Left: radial and vertical vorticity for κ = 1 and with c1 (first root).Right: same as left with κ = 5 and c5(fifth root).

Figure 7 :

 7 Figure 7: Relative errors (as defined in (158)) in the computed coefficients X(1) , . . . , X(9) , which belong to the flow in eq. (156) with κ = 1 and c1. The Poisson problems were solved in double precision via the Chebyshev-Galerkin solver (left) and in quadruple precision via the Bessel-Chebyshev solver (right). The true coefficients were obtained via symbolic differentiation inside a computer algebra system by the method described in Section B.

Figure 8 :

 8 Figure 8: A visualization of relative errors in the velocity (left) and vorticity (right) fields due to different resolutions. The fixed run-time parameters are: S = 25, time-step 0.05, 5-th order B-spline cascade interpolation.

Figure 9 :

 9 Figure 9: The effect of different time-steps (left) and orders in the B-Spline cascade interpolation (right).Larger time-steps and lower interpolation orders yield the best results. The run-time parameters are as in Figure8.

Figure 10 :

 10 Figure 10: Test flow and standard conserved quantities, computed on 513 × 1024 points, with S = 40, timestep t = 0.1, Chebyshev-Galerkin solver and 5th order B-spline cascade interpolation.

Figure 11 :

 11 Figure 11: Test flow and conserved quantities that are particular to the CLA; same run-time parameters as in Figure 10.

Figure 12 :

 12 Figure 12: Test results for the given swirl-free flow with a = 1 and n = 2, and run-time parameters: S = 20, t = 0.1, Chebyshev-Galerkin solver, 5th order B-spline cascade interpolation, no dealiasing, on a 513 × 1024 Chebyshev-Fourier grid. Left: kinetic energy conservation and movement of the boundary particles in radial direction, all the conserved quantities, which are not plotted here, stayed identically zero in the simulation; Right: verification of the min/max principle for ω α /r | r =0 , and verification of eq. (160).

• 10 - 3 ω

 3 ∞ 3.76999 • 10 3 9.0847 • 10 4 1.54276901 • 10 5 4.3127 • 10 6 5.8413 • 10 9 1.2401 • 10 12

Figure 13 :

 13 Figure 13: Resolution study with fixed parameters: S = 15, t = 10 -5 , Bessel-Chebyshev solver, 5th order B-spline interpolation method. The legend in the vorticity growth plots applies to all other panels.Here and in the following plots, the black crosses in the vorticity plot are the exact coordinates from Table1at times t = 0.003, t = 0.0031 and t = 0.0034. Also the helicity panel displays absolute errors, because the initial helicity vanishes for this flow.

Figure 14 :

 14 Figure 14: Several fixed time-steps are compared on a 513 × 1024 grid; same run-time parameters, other than time-step, as in Figure 13. Larger time-steps deliver good and fast results. The legend in the first panel is valid for the others too. Note that the fourth panel does not show the whole integration period.

Figure 15 :

 15 Figure 15: The effect of different truncation orders in the adaptive time-stepping. Higher orders imply larger time-steps and fewer iterations. These runs were executed on a 513 × 3072 grid, the other parameters are the same as in the plots from above. Sudden drops in the size of the time-steps in the lower left panel at t = 0.003 (and selected later times) are forced time-steps to arrive exactly at that time for a comparison to the values in Table1.

Figure 16 :

 16 Figure 16: Energy spectrum at the boundary r = 1 for various time instances. Small modes fill up rapidly and the flow can soon not be resolved accurately anymore. The underlying grid size was fixed to 513 × 4096 Chebyshev-Fourier points. All other fixed run time parameters are given in the text.

Figure 17 :

 17 Figure 17: Higher resolution runs are shown. Here, the truncation order was set to S = 18 to maintain the initial time-step 10 -5 longer. The interpolation was done by the 5th order B-spline cascade; no dealiasing is used here. Sudden breaks in the conserved quantities arrive much later than in the previous plots, indicating a lasting stability with increased resolution.

Figure 18 :

 18 Figure 18: Upper row: Plot of the velocity components at r = 1 until t = 3.51 × 10 -3 . The values were obtained from a simulation with 513 × 4096 grid points, the other run time parameters were the same as in the preceding figures. The velocity component v r in the top left panel should stay zero at r = 1 due to the no-flow boundary condition, but large errors appear eventually at the symmetry axes. Lower row: Plot of the vorticity components at r = 1.Here, the evolution of the components is only displayed until t = 3.0 × 10 -3 , the growth of the local extrema is too strong to be conveniently visualized. For surface plots of the flow field components, we would like to refer the reader to[START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF].

T

  n (x) = cos(n arccos(x)), x ∈ I, n ∈ N 0 ,(167)whereI := [-1, 1] and N 0 := N ∪ {0} = {0, 1, 2, . . .}. Setting x = cos(ϑ), ϑ ∈ [0, π], we find T n (cos(ϑ)) = cos(nϑ),(168)from which all the properties in the following may be derived. The Chebyshev polynomials satisfy the three term recurrence relationT 0 (x) = 1,T n+1 (x) = 2xT n (x) -T |n-1| (x), ∀n ∈ N 0 .

N i=0 1 /

 1 (c i c N -i ), indicating the halving of the first and last terms, is often denoted by and N i=0 1/c i by . Here we intend to avoid primed sums, the appearing primes in the text represent derivatives with respect to the sole argument of a function. A multiplication of two Chebyshev polynomials entails 2T n T m = T n+m + T |n-m| , ∀n, m ∈ N 0 ,

4 j

 4 shifted Chebyshev polynomial of the first kind of order n ∈ N 0 is obtained from the ordinary Chebyshev polynomial by a change of the variable, namely t = 2r -1 with r ∈ [0, 1], and is denoted byT * n (r) := T n (2r -1), r ∈ [0, 1] . (180)The shifted Chebyshev polynomials satisfy the three term recurrence relationT * 0 (r) = 1, T * n (r) = 2(2r -1)T * n-1 -T * |n-2| , n ∈ N, + 2T * n + T * n+1 , ∀n ∈ N 0 .(182)An explicit representation of the shifted Chebyshev polynomial of order n > 0 is given by the polynomial sum (-1) n-j (j + n -1)! (2j)!(n -j)! r j .

  up to an additive constant of integration. The next integral relation follows readily from eqs. (182) and (185), τ (n, r) := rT * n (

n- 1 j=0(x j ) k w j = 1 - 1 x- 1 j=0p(x j )w j = 1 - 1 p

 111111 k ω(x) dx, 0 ≤ k ≤ n -1,where x0 < x1 . . . < xn-1 are the roots of the n-th Chebyshev polynomial from eq. (171). Then there holds w j = π n , for j = 0, . . . , n -1, and n(x)w(x) dx, for all p ∈ P 2n-1 .

( 1 )

 1 l } l=0,...,N -1 . The latter coefficients are obtained via the relation c l-1 g

l

  , l = N, N -1, . . . , 1,

Table 1 :

 1 The approximate maximal vorticity values in this table are taken from

Table 2 :

 2 The body of this table states the computed maximal vorticity at selected times and grid sizes. The values in the last row are those of Table1and come from the given references. 10 4 1.54167635 • 10 5 1.60528443 • 10 6 513 × 4096 9.08466826 • 10 4 1.54277249 • 10 5 3.13597705 • 10 6 1025 × 8192 9.08466912 • 10 4 1.54277317 • 10 5 3.82533611 • 10 6 1025 × 10240 9.08466912 • 10 4 1.54277318 • 10 5 3.94749487 • 10 6Our maximal vorticity values show good agreement with those from the references for t = 3 • 10 -3 and t = 3.1 • 10 -3 . The value at 3.4 • 10 -3 is approached but only agrees in magnitude. Those are exceptional results, in particular, if the execution time of the simulations is considered. All of the computations presented in Table2and fig.17were carried out on a single CPU12 and finished within one day. Typically, implementations of CFL bound algorithms need multiples of that time span on a comparable CPU and execute thousands of time-steps. The run of the CLA on the highest resolution in Table2reached t = 3.1•10 -3 in only 310 iterations and t = 3.4•10 -3 in 344.

	time	3.0 • 10 -3	3.1 • 10 -3	3.4 • 10 -3
	513 × 1024 9.08451276 • ref. [45, 1] 9.0847 • 10 4	1.54276901 • 10 5 4.3127 • 10 6

It is explicitly given in eqs. (5) and (6) in beginning of section

[START_REF] Bateman | Tables of integral transforms[END_REF] See Section 3 for a precise definition of the modeled cylindrical domain.

as far as the computational limits, i.e. working precision, storage, etc., allow.

Other sets of basis functions are possible, e.g. B 2 := {(1 -t 2 )T l (t)} l=0,...,N -2 , (see[START_REF] Hertel | The Cauchy-Lagrange algorithm for axisymmetric incrompressible Euler flow in a wall bounded cylindrical domain[END_REF] for more details on this choice).

e.g. https://people.sc.fsu.edu/~jburkardt/f_src/lagrange_interp_1d/lagrange_interp_1d.html

The author of this module is Jacob Williams, it can be obtained from https://github.com/jacobwilliams/ bspline-fortran/blob/master/src/bspline_sub_module.f90.

In the sequel of this article we shall refer to this quantity as the total energy, or simply energy.

i.e. when those conserved quantities show an error above 10 -5 .

We use |H(0) -H(t)| instead of |H(0) -H(t)|/|H(0)|, since the initial helicity H(0) is zero for this flow.

Here, we mean concerning a larger integration period together with improved conservation of the discussed quantities.

[START_REF] Canuto | Spectral methods: fundamentals in single domains[END_REF] mainly the 2/3-rule and double length methods

The cluster in use provides Intel Xeon Gold 6148 CPUs with 20 cores and 40 possible threads, which are individually supplemented by 192GB of memory (RAM).

Here, we refer to double (8 Byte floats), extended double (10 Byte floats), or quadruple (16 Byte floats) precision computations

The appearing linear systems in the CG method are internally solved with the LAPACK package, which is restricted to double precision computations.

For the results in this work, we have not applied the mentioned refinement methods in the interpolations.

We use FFTW 3.3.9, see http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html for a definition of their implemented transforms.

Acknowledgements

We are very grateful to Claude Bardos, Sai Swetha Venkata Kolluru and Rahul Pandit for fruitful discussions. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Also the division g(r)/r can be given by a shifted Chebyshev series with coefficients {a l } l=0,...,N , in solving the linear system of equations

where

B. Explicit Computation of the Time-Taylor Coefficients for Stationary Solutions

We depart from the characteristic equation and leave the dependence on the arguments out for convenience (the following is true in any coordinate system),

for a classical stationary solution V = V (a) to the three dimensional incompressible axisymmetric Euler equations. If X(t) is analytic in time, then so is the velocity along the trajectories V (X(t)) due to eq. ( 194). Let g(t) := V (X(t)), then

One finds

and after collecting powers of t, X (s) = g (s-1) , s ≥ 1 .

The first coefficients are, thus, given by

.. Despite the fact that the general form may be given by a recursion formula and the formula of Faa di Bruno [START_REF] Bruno | Note sur une nouvelle formule de calcul différentiel[END_REF][START_REF] Gzyl | Multidimensional extension of Faa di Bruno's formula[END_REF], here, we are mostly interested in the explicit calculation of the time Taylor coefficients in cylindrical coordinates. This is possible by hand in the case of the vortex, but rather exhausting in the case of more intricate solutions as presented in [START_REF] Leprovost | Dynamics and thermodynamics of axisymmetric flows: Theory[END_REF]. We will therefore make use of computer algebra systems, but first state the recursion in cylindrical coordinates.