{'fraction_non_alphanumeric': 0.04053036190385098, 'fraction_numerical': 0.010927655097380471, 'mean_word_length': 4.440837832563585, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Cauchy-Lagrangian method for numerical analysis of Euler flow

A novel semi-Lagrangian method is introduced to solve numerically the Euler equation for ideal incompressible flow in arbitrary space dimension. It exploits the time-analyticity of fluid particle trajectories and requires, in principle, only limited spatial smoothness of the initial data. Efficient generation of high-order time-Taylor coefficients is made possible by a recurrence relation that follows from the Cauchy invariants formulation of the Euler equation (Zheligovsky&Frisch, J. Fluid Mech. 2014, 749, 404-430). Truncated time-Taylor series of very high order allow the use of time steps vastly exceeding the Courant-Friedrichs-Lewy limit, without compromising the accuracy of the solution. Tests performed on the two-dimensional Euler equation indicate that the Cauchy-Lagrangian method is more - and occasionally much more - efficient and less prone to instability than Eulerian Runge-Kutta methods, and less prone to rapid growth of rounding errors than the high-order Eulerian time-Taylor algorithm. We also develop tools of analysis adapted to the Cauchy-Lagrangian method, such as the monitoring of the radius of convergence of the time-Taylor series. Certain other fluid equations can be handled similarly.


Introduction
As is well known, fluid flow can be characterised in terms of the current positions of fluid particles (Eulerian coordinates), or in terms of their initial positions (Lagrangian coordinates). For ideal (inviscid) incompressible fluid both formulations were introduced in the 18th century [15,28]. Lagrange himself was not particularly fond of the "Lagrangian" form of his equation (obtained from a variational formulation) and quickly shifted to what is now known as the Eulerian form. From a numerical point of view, the Eulerian formulation seems to have a significant edge because it looks much simpler; the Lagrangian formulation has even been qualified "an agony", because of its complexity [36].
However, the Eulerian time-stepping methods have one serious well-known drawback, the Courant-Friedrichs-Lewy (CFL) condition [12], which constrains the time step to be less than a dimensionless constant multiplied by the time needed to sweep across the spatial mesh at the maximum flow speed. As a consequence, the complexity of computations with N collocation points in each spatial direction is roughly O(N 4 ) in three dimensions. Hence, progress in high-resolution numerical simulation is creepingly slow, in spite of Moore's law. (As we will see in Section 5, this is not the only drawback of Eulerian schemes.) Foremost because of the CFL constraint, which may conflict, for example, with the desire to quickly produce well-resolved numerical weather forecast, there has been a strong incentive to develop semi-Lagrangian (SL) schemes, in which some form of Lagrangian integration -with time steps not constrained by CFL -alternates with a reversion to a Eulerian grid (see, e.g., [40] and references therein). The SL algorithms used so far, e.g., in geophysical fluid dynamics, engineering, mechanics and plasma physics, were designed for situations where satisfactory results can be obtained with rather low-order temporal schemes. They are thus not appropriate for numerical experimentation on delicate questions, such as the issue of blow up in three-dimensional flow (see, e.g., [22]).
A mostly forgotten Lagrangian formulation of the equations of ideal incompressible flow, due to Cauchy [9,20], combined with some recent advances regarding the time-analyticity of Lagrangian trajectories [21,44], allows us to propose a new SL algorithm, called the Cauchy-Lagrangian algorithm, that is particularly well-suited for problems in which high precision is a prerequisite, and actually better suited than Eulerian schemes.
Cauchy's 1815 Lagrangian equations are formulated in terms of the Lagrangian map from Lagrangian to Eulerian coordinates, a → x = x(a, t), defined as the solution of the ordinary differential equation for fluid particle trajectory,ẋ = v(x, t), with the initial condition x(a, 0) = a. Here, the dot denotes a (Lagrangian) time derivative. The three-dimensional Cauchy invariants equations are (see Section 2 for a simplified derivation): where ω (init) = ∇ L × v (init) denotes the initial vorticity and ∇ L is the Lagrangian gradient, i.e., the space derivatives in a. Since the r.h.s. of the first equation in (1) does not depend on time, its l.h.s. is obviously a Lagrangian invariant, whose scalar components are called the "Cauchy invariants". These invariants were much later interpreted as a consequence of Noether's theorem applied to a continuous symmetry of the Euler equation, called the relabelling invariance (see [20], Section 5.2). It was shown in [21,44] that the Cauchy invariants equations (1) imply analyticity of fluid trajectories for initial conditions that have certain rather weak regularity. The proof of analyticity is derived from explicit recurrence relations for coefficients of the time-Taylor series for the Lagrangian displacement ξ = x − a. These can also be used to construct a numerical Lagrangian method of very high order in time in both two dimensions (2D) and three dimensions (3D). Its time steps are only constrained by the radius of convergence of the Taylor series (typically, the inverse of the largest initial velocity gradient). For more details, see Section 2 and [21,44].
The paper is organised as indicated hereafter. In Section 2 we recall some of the known results about Cauchy's Lagrangian formulation and its application to the time-analyticity of the Lagrangian map. In Section 3 we describe the Cauchy-Lagrangian (CL) algorithm in detail: we begin with an overview and show that CL may be considered as a semi-Lagrangian algorithm of arbitrary high order (Section 3.1), present the interpolation technique for reverting to Eulerian coordinates at the end of each time step (Section 3.4) and show how to optimise the choices of the time step and of the order of the Taylor expansion truncation to minimise computational complexity (Section 3.5). Section 4 is devoted to testing the CL algorithm in the 2D case against various numerical methods, and to CPU benchmarks. Section 5 is a comparison of high-order time-Taylor expansions in Eulerian and Lagrangian coordinates: the Lagrangian method suffers much less from the rounding errors. In Section 6 we determine how quickly we have to decrease the time step in the CL method due to the fact, that radius of convergence R(t) of the time-Taylor series around time t generally shrinks as t increases, as more and more small-scale eddies (high spatial Fourier harmonics) are generated. Finally, in Section 7 we present concluding remarks and point out that the Cauchy-Lagrangian method is well adapted for certain other problems for ideal fluid flow.

Mathematical background
Here we recall some of the background material which is used to develop the Cauchy-Lagrangian method: derivation of the Cauchy invariants equations and recurrence relations for time-Taylor coefficients from which follows the analyticity of the time-Taylor series (Section 2.1). We then present, in a heuristic way, a result allowing the determination of the radius of convergence of such Taylor series, an important new tool for analyzing the output of Cauchy-Lagrangian computations (Section 2.2).

Cauchy's Lagrangian formalism for ideal incompressible fluid flow
The flow of an ideal incompressible fluid, described in Eulerian coordinates, is governed by the Euler equation where v(x, t) is the velocity and p(x, t) the pressure (divided by the density, which in the incompressible case is constant). Here the spatial differentiation ∇ is performed in the Eulerian coordinates x. The Eulerian equation (2) is notorious for difficulties in its investigation, both analytical and numerical. It turns out that its Lagrangian analogue can be once integrated in time, which yields the Cauchy invariants equations (1). Since these equations are central to our numerical method, we briefly recall how they are derived.
In Lagrangian coordinates, fluid particles are characterised by their initial positions a and the time t. The subsequent positions of fluid particles are then described by the Lagrangian map a → x(a, t). The spatial gradient in the Lagrangian coordinate a is denoted ∇ L .
Cauchy's derivation of (1) is now briefly recalled [9]. Observe that the l.h.s. of (2) is the accelerationẍ of the fluid particle; thus,ẍ = −∇p. Multiplying this equation by the transpose of the Jacobian matrix ∇ L x, we transform the Eulerian gradient in the r.h.s. into the Lagrangian one and obtain, in the component form: Cauchy then applies a Lagrangian curl to (3), and finds that the l.h.s. can be exactly integrated in time to yield the first equation of the set (1). The second one, the Jacobian equation, expresses the conservation of volume, equivalent to incompressibility. The Cauchy invariants equations are quadratically nonlinear equations, not resolved in the time derivatives; the Jacobian equation is also quadratic (in 2D) or cubic (in 3D). At first glance, it is unclear, in what respect (1) are advantageous compared to the Eulerian equations (2). To see what is gained, following [21,44], we consider the equations for the displacement ξ = x − a, that ensue from (1) and are in 3D as follows: and expand the displacement in the time-Taylor series The structure of equations (4)-(5) enables us to derive recurrence relations for the coefficients ξ (s) (a), which in 3D take the following form [21,44]: Here δ s 1 is the Kronecker delta and ε ijk is the unit antisymmetric tensor. To determine ξ (s) , knowing its (Lagrangian) curl and divergence, is a Helmholtz-Hodge problem, with a unique solution for suitable boundary conditions. For example, for space-periodic flow we demand that the averages of all ξ (s) over the periodic box vanish, and construct the solution of (7)-(8) by taking the gradient of (8), subtracting from it the curl of (7) and thus obtaining a Poisson equation for ξ (s) . In particular, for s = 1, (7)-(8) are equivalent to the relation ξ (1) The explicit character of the recurrence relations (7)-(8) permits to easily derive bounds for the time-Taylor coefficients, provided the initial condition possesses a certain minimum smoothness, e.g., if the initial vorticity is Hölder-continuous (i.e., belongs to the space C α for some 0 < α < 1) or belongs to the space of absolutely summable Fourier series (henceforth, ASFS). This way one can prove various analyticity theorems [21,44], which state that the Lagrangian gradient of the Taylor series (6) converges in the space to which the initial vorticity belongs, provided |t| < R bound . Here, R bound denotes a lower bound for the radius of analyticity, which is typically order of the inverse C α or ASFS norm of the initial vorticity. Within the disk of radius R bound , convergence of the time-Taylor series (6) and of its spatial gradient are guaranteed and so is analyticity in time. The interested reader will find elementary proofs in [44]. Of course, this result is key for the construction of numerical solutions. Note that in 2D, by a theorem of Wolibner [42], when the initial vorticity is Hölder-continuous, the vorticity stays smooth for all positive times. In 3D this problem is open.
Suppose that, by using the Taylor series around t = 0 we have solved (1) in some real time interval [0, t 1 ], where t 1 < R bound . If the solution at time t 1 has the required smoothness, it can be used as a new initial condition for (1) at time t 1 . For this, Lagrangian coordinates must be reintroduced for fluid particles starting at t 1 . At t = t 1 the new Lagrangian coordinates coincide, of course, with Eulerian coordinates. Hence we need to revert to Eulerian coordinates at t 1 . In principle, this can be done by composing the Lagrangian fields in the original coordinates with the inverse of the Lagrangian map, but there are other ways, as we will see. This procedure can be extended into a Weierstrass-type analytic continuation: one repeats the process at times t 1 < t 2 . . ., as long as the successively constructed solutions have enough spatial smoothness to give time-Taylor series with a finite non-vanishing radius of convergence. In 2D, when the initial vorticity is Hölder-continuous, the required smoothness will persist forever. In 3D loss of smoothness (blow up) cannot be ruled out.

The radius of convergence of the time-Taylor series
When the initial vorticity is bounded in a suitable function space, such as C α or ASFS, we know that, for |t| < R bound , the gradient of the Taylor series (6) is guaranteed to converge in that space for all a. Convergence may hold over a larger time interval, but without any proven control over the spatial smoothness of the Lagrangian map thus constructed. Knowing the radius of convergence of this Taylor series is important for numerical applications, even if |t| does not approach the edge of the disk of convergence: this parameter controls how fast the terms in the series tend to zero at large orders s.
For a given starting point a, there exists a non-negative number R(a, 0), called the radius of convergence, such that the series (6) is absolutely convergent for any 0 ≤ |t| < R(a, 0). Here, |t| denotes the modulus of t, which is allowed to take complex values. The second argument of the radius R indicates the time around which the Taylor expansion is performed. Note that the radius of convergence can be infinite. It is given by the Cauchy-Hadamard formula (see, e.g., [29] It might be natural to define a radius of convergence R(0) by taking the infimum of R(a, 0) over all initial points a. How can we measure this infimum in practice, given a set of Taylor coefficients ξ (s) (a) up to a large value of s, and assuming that the a values are on a large set of collocation points? Following the usual numerical analyst's approach, one might try to just evaluate the local radii of convergence R(a, 0) on a grid in the fluid domain and then find the minimum. This approach would have the advantage that it might indicate, what region in space is characterised by small values of the local radii of convergence, eventually perhaps leading to blow up. However, the techniques for estimating the radius of convergence are numerically quite demanding, as discussed in Appendix B, suggesting that this approach is not numerically optimal. Furthermore, when the radius of convergence R(a, 0) is discontinuous, such an approach is clearly inappropriate.
Even if ξ (s) (a) are analytic in a, we cannot guarantee the continuity of R(a, 0) (the series s≥0 2 −|a−a⋆| 2 s 2 t s , unrelated to the Euler equation, illustrates this difficulty). If R(a, 0) takes small values at exceptional off-grid points, we will miss them when measuring R(0) numerically on a discrete grid. We find it therefore more suitable to employ the quantity where "ess" stands for "essential", namely, small outlier values on a zero-measure set are disregarded.
A simple way to evaluate R inf (0) makes use of the ordinary power series obtained by replacing ξ (s) (a) by its L p norm, namely, We denote by R p the radius of convergence of the series (11). We assume that the initial vorticity is space-periodic and is either Hölder-continuous or has an absolutely summable Fourier series. The same properties hold then for all the gradients of the Taylor coefficients ∇ L ξ (s) (a). As a consequence, all the Taylor coefficients ξ (s) (a) are in the Lebesgue spaces L p of functions, such that the pth power of their absolute value is integrable. In Appendix A we prove a theorem stating, roughly, that the radius R inf (0) is also the radius of convergence of the series (11). More precisely, 1. For any p ≥ 1, we have R inf (0) ≥ R p . 2. If the displacement ξ(a, t), defined by (6), is in L p for any 0 ≤ t < R inf (0), then R inf (0) = R p .
(The proof of this theorem does not involve the Euler equation.) For certain solutions of the Euler equation, the assumption in the second part of the theorem is known to be satisfied. For example, for 2D flow whose initial vorticity is Hölder-continuous, Wolibner's theorem [42] implies that the Lagrangian displacement is always in L p . In general, for 3D Euler flow, we do not know the spatial regularity properties of the displacement beyond the time of guaranteed analyticity and thus we only know that R inf (0) ≥ R p .
To gain heuristic insight into why the theorem holds, we observe that since the radius of analyticity depends only on the modulus of the Taylor coefficients, no generality is lost by taking them to be scalar and non-negative functions of a. For such a Taylor series, with radius of convergence R(a, 0), the coefficient of order s is roughly proportional to R(a, 0) −s (this would be exact for a geometric series). To calculate the L p norm of the coefficient of order s, we must integrate the pth power of this Taylor coefficient over a and take the pth root. It seems plausible that, for large s, these integrals will be dominated by those a where R(a, 0) is close to its infimum. We thus expect that the L p norm of the Taylor coefficient of order s is roughly proportional to (ess inf a R(a, 0)) −s , thereby yielding a radius of convergence given by R inf , namely, (10).

Description of the Cauchy-Lagrangian algorithm
In Section 3.1 we show that the basic ideas of the Cauchy-Lagrangian (CL) algorithm are quite simple, given our understanding of the analytical properties of the Lagrangian map, which were recalled in Section 2.1. Then, in Sections 3.2-3.5 we successively examine the constraints on the initial conditions (in connection with the spatial truncation errors), the temporal truncation errors, the reversion from the Lagrangian to the Eulerian representations and the minimisation of computational complexity.

Overview: CL as a high-order semi-Lagrangian method
Our purpose is to solve the Euler equation (2) with the initial condition v( In principle, we can just require that v (init) or ω (init) be in a function space that guarantees finite-time analyticity of the Lagrangian map, such as (for the vorticity) C α or the space of absolutely summable Fourier series ASFS. For practical reasons, more regularity is desirable (see Section 3.2).
The simplest case is when T is less than the radius of analyticity bound R bound (0) of the time-Taylor series of the Lagrangian map, starting at t = 0. We then use the time-Taylor expansion of the displacement (6), together with the recurrence relations (7)-(8) to calculate as many Taylor coefficients as needed (see Sections 3.3 and 3.5). Here we just observe that, for spatial periodicity conditions, it is straightforward to calculate the Taylor coefficients using a pseudospectral method with dealiasing. In 3D, we can then obtain the vorticity in Lagrangian variables at time T by applying Cauchy's vorticity formula [9]: (In 2D, this formula reduces to the constancy of the vorticity ω L (a, T ) = ω (init) (a).) Finally we can revert to Eulerian coordinates by composing this Lagrangian vorticity with the inverse Lagrangian map a(T ) : x → a(x, T ), which maps the fluid particle position x at time T to its initial position a, to obtain Numerically, this is best done by interpolating from a Lagrangian to a Eulerian grid (see Section 3.4). Finally, the Eulerian velocity at time T , if needed, can be obtained either from the Eulerian vorticity at time T (by inverting the curl, e.g., applying Fourier methods, or by interpolating to a Eulerian grid the time derivative of the Lagrangian map, obtained by time differentiation of the Taylor expansion (6)). Now, consider the case T ≥ R bound (0) when convergence of the Taylor series is not guaranteed. Lack of convergence does not imply loss of regularity at some real time t between 0 and T . It may just be due to complex-time singularities whose moduli are less than T . Even if the Taylor series still converges, the sum might be insufficiently smooth to allow starting a new Taylor expansion (see end of Section 2.1). This is impossible in 2D with a Hölder-continuous initial vorticity, but cannot be ruled out in 3D. The way to handle times T that are beyond the radius of analyticity bound R bound (0), is the multistep method, which essentially implements numerically the analytic continuation explained in Section 2.1. We pick a time t 1 such that 0 < t 1 < R bound (0) and start the process all over at time t 1 using the Eulerian field v(x, t 1 ) as the initial velocity (or ω(x, t 1 ) as the initial vorticity). We then have, in the complex time plane, a new disk of guaranteed analyticity centred at t 1 with radius R bound (t 1 ), which may well extend on the real positive time axis beyond the furthest point reached by the first disk of guaranteed analyticity. The process may be continued (see Fig. 1) with analyticity disks centred at various times 0 < t 1 < t 2 < . . ., as long as the solution remains in a space such as C α or ASFS which guarantees that R bound > 0.
The algorithm outlined above is illustrated by the flow chart in Fig. 2.  To finish this overview, we observe that our Cauchy-Lagrangian method belongs to the family of semi-Lagrangian (SL) methods, often employed, e.g., in meteorology, engineering, mechanics and plasma physics [5,40]. An SL method solves an evolution equation along trajectories of particles for one time step and then reverts back to the static regular Eulerian grid, usually by interpolation. Unlike pure Eulerian ones, SL methods are not subject to the CFL constraint [12]; consequently, their time steps are independent of the spatial resolution. For time-analytic Lagrangian fluid particle trajectories, the time step must just be smaller than the radius of convergence of the time-Taylor series, which is typically the inverse of the maximum velocity gradient.
SL methods fall into two categories: upstream methods compute the so-called departure points, i.e., the positions of Lagrangian particles which are advected onto the regular grid during the time step, whereas for downstream SL methods the departure points of trajectories are on the regular grid, and the arrival points are on the distorted grid. The latter methods are simpler, because they do not require solving auxiliary problems to find departure points [30]. In this terminology, our CL method is downstream. Standard SL methods are typically of second-order accuracy in time; numerical schemes of higher orders are not widely used, being too complicated. Indeed, to advance by one time step, high-order SL method of the traditional type would need information about the solution at many time levels; this would require many additional interpolations. In contrast, in the Cauchy-Lagrangian approach, the recurrence relations (7)-(8) allow us to construct very high-order schemes requiring only one interpolation per time step.

Spatial truncation errors and the choice of initial conditions
As we will see in Sections 4 and 5, our CL method is particularly well adapted and efficient for problems where high accuracy is important. Thus we want to keep all errors as (consistently) small as possible. This includes of course spatial truncation errors, stemming from the use of a finite number of modes. Actually, this is much easier if we limit ourselves to initial conditions that are analytic functions of the space variable. Then, the Eulerian solution remains analytic in space for at least a finite time (all times in 2D) [3]. During that time interval, the spatial Fourier transform of the solution falls off roughly as e −δ(t)k , where k is the wave number and δ(t) is the radius of the tube of analyticity, the distance between the real spatial domain and the nearest complex-space singularity for the complexified spatial variables. Hence, truncation errors fall off then as e −δ(t)kmax , where k max is the maximum wave number present in the numerical truncation; thus they remain small as long as δ(t) is significantly larger than the spatial numerical mesh [41].
In contrast, consider what happens when an initial condition has merely the minimum regularity required for ensuring the time-analyticity of the Lagrangian map, e.g., an initial vorticity is Höldercontinuous. The spatial Fourier transforms of the velocity field during the time interval [0, T ] will then fall off in a roughly algebraic way with the wave number k. Thus spatial truncation errors will fall off only algebraically with k max . Achieving high accuracy may then require spatial resolutions that are beyond the reach of any existing computer. To avoid this, in our test computations we limit ourselves to spatially analytic initial data.
One of the simplest space-analytic initial conditions that can be used for testing a CL code is the AB flow, discussed in Section 2.6.1 of [44], a steady-state solution of the 2D Euler equations. Its vorticity is ω (init) = sin x cos y. Although its Eulerian-coordinates dynamics is trivial, its Lagrangian dynamics is not completely trivial, but can still be determined analytically using elliptic integrals. More examples of flows with analytic initial data and truly non-trivial dynamics will be given in Section 4.1.

Temporal truncation errors
In numerical applications the Lagrangian displacement is approximated by the time-Taylor series (6) truncated to some finite order S: Which constraints on ∆t and on S guarantee that the remainder of the series is sufficiently small? First, we demand that the series (6) be convergent, i.e., ∆t should be less than the radius of convergence for all starting points a, and hence less than R inf , the infimum over a of these radii (ignoring sets of zero measure).
The numerical determination of R inf proceeds as follows. Using the material of Section 2.2 and of Appendix A, we replace the full functional Taylor series by the series of the L p norms of the Taylor coefficients. This gives a series with radius of convergence R p . In 2D, R p = R inf and, in principle, in 3D we just have R p ≤ R inf -however, since in 3D a blow up of solutions to the Euler equation in the absence of solid boundaries perhaps never takes place, it is conceivable, at least for space-analytic initial conditions, that R p = R inf . The most convenient is to use the L 2 norm, which can be evaluated using Parseval's theorem. We thus have to determine the radius of convergence R 2 of a power series with positive coefficients. How to do this in practice is discussed in Appendix B.
Second, we demand that the error due to truncation of the time-Taylor series be small. To see what this implies, we consider the remainder of the Taylor series and majorise its L 2 norm as follows: where · denotes the L 2 norm. The series in the r.h.s. has positive coefficients and its sum has a singularity at ∆t = R 2 . Thus the sum is typically well approximated (at large orders and up to logarithmic factors) by a geometric series of ratio r ≡ ∆t/R 2 , and the L 2 norm of the remainder will typically be less than ξ (S) (∆t) S r/(1 − r). For the values of the ratio r that minimise the complexity of the computation (see Section 3.5), r/(1 − r) ≈ 0.15 . Hence, our criterion for keeping the truncation error small is That is, we demand that, in the time-Taylor series, the last term kept before truncation be less, in L 2 norm, than a prescribed accuracy ε. Note that taking larger values of S allows to advance forward with larger time steps. The choice of the optimal truncation order is discussed in Section 3.5.

Reversion to a Eulerian grid by interpolation
In our implementation of the Cauchy-Lagrangian method, the calculation of all the time-Taylor coefficients is done by pseudospectral techniques. For this to be done efficiently, the data at the beginning of a new time step must be known on a regular grid of N d collocation points, where N is the number of grid points per spatial period and d the dimension of space. This will be called here the uniform grid. Note that after the Lagrangian time-stepping is done, the new field, say, the vorticity, can either be determined in Lagrangian coordinates on the uniform grid or in Eulerian coordinates on a distorted grid, which is the image of the original grid by the Lagrangian map during the time interval ∆t, shown in black in Fig. 3. Henceforth, we write just "image" for "image by the Lagrangian map". To repeat the time-stepping process, we need to know the new Eulerian field, say, the vorticity, on a uniform Eulerian grid. We thus need to interpolate the Eulerian field from the distorted grid to the uniform grid. As we will see, the distortion of the grid remains moderate, so that the interpolation can be done using the cascade interpolation introduced in [37], which, in dimension d, is a superposition of d one-dimensional interpolations. We now briefly describe it.
We begin with the 2D interpolation procedure, assuming 2πperiodicity and a time step from t = 0 to t = ∆t. Positions of fluid particles at times 0 and ∆t are a = (a, b) and x = (x, y) in the Lagrangian and Eulerian coordinates, respectively. The Lagrangian map takes a to x(a) = (x(a, b), y(a, b)). The uniform grid is composed of the points a ij = (a i , b j ). Here, a i = 2πi/N and b j = 2πj/N, where the indices i and j run from 1 to the number of grid points per spatial period N. The same uniform grid is used for the Eulerian description after one time step. The N lines (a, b j ), where a varies continuously in the interval [0, 2π], are called the horizontal grid lines. Similarly, the N lines (a i , b) are called the vertical grid lines.
The interpolation employs an intermediate hybrid grid, formed of intersections of the horizontal grid lines with the images of the vertical grid lines. This hybrid grid is shown by stars in Fig. 3. For the cascade interpolation, it is essential that each line of the former and latter type have a single intersection. This follows from a monotonicity result, discussed below.
First, we determine the x-coordinates of hybrid grid points. For each given pair of a horizontal and a vertical line this is done as follows. The images of the discrete grid points on the vertical line have coordinates (x n , y n ). Along the image of the vertical line, in the neighbourhood of the intersection, a 1D interpolation of x as a function of y is done; this gives the x value where y takes the value y j corresponding to the given horizontal line. Second, the vorticity known at the same discrete locations (x n , y n ) along the image of the given vertical line, is also 1D-interpolated in y to the same values y j . Finally, along each horizontal line, the vorticity is 1D-interpolated from the points on the hybrid grid to the points on the uniform grid. It is straightforward to extend the procedure to 3D with grid lines and grid planes and again only 1D interpolations (see [37]). Now we turn to the issue of monotonicity of the map b → y(a i , b). First a qualitative observation: given the time-analyticity of the Lagrangian map and the fact that it is slightly smoother in space than C 1 , for small ∆t the map is close to the identity map and its Jacobian matrix is close to the identity matrix. Hence, ∂y(a, b)/∂b > 0 for sufficiently small ∆t, which implies the required monotonicity. This observation can be made quantitative by using the tools of [44] (Section 2.3). It is shown there (in 3D) that analyticity holds if Γ∆t < T c = (8 − 5 √ 2)/3, where Γ is the sum of the moduli of the Fourier coefficients of the initial vorticity. It is easy to show that the required positivity then holds at least to the same time. The positivity is unlikely to hold if ∆t is just restricted to be less than the radius of convergence R inf of the Taylor series, but in practice, we use values of ∆t significantly smaller (see Section 3.5); in none of the simulations reported in Section 4 have we detected any lack of monotonicity.
In the current version of the code, the 1D interpolations are polynomial and use eight points. The interpolation procedure can be modified in a number of ways. We can interpolate first in the x-variable and then in the y-variable, and compare the results to estimate the accuracy of the interpolation. Instead of polynomial interpolation, we can use rational functions, splines, or any other 1D interpolation algorithm. We can use different numbers of interpolation nodes in different regions of the periodicity square, depending, e.g., on how large is the vorticity gradient (note that in Fig. 7 of Section 4 the areas of large gradients are small for t ≥ 3). For the interpolation we can use a mesh finer, than the one used to compute the Lagrangian map. Spatial derivatives of the vorticity can also be used for interpolation. We will experiment with these modifications in future studies.
Errors of interpolations can be estimated on the fly as follows: the interpolation yields the vorticity on the target regular grid; from this, by FFT we can calculate its Fourier coefficients. By direct summation of the Fourier series ("slow Fourier transform"), we can then compute the vorticity on the distorted grid and find how closely it matches the one we employed for the interpolation. Such computations are CPU-intensive, but it suffices to perform them at those points where the largest interpolation errors are suspected, for example, e.g., where high-order spatial derivatives of the vorticity are large.

Choosing the optimal truncation order and time step
In computations, the displacement is given by the time-Taylor series (14) truncated at order S. For a given accuracy ε, how should we choose S and the time step ∆t to minimise the complexity of computations when integrating in a time interval [0, T ]? We assume that the spatial resolution and thus the complexity of the spatial FFTs are fixed. The order of the interpolations needed for reversion to Eulerian coordinates is held fixed as well. Furthermore, we assume that the radius of convergence R inf (t) of the time-Taylor series around time t does not vary significantly in the interval [0, T ] and we denote it R. As we will see in Section 6, at least in two dimensions, R inf (t) does not decrease in time dramatically, and hence this is a reasonable approximation.
We prescribe the accuracy ε and apply the condition (16). The L 2 -norm series ∞ s=0 ξ (s) (∆t) s , which has positive coefficients and a radius of convergence R, can be approximated, at least for large orders, by the geometric series C ∞ s=0 (∆t/R) s , where C is a positive constant. This exponential dependence on the order, is well supported by the numerical results on high-order time-Taylor expansions of Section 5 (see, in particular, Fig. 11). For large truncation orders S, the number of operations is O(S 2 ) and the number of time steps for integration to time T is T /∆t; (16) and (17) yield ∆t = R(ε/C) 1/S . Hence, the number of operations needed to reach a given time is proportional to S 2 (ε/C) −1/S . Minimising this expression over S, we find the optimal order S ≈ −(1/2) ln(ε/C) and ∆t = Re −2 . For small S and large N, the complexity SN D ln N of FFTs needed to solve the Helmholtz-Hodge problems for ξ (s) becomes dominant. By similar arguments, the optimal order is then S ≈ − ln ε. When both complexities are essential, the optimal order of expansion, S, satisfies the inequality

Testing the Cauchy-Lagrangian numerical method in two dimensions: Comparison with Eulerian simulations
In Section 4.1 we specify the test flows and describe three Eulerian algorithms that have been used for comparison with the Cauchy-Lagrangian algorithm. Validation, with emphasis on accuracy, is presented in Section 4.2. Efficiency of the CL algorithm is discussed in Section 4.3. In Section 4.4 we discuss spatial truncation artefacts. All computations reported in this section are in double precision.

Flows and numerical methods used for validation
For the reason explained in Section 3.2, all our tests of the CL method have been done using 2D flows with analytic initial data having non-trivial dynamics in both Eulerian and Lagrangian coordinates. Two different flows, called the test flows, were used as initial condition. The first one is a very simple flow, here called the "4-mode" flow, with the initial vorticity ω (init) = cos x + cos y + 0.6 cos 2x + 0.2 cos 3x.
The second one is a particular realisation of the so-called "random" flow [45], used in Section II.C of [39], where the time evolution of the latter flow is presented in detail. The characterisation of the random flow is best done in the Fourier space, where each wave vector consists of a pair of signed integers k ≡ (k 1 , k 2 ), conventionally combined into shells K ≤ |k| < K + 1 where K is integer. Each such shell involves N(K) Fourier harmonics. For all k in the Kth shell, the Fourier coefficients ω k of the initial vorticity are assigned the same modulus 2K 7/2 exp(−K 2 /4)/N(K) and phases that are uniformly and independently distributed in the interval [0, 2π], except that the coefficients for opposite wave vectors are complex conjugate, so that the flow is real. We have used four programs to solve the 2D Euler equation, realising the CL algorithm with the Taylor series truncated to order S (denoted CLS), described in Section 3, and three algorithms working exclusively in Eulerian variables, namely, the Runge-Kutta algorithms of orders two and four (denoted RK2 and RK4, respectively) and the algorithm relying on the Eulerian time-Taylor expansion truncated to order S (denoted ETS). The CL algorithms used in this section are CL8, CL16 and CL24 with time steps, chosen as explained in Section 3.5, and with the accuracy in (16) set to ε = 10 −12 . Fig. 4 shows the evolution of the time steps for the test flows. In principle, they are allowed to vary in time, but, as seen, they change very moderately. This is apparently related to the very slow temporal variations of the radii of convergence of the time-Taylor series, discussed in Section 6.
Our calculations have mostly been performed with spatial resolutions of 1024 2 and 2048 2 Fourier harmonics before dealiasing (performed using the standard 2/3-rule applying to quadratic problems). Occasional calculations had resolutions of up to 8192 2 harmonics.
Concerning the standard RK2 and RK4 algorithms, we only comment on the values of the time steps, taken constant in our calculations. For RK2, the smallest time step that we were able to run in reasonable CPU time was ∆t = 10 −4 . As we will see, this has yielded accuracies in the range 10 −6 − 10 −9 , significantly worse than the ones obtained with the other algorithms, requiring comparable or smaller CPU time. For RK4, the time steps used are listed in Table 1. They have been chosen as the largest value that gives stable integration, so that doubling of the time step causes numerical overflow. The corresponding Courant numbers (based on the maximum velocity and the maximum wave number after dealiasing) are indeed slightly in excess of unity.
The Eulerian Taylor (ET) algorithm deserves more explanations. By analogy with time integration exploiting analyticity of Lagrangian fluid particle trajectories, one can perform Eulerian time integration exploiting the analyticity in time of the velocity and the vorticity when the initial data are analytic in space and remain so for at least a finite time. The ET algorithm was actually applied long ago to study real and complex time singularities of the Taylor-Green 3D vortex flow [6]. This was done with a single time step from 0 to t, followed by attempted analytic continuation in t beyond the disk of convergence using, e.g., Padé approximants. Unfortunately, this method failed to give clear evidence for or against finite-time blow up in 3D. Here, we are interested in the 2D multistep form of the ET algorithm applied to the vorticity formulation of the Euler equation where ω denotes the only component of the vorticity and e z is the unit vector in the z-direction. We assume that the vorticity field ω is known at time t, and Taylor-expand ω(t + ∆t) in powers of ∆t: After substitution in (20), we obtain the following recurrence relations (where ω 0 ≡ ω(t)): Time-stepping is done, as in the CL method, by applying the Taylor series truncated to a finite order S and choosing an appropriate time step ∆t. Here, we report results obtained with an 8-term truncation, called ET8. Higher-order truncations have special problems, that will be discussed in Section 5. Given the analyticity in time of the Eulerian solution (for all times in 2D and a finite time in 3D [3]), we certainly must require that the time step be smaller than the (Eulerian) radius of convergence of the Taylor series. However, no Eulerian analogue of the condition (16) is available, and much smaller time steps must be chosen; typically Courant numbers order unity are needed. The reason, however, is not stability, as in low-order schemes, but the rounding errors, as explained in Section 5. In practice, for ET, we use the time steps of the same length as for RK4, given in Table 1.
Finally, a few words about spatial truncation errors. Since we restrict ourselves to analytic initial data, the solutions are also analytic in both Lagrangian and Eulerian coordinates. The Fourier transforms of such functions decrease exponentially with the wave number k as e −δk , where δ is the radius of the cylinder of analyticity (also sometimes called "width of the analyticity strip"). As shown in [41], δ can be estimated from the high-K asymptotics of the vorticity spectrum (also called "enstrophy spectrum"),  (a) (b) where ω k are the spatial Fourier coefficients of the vorticity ω. The large-wave-number behaviour of the vorticity spectrum is typically Here, the distance δ tells us how far away from the real domain the nearest complex singularities are, while the exponent α contains information about the type of singularities [35]. As a consequence of (24), the relative spatial truncation errors are typically estimated as e −δkmax , where k max is the maximum wave number after dealiasing, k max = N/3, where N is the number of grid points per spatial period. The numerical determination of δ from the vorticity spectrum is done by the same fitting method as for the numerical determination of radii of convergence (see Appendix B). When using the Cauchy-Lagrangian method, at each new time step, the Lagrangian and Eulerian descriptions coincide and, as a consequence, the Lagrangian and Eulerian δ are very close. We determine from Fig. 5 that, at the largest times shown for both flows, δ is about 10 −2 . Since k max = 8192/3 ≈ 2731, the resulting relative truncation error is about e −δ(t)kmax ≈ 2 × 10 −12 . (The product δ(t)k max is sufficiently large to use this asymptotic wave number formula.) We also note that for the lowest resolution used in our simulations, namely, 1024 2 harmonics, the 4-mode flow achieves comparably small spatial truncation errors only up to t = 2. If, however, we just request a level of spatial truncation error of about 10 −4 , which is not visually detectable on contours of the vorticity and of its Laplacian, we can extend the computations to about t = 4. We will come back to this in Section 4.4, when discussing truncation artefacts.

Validation of the CL algorithm and accuracy of agreement
Before comparing the results for the different methods, it is of interest to display the evolution of the test flows graphically. For the random flow, this was done in Section II.C of [39]. Here we focus on the 4-mode flow.
The evolutions of the vorticity and of its Laplacian for the 4-mode flow are shown in Figs. 7 and 8, respectively. The flow organises itself into a periodic array of large eddies centred at the centre and corners of the displayed periodicity box. Fine structures involving large vorticity gradients develop at the boundaries between these large eddies.
At t = 2, we see in the plot of vorticity isolines, and more clearly in the isolines of the Laplacian two rather fine-scale structures (related by the central symmetry of the flow) and two other ones (roughly perpendicular to the former), which display only somewhat larger scales. By t = 4, these structures have thinned to comparable scales. The thickness of such structures can be interpreted as a rough measure of the distance to complex-space singularities, because it determines the strength of spatial derivatives in the real-space domain. So, here we have a pair of singularities at complex-conjugate spatial locations that competes with another similar pair, with eventually the one pair furthest from the real domain catching up with the one which was closer. Such a phenomenon naturally leads to interference patterns in the vorticity spectrum. This apparently explains the bending of the plot in Fig. 5 and the rippling seen in Fig. 6 at t = 4.
We now show that the CL8 results agree with those obtained by traditional methods for both test flows within the time intervals where spatial truncations effects are negligibly small. Table 2 shows the vorticity discrepancies between results of different algorithms, namely, the maximum over space of the difference of the vorticity calculated with two different codes. This is done for the two test flows at different spatial resolutions and different output times. We also performed energy and enstrophy conservation tests, whose results are presented in Tables 3 and 4, respectively. (We will return to conservation issues in Section 4.4.) The main result is that the vorticity discrepancies between the CL results and those obtained by the standard Eulerian RK4 method are very small, around 10 −11 -10 −12 , except for later times, when spatial truncation errors become important. Note that the discrepancy with the RK4 results is roughly the same for CL8, CL16 and CL24. Choosing between these methods is thus only a matter of efficiency, as discussed in the next section.

Efficiency of the CL method
Before discussing the relative performances of the various algorithms, let us note that most of the computations were made on a computer -marketed in 2011 -that has an Intel Core i7-3960X processor (64 bit data width, frequency 3.3 GHz up to 3.9 GHz in the Turbo regime, 15 Mbyte L3 cache). It has 6 cores with hyper-threading, which does not actually matter for our application, since we have so far experimented only with sequential versions of the codes. Parallelising such a code is standard and requires the same tools as for the parallelisation of any spectral code.
The CPU timings for the integration to various output times for our codes are shown in Table 5. The timings for the Lagrangian CL codes are way below those for the best Eulerian code: CL24 needs at least an order of magnitude less CPU time than the commonly used RK4; the ET8 algorithm is at least 2 times slower than RK4, and RK2 is much slower.
We also observe that CL16 and CL24 have comparable efficiency and are significantly faster than CL8. By the arguments of Section 3.5, the optimal number of terms in the time-Taylor series is      around S = 14. But these arguments did not account for the fact that the operation count for FFTs has an extra logarithmic factor (of the size of the arrays). An improved estimate pushes the optimal S to somewhere between 16 and 24. Let us also discuss the durations of the runs for the extreme resolution of 8192 2 harmonics, used to calculate the radii δ of the cylinder of analyticity, shown in Fig. 5. Since in a Lagrangian code the time step does not depend on the spatial resolution, we expect even larger CPU gains than reported above. This is indeed the case. For the 4-mode initial condition, CL20 runs about 120 times faster than RK4, and for the random initial condition about 50 times faster than RK4. Actually, our CL runs took a few days on the machine described at the beginning of this section. We did not attempt to carry out the RK4 computations, since this would have taken many months. Of course, the huge increase in speed, using CL method, is taking full advantage of the spatial 2π-periodicity of our problem.

Galerkin truncation artefacts
In spectral simulations of ideal flow, one typically runs out of spatial resolution after some time of integration. As pointed out in Section 3.2, for analytic initial data, spatial truncation errors roughly behave as e −δ(t)kmax , which ceases to be small when complex-space singularities approach the real domain to within a few spatial grid mesh sizes. In Eulerian simulations, localised artefacts, in the form of spurious oscillations of contour lines, called "tygers," are then appearing in unexpected locations [39]. These are produced by truncation-generated waves interacting resonantly at places where the fluid velocity matches that of developing fine-scale structures. For the 4-mode initial condition, the phenomenon is illustrated in Fig. 9, which shows Eulerian RK4 simulations with the resolution of 1024 2 harmonics at the earliest times when tygers become visible. As is typical for small-scale objects, the tygers become at first conspicuous in contour levels of the Laplacian of the vorticity, around t = 4 when δ is about 5 spatial mesh sizes, before they are seen in the vorticity itself, around t = 5 when δ is about 1.5 spatial mesh sizes. For the random initial condition at the same resolution (not shown) similar phenomena happen around t = 0.71 and t = 1.
Do analogous artefacts emerge in our Cauchy-Lagrangian simulations? Contours of the Laplacian of the vorticity show similar spurious oscillations (for the same flow and the same spatial resolution) around t = 4.1 (see the last panel in Fig. 8). However, the vorticity shows nothing comparable until δ drops well below one spatial mesh.
Hence, Galerkin truncation artefacts are sharply reduced in our CL simulations, compared to Eulerian ones with the same resolution. This is surprising given that the Euler and Lagrangian coordinates are resynchronised at each time step, so that not much difference in behaviour and deterioration in time of solutions would be expected. Part of this reduction could be due to the smoothing effect of the interpolation stage, also observed in the semi-Lagrangian approach of [32]. There is, however, another possibility: solutions truncated to a finite number of spatial Fourier modes (generally referred to as Galerkin truncated) behave differently in Eulerian and Lagrangian coordinates. In Eulerian coordinates, it is easily shown that quadratic invariants (such as energy and enstrophy, in 2D, or energy and helicity, in 3D) remain invariant after Galerkin truncation. Hence, as observed in [39], energy and enstrophy, which in the absence of truncation would be transferred to smaller scales, have to be transferred somewhere else and get pushed into the tygers. In Lagrangian  coordinates, the behaviour is different: The energy is E = (1/2) |ẋ(a, t)| 2 d 2 a, where the integral is over the [0, 2π] × [0, 2π] periodicity domain. Before truncation, the energy is exactly conserved in any time interval ∆t. After Galerkin truncation, it is only conserved up to second order in ∆t, and the same holds for enstrophy.

Rounding noise in high-order Lagrangian vs Eulerian time-marching
Some challenging open problems in hydrodynamics, such as the issue of blow up for 3D flow [22] or for axisymmetric 2D flow [31], may require extremely accurate simulations in order to distinguish genuine blow up from simulation artefacts. This means that spatial and temporal truncation errors must be kept very small and the precision very high.
For flow in a domain of sufficiently simple geometry (foremost space-periodic flow), the most accurate methods for handling the spatial discretisation are of the spectral or pseudospectral type [23]. The reason is that, for spatially analytic flow, the Fourier coefficients fall off exponentially at large wave numbers [17,41] and, thus, the truncation error drops exponentially when increasing the resolution. Unfortunately, (pseudo-)spectral simulations hardly ever use temporal discretisation (timemarching) of order higher than 4, the most frequent being Runge-Kutta of fourth order (RK4). We have therefore a striking imbalance between temporal and spatial discretisation. Of course, when using Eulerian methods, one generally has to satisfy the Courant-Friedrichs-Lewy [12] condition: where k max is the maximum wave number (with dealiasing taken into account), U max is the maximum velocity in physical space, and C is an order unity constant, depending on the numerical scheme. As a consequence, high spatial resolution requires a very small time step, so that a fourth-order scheme might be deemed sufficiently accurate. The situation is drastically different when working in Lagrangian coordinates. The time step is not controlled by a CFL condition any more, but solely by the radius of convergence of the time-Taylor expansion for the Lagrangian map. This radius is roughly the inverse of the largest velocity gradient present and does not depend on spatial resolution. Hence, in principle, there is no reason to limit oneself to low-order time-marching.
We now report results of experimenting with very high-order time-Taylor expansions in both Lagrangian and Eulerian coordinates. We will see that the severest limitation stems from the growth of rounding errors.
In Fig. 10, the L 2 norms of time-Taylor coefficients in Lagrangian coordinates, obtained by the Cauchy-Lagrangian method, and in Eulerian coordinates, obtained by the Eulerian Taylor method, are shown for the 4-mode initial condition, computed in double and quadruple precision up to the orders beyond which an unnatural behaviour, dominated by rounding errors, sets in. Taylor coefficients of the displacement, ξ (s) , and of the vorticity, ω s , in the expansions of the Lagrangian and Eulerian solutions, respectively, have been computed around t = 0 with two resolutions of 512 2 and 1024 2 harmonics.
The most striking feature is that the L 2 norms show a range of orders over which the results display very little dependence on precision and then a rapid transition to a regime in which drastic differences of many orders of magnitude are observed. This happens around a "transition order" that depends very much on precision and which is considerably larger (by about a factor 4) in the Lagrangian case than in the Eulerian case. (There is also some dependence on the resolution.) The exponential behaviour (straight line in the lin-log coordinate system used) is not observed at the smallest values of the order s. Indeed, such exponential behaviour ∼ R −s of the coefficient of order s of a power series with radius of convergence R is typically only an asymptotic property at large s.
The large difference between the Lagrangian and the Eulerian computations in the transition order, where rounding noise takes over, can be seen as a consequence of a phenomenon discovered by Ebin and Marsden [14], which is called the loss of derivatives in the Eulerian formulation. In this formulation, each time derivative ∂ t is accompanied by a space derivative v · ∇. As a consequence, if the initial velocity has only a finite number of space derivatives, it will also have a finite (generally, the same) number of time derivatives and, thus, cannot be analytic in time. In Lagrangian coordinates, this is not the case: time-analyticity of the Lagrangian map, for at least a finite time, holds if the initial velocity is slightly smoother than once differentiable (see Section 2). Ebin and Marsden, although they did not prove time-analyticity in the Lagrangian formulation, did show that the variational formulation of the Euler equation, which is essentially Lagrangian, allows interpreting it as an ordinary differential equation in function spaces, such as Sobolev spaces, possessing a finite number of space derivatives. The proof was given using the language of differential geometry on infinite-dimensional Riemannian manifolds. Although somewhat simplified by Bourguignon and Brezis [4], the proof did not suggest any numerical implementation. This became possible with the rediscovery of Cauchy's Lagrangian formulation, enabling us to generate time-Taylor coefficients of arbitrary orders. Now, we comment on why loss of derivatives in Eulerian coordinates is responsible for rounding noise that grows very fast with the order of the time-Taylor coefficient. The 4-mode initial condition used for the present study has of course space derivatives of arbitrary order (it is an entire function). Hence time derivatives of arbitrary orders at t = 0 also exist. It is easy to see that spatial Fourier harmonics of the time-Taylor coefficient of order s (not subjected to Galerkin truncation) vanish beyond a wave number growing linearly with s. (This is true both in Lagrangian and Eulerian coordinates.) Once the problem is discretised in space and FFTs are used, those wave vectors which should not be excited will be populated by rounding noise. In the Eulerian case, this noise is amplified by the presence of the space derivatives in the recurrence relations, each yielding a factor ik on the Fourier coefficients. In the Lagrangian case, the recurrence relations are all written in terms of gradients of time-Taylor coefficients, but otherwise we only perform additions, multiplications and apply bounded Calderon-Zygmund operators. Thus, amplification of noise at high wave numbers is not so strong as in the Eulerian recurrence relations. A further confirmation of this mechanism is obtained by considering the relative change with the resolution of the order s at which the transition to noise takes place. When changing from resolution of 512 2 to 1024 2 harmonics, the transition orders change by just a few percent in the Lagrangian calculations and by nearly 20% in the Eulerian calculations.
The advantage of Lagrangian time-marching over Eulerian time-marching from the point of view of rounding noise is not limited to initial conditions comprised of a finite number of non-vanishing harmonics. To show this, we use a multistep method with the same initial data as above, but where a very broad band of Fourier harmonics become excited beyond the first time step. We begin by explaining how we choose our time steps in the Eulerian calculations. The initial data considered here are analytic in the space variables and remain analytic in space and time as long as there is no blow up (for ever in 2D) [3]. At any time t > 0 there is a finite radius of convergence R(t) of the Eulerian time-Taylor series giving the solution (say, the vorticity) at time t + ∆t in terms of the solution at time t. Certainly, we have to take ∆t < R(t), but this is not enough, because of a rounding noise problem, closely linked to what we discussed above.
In Eulerian coordinates, the fastest temporal variation of small-scale eddies comes from them being swept by the large-scale energy-carrying eddies. Denoting by ϑ(t) the amplitude (say, the vorticity) associated to such an eddy, the sweeping can be modelled by the equation ∂ t ϑ + (U · ∇)ϑ = 0, where U is a large-scale velocity which can be taken uniform to leading order (the expansion parameter being the ratio of scales of the swept eddies and the sweeping eddies). Let us denote by ϑ k the Fourier coefficients of the small-scale vorticity. The sweeping equation has the solution ϑ k (t + δt) = exp(−ik · U ∆t)ϑ k (t). Expanding the exponential in a Taylor series in ∆t, we find: We define C(k, U ) ≡ |k · U ∆t|. Its maximum over all wave vectors and velocities present in the flow is the Courant number Co ≡ k max U max ∆t. Here, k max is the maximum wave number, determined by the number of grid points per spatial period, N, and the dealiasing (typically, k max = N/3), and U max is the maximum velocity. We observe that, in the series (26), the moduli of the various terms decrease with s when Co < 1. Otherwise, for sufficiently small eddies, swept by sufficiently fast large eddies, i.e., for C(k, U ) > 1, the moduli increase at first with the order s while s < C(k, U ). Of course, the whole sum, being an imaginary exponential, has a unit modulus, but this is a sum of a large number of terms, some of which have moduli that are exponentially large in the Courant number. The accuracy of such a summation for large Co is low, due to the finite precision of the computations. This can result in catastrophic loss of precision -a well-known phenomenon that can mar numerical summation of real series with alternating signs (see, e.g., [1]). To avoid this difficulty, one can require Co < 1, which is precisely the standard Courant-Friedrichs-Lewy condition [12]. Note that the CFL criterion for the Eulerian Taylor method does not stem from a stability condition, but from the finite precision of the computations. Fig. 11 shows the L 2 norms of time-Taylor coefficients of various orders, in Lagrangian and Eulerian coordinates, respectively. All computations are done in quadruple precision with the resolution of 1024 2 Fourier harmonics. The various graphs are plotted for Taylor expansions around different output times. Because of the constraints on time steps that are quite different in Lagrangian and Eulerian coordinates, the number of time steps between successive output times is quite different in Lagrangian coordinates (typically of the order of tens) and in Eulerian coordinates (typically of the order of hundreds). While the Lagrangian and Eulerian time-stepping is performed using Taylor expansions truncated at the eighth order, at output times they are pushed to much higher orders. For the Eulerian integrations we observe again a transition to rounding-noise dominated coefficients at orders roughly in the range 15-20 (a little beyond order 20 for t = 3). In the Lagrangian case this also happens but, at orders about four times larger (not shown in the figure). Thus, we observe essentially the same phenomenon as for the Taylor expansions around t = 0. The explanation is now somewhat different: in the Eulerian case it is not just the noise that is amplified by the ik factors stemming from spatial differentiation, but the actual signal, which can now be much stronger then the direct rounding noise. However, when we use FFTs, there is indirect rounding noise, stemming from much lower wave numbers that have much higher amplitude. Furthermore, in the Eulerian simulations, this indirect rounding noise gets again amplified by the ik factors stemming from the space derivatives.

The radius of convergence and depletion
A well-known advantage of Lagrangian and semi-Lagrangian methods is that the time step is constrained not by the usual CFL condition but, roughly, by the inverse of the largest velocity gradient. For our Cauchy-Lagrangian method, the optimal time step ∆t for a given accuracy, when starting from time t, is a fixed fraction of the radius of convergence R inf (t) of the time-Taylor series for the Lagrangian map (see Section 3.5). It is of interest to know how R inf (t) changes in time, since this determines how the time step should be changed and also the total number of time steps that will be needed to reach a given termination time t fin . When using initial conditions that have mostly low-wave-number harmonics, such as the 4-mode flow, the subsequent transfer of excitation to higher harmonics will be accompanied by some growth in the velocity gradient: substantially so in 3D and more moderately in 2D, because of the conservation of vorticity.
As a consequence, we expect a decrease of the radius R inf (t). Clearly, in 3D, when R inf (t) → 0 as t → t ⋆ from below, a blow up of the solution takes place at the finite time t ⋆ . Whether the converse holds is so far unknown. Anyway, in the present numerical studies, we are only concerned with the 2D case. Before turning to numerical results, let us recall the proven constraints. A theorem of Wolibner [42] then guarantees regularity for all t > 0 when the initial vorticity is Hölder-continuous of exponent α 0 > 0. More precisely, Wolibner showed that, at time t > 0, the vorticity is Höldercontinuous of exponent α(t) ≥ α 0 e −t sup |ω 0 | , where sup |ω 0 | is the supremum of the modulus of the initial vorticity. The constraint on vorticity conservation does not prevent the velocity gradient from growing, but it cannot grow faster than e t sup |ω 0 | , a bound which is actually sharp [27]. Consequently, for such 2D flow, the radius of convergence R inf (t) has a lower bound proportional to e −t sup |ω 0 | (see Section 2). The actual numerical results, shown in Fig. 12, indicate that the decrease in time of the radius of convergence is much slower than this bound. Consider, for instance, the 4-mode initial condition, for which the supremum of the initial vorticity is 2.8 . The aforementioned lower bound could have the radius of convergence decreasing by a factor e −5×2.8 ≈ 8 × 10 −7 over the time span from 0 to 5. Actually, the radius of convergence decreases only by about 30%. As to the random initial condition, its radius does not decrease at all. For comparison, the radii of convergence for the Eulerian Taylor method are also shown: they decrease much faster. This huge discrepancy between the proven bounds and the actual results is most likely related to the phenomenon of depletion: in both 2D and 3D incompressible flow it is frequently found that the nonlinear effects are growing in time much slower than permitted by the best proven bounds [19], probably because such flow tends to organise itself into structures that have an almost vanishing nonlinearity, e.g., quasi-one-dimensional flow depending mostly on a single coordinate. Depletion is definitely why 2D flow is often found numerically to be even more regular than implied by Wolibner's lower bound, and also why 3D simulated flow show little reliable evidence for blow up (one possible exception is a recently studied swirling axisymmetric flow with a solid boundary [31]).

Conclusion
We have developed a novel Cauchy-Lagrangian (CL) method to solve the initial-value problem for ideal incompressible flow, governed by the Euler equation. The method operates in Lagrangian coordinates, relies on temporal analyticity of particle trajectories in time and employs simple recurrence relations for coefficients of the Lagrangian time-Taylor series, stemming from the Cauchy invariants. Our experiments in 2D show that this method is highly efficient (like many semi-Lagrangian methods, it does not suffer from the Courant-Friedrichs-Lewy constraint) and enables typically an order of magnitude faster time-marching than the fourth-order Runge-Kutta method and the Eulerian Taylor method, at the same time being significantly more accurate. As shown in Section 4.3, the gain in speed can even be higher, in the range 50-100, at the highest spatial resolutions. Last but not least, the Cauchy-Lagrangian method allows very-high-order time marching without the severe rounding error problems that affect high-order Eulerian time marching. Of course, these impressive results are so far obtained for just a narrow niche of problems: at present, incompressible space-periodic two-dimensional flow.
The CL algorithm has been presented in Section 3 for flow of arbitrary dimension. Since the CL method allows efficient high-order and high-precision simulations with only modest growth of rounding errors, it is a strong candidate for the reliable numerical exploration of the blow up problem for the 3D Euler equation (see, e.g., [22]).
The method, here described for incompressible Euler flow, can be readily extended for solving certain other problems arising in various areas of physics. These are described by equations for which known theorems guarantee the time-analyticity of Lagrangian fluid particle trajectories. Analyticity in time of the Lagrangian map for at least a finite time, under the assumption that the initial velocity gradient is, e.g., Hölder-continuous, holds for the Euler equation (which applies to inviscid incompressible fluid flow and to the 2D drift-Poisson equation for an electron plasma). Furthermore, it also holds for the surface quasi-geostrophic equation, the incompressible porous medium equation and the Boussinesq equation. So far, the latter three have been handled not by a recurrence relation approach, using Cauchy's Lagrangian formulation, but by applying the Faà di Bruno formula [11] which gives an explicit but somewhat cumbersome representation of all the time-Taylor coefficients [10]. Similar analyticity results with simple recurrence relations were also proved for compressible flow of cosmological interest, such as the Euler-Poisson equations for an Einstein-de Sitter universe [44], recently extended to the more realistic case of a ΛCDM universe [38]. The Euler-Poisson equations apply to the collisionless flow of dark matter only as long as multi-streaming (development of caustics with more than one velocity) is absent or negligible. In such a framework, our methods are useful in at least two cases: the blow up problem (called shell-crossing by cosmologists) and the reconstruction problem, in which the whole dynamical history of the universe is obtained from just the spatial distribution of matter at the present epoch [18,7].
But a warning is necessary: not all ideal fluid flow problems are endowed with time-analytic Lagrangian fluid particle trajectories. For example, this is apparently not the case of magnetohydrodynamic (MHD) ideal flow and of barotropic flow. For non-diffusive MHD flow, one can still write recurrence relations for the time-Taylor coefficients of the Lagrangian map. However, the Lagrangian gradient of a time-Taylor coefficient of a higher index does not involve just the gradients of those of lower index, but also their space derivatives of order up to three. In other words, in contrast with the ordinary Euler equations which do not lose derivatives in Lagrangian coordinates (see Section 5), the MHD equations do, and thus MHD flow with only limited initial spatial smoothness lacks analyticity of fluid particle trajectories. Such a loss of derivatives may be viewed physically as due to the frozenness of the magnetic field. The barotropic equations govern compressible fluids, in which the pressure is a function of density only. They play an important role in geophysical fluid dynamics [33]. The standard derivation of the Cauchy invariants applies also to the barotropic equations. However, the second equation, stemming from mass conservation, is profoundly modified and, in it, loss of derivatives and hence of analyticity occurs again.
Another class of flows which might be amenable to some variant of the Cauchy-Lagrangian method is viscous flow (e.g., Navier-Stokes flow) that is initially analytic in the space variable (if it is not space-analytic initially, it becomes so at an arbitrary short positive time [16]). Such flow is also analytic in time in both Eulerian and Lagrangian coordinates. There is however a new practical problem: when the equations are written in Lagrangian coordinates, if we try following Cauchy's approach, we find that viscous terms involve nonlinearities in the Lagrangian map of degree four in 2D and of degree 6 in 3D; thus the recurrence relations become somewhat cumbersome. An alternative and much simpler strategy for viscous flow is a fractional step method [43] in which the reversion to Eulerian coordinates by interpolation is followed by a solution of the viscous heat equation. At high Reynolds numbers, this will result in a low-order scheme for the smallest viscosity-dominated scales but the larger scales will still benefit from the high accuracy permitted by the CL method.
In this first paper on the Cauchy-Lagrangian method we have concentrated deliberately on explaining and validating the technique. We expect that, over the years, a number of researchers will apply and improve the technique for a range of problems. One of them is of course the outstanding and very difficult issue of blow up for Euler flow. Real blow up cannot take place in 2D, but there is already strong evidence that in 2D the Lagrangian time-Taylor series (around t = 0) has a finite radius of convergence. Hence, Lagrangian trajectories have complex-time singularities. It may be important to determine what happens at such complex blow ups.
We consider an abstract infinite function series where a is in the d-dimensional periodicity domain (torus) T d . This need not be a Lagrangian solution of the Euler equation; in particular, it suffices for our purposes to provide arguments for a scalar-valued series. Let R(a) denote the radius of convergence of the time series at point a. We define R inf ≡ ess inf a R(a); we denote the L p norm of a function ξ(a) by ξ(a) p and the radius of convergence of the series ∞ s=1 ξ (s) (a) p t s by R p . The Lebesgue space L p (T d ) under consideration is that of scalar-valued functions of the spatial variable a of an arbitrary dimension d; we assume that all ξ (s) (a) and the sum ξ(a, t) are space-periodic.
We prove now the following Theorem. For some p ≥ 1, let all ξ (s) (a) be in L p . Then 1. For all t < R p , the sum ξ(a, t) belongs to L p , and R inf ≥ R p . 2. If ξ(a, t) is also in L p for any t < R, then R p ≥ R inf . Proof. The radius of convergence of the series (A.1), R(a), does not change, when the coefficients ξ (s) (a) are replaced by their absolute values, |ξ (s) (a)|. Similarly, applying the absolute value does not modify the Lebesgue space norms. Consequently, without any loss of generality, we assume ξ (s) (a) ≥ 0 for all s and a, otherwise replacing the coefficients by their absolute values.
To prove statement 1 of the theorem, we consider the truncated series ξ S (a, t) = S s=1 ξ (s) (a)t s for a fixed t in the interval 0 < t < R p . Since the partial sums ξ S (a, t) constitute a Cauchy sequence in L p , it converges in L p to a limit function ξ(a, t) also belonging to L p . By a well-known theorem of functional analysis (see, e.g., [2]), one can choose a subsequence ξ S k (a, t) that converges almost everywhere to ξ(a, t). (In other words, for our fixed t, there exists a set of measure zero such that, at any point a in the complement of this set, lim k→∞ ξ S k (a, t) = ξ(a, t), the limit being here understood as the usual limit for infinite sequences of real numbers.) However, since ξ (s) (a) ≥ 0 for all s ≥ 1, the partial sums ξ S (a, t) grow monotonically at each point a (for the fixed t) on increasing the index S, and hence the entire sequence ξ S (a, t) converges for S → ∞ almost everywhere to ξ(a, t). Therefore, any fixed time t < R p is within the disk of convergence of the series (A.1), i.e., R(a) ≥ R p for almost all a, which implies the first statement of the theorem.
We prove now statement 2. For any p ≥ 1 and any quantities α s ≥ 0 such that ∞ s=1 α s = 1, clearly, ∞ s=1 α p s ≤ 1. Taking α s = ξ (s) (a)t s /ξ(a, t) (recall ξ (s) (a) ≥ 0 for all s), we thus find for almost all a and any t such that 0 < t < R(a), and hence for any 0 < t < R inf . Integrating this inequality over the periodicity domain yields s≥1 ξ (s) (a) p p t sp ≤ ξ(a, t) p p . For t > R p = (lim sup s→∞ ξ (s) (a) 1/s p ) −1 , the series in the l.h.s. of the latter inequality diverges, but it converges for all t < R inf , for which the inequality holds. Consequently, R inf ≤ R p . This concludes the proof of the theorem.
It remains to comment on the significance of the condition in statement 2 that ξ(a, t) belongs to L p for all t < R. One can indeed construct an abstract counterexample showing that the radii are not necessarily equal, if this condition is not satisfied.
Remark. The Lebesgue space L p (T d ) in the statement of the Theorem can be replaced by a Sobolev function space W p q (T d ) for p ≥ 1 and q > 0. Its norm is defined as ξ(a) q,p = (I − ∇ 2 ) q/2 ξ(a) p , where I is the identity operator and ∇ 2 the Laplacian (the spaces for p = 2 are often denoted by H q (T d )). The quantity R p is then replaced in the modified theorem by R q,p , the radius of convergence of the series ∞ s=1 ξ (s) (a) q,p t s , and instead of R inf we use R q,inf , the essential infimum over space (i.e., over the variable a) and over all q ′ , such that 0 ≤ q ′ ≤ q, of the radii of convergence of the series Convergence of the series ∞ s=1 ξ (s) (a) q,p t s implies convergence of (A.3) in L p for all q ′ ≤ q and, by the arguments in the proof of the Theorem, the inequality R q,inf ≥ R q,p . Conversely, convergence of (A.1) in W p q implies convergence of (A.3) in L p for q ′ = 0 and q ′ = q, and by the second part of the Theorem the radii of convergence of the two series, ∞ s=1 ξ (s) p t s , and ∞ s=1 (−∇ 2 ) q/2 ξ (s) p t s , are both bounded from below by R q,inf , whereby R q,inf ≤ R q,p . Note that R q,p and R q,inf are monotonically decreasing functions of q for any p ≥ 1. The inequality R q,p < R q ′ ,p for some q > q ′ signals existence of a ring within the disk of convergence in W p q ′ , where the series (A.1) does not converge in W p q . An intriguing problem is to establish whether for a solution of the Euler equation R q,inf = R q ′ ,inf and R q,p = R q ′ ,p hold for all sufficiently small q = q ′ .
Appendix B. Computation of the radius of convergence of a power series with positive coefficients The implementation of the Cauchy-Lagrangian method does not require the evaluation of the radius of convergence of the time-Taylor series employed. However, to test the method, to analyse its performance and interpret some of the results (Sections 4-6), we need to compute the radii of convergence of the time-Taylor series for both CL and ET algorithms. For estimating these, using the result of Appendix A, we can replace the full functional Taylor series ∞ s=1 ξ (s) (a) t s by the ordinary power series ∞ s=1 ξ (s) (a) t s , where · denotes the L 2 norm. To obtain reasonable precision on the radii, say, a few percent, this computation actually requires the use of quadruple precision and the knowledge of enough Taylor coefficients, usually more than is required for the time-stepping. Due to this added complexity, we do not compute the radii at every time step, but just frequently enough to be able to monitor how the radii evolve in time, e.g., every tenth time step.
In this Appendix we describe how we calculate the radii for such power series.
Here is the problem that we want to solve: Given a finite number, S, of coefficients of the power series where f s ≥ 0, we want to estimate its radius of convergence, which we denote by R. For our purposes, an accuracy of a few percent suffices.
To begin with, we recall some classical results. By the Cauchy-Hadamard theorem, when the limit exists. With less than a hundred coefficients of the time-Taylor series at hand, we can obtain the radius of convergence with a reasonable accuracy using such formulae only, if we know several terms in the large-s asymptotic expansion of the coefficients f s . We should derive the asymptotics of f s = ξ (s) (a) from the recurrence relations (7)-(8) for ξ (s) . Unfortunately, no such results are available so far. To proceed, we look at our problem from an entirely different angle and ask ourselves the question: Does the analytic function defined by the series ∞ s=1 ξ (s) (a) t s have the "generic" simple structure near the disk of convergence, such as many analytic functions (B.1) with real positive coefficients f s have?
A standard result is that the radius of convergence of a Taylor series is the distance in the complex t-plane to the nearest singularity(ies). If there is a single such singularity and the Taylor coefficients are all positive, the singularity is on the real positive axis at t ⋆ = R. The simplest case is when the singularity has a simple structure, e.g., a pole or a power-law branch point that can be characterised by a local behaviour ∝ (t ⋆ − t) ρ , where the exponent ρ is a negative real number. It may then be shown that the Taylor coefficients of large orders s have the following behaviour: where β = − ln R and α = −ρ − 1. (The same kind of asymptotics applies to Fourier series; see, e.g., [8].) In such a case, if we are able to extract the leading and subleading behaviour of f s , we may then infer not only the radius of convergence, but also the type of the nearest singularity on the real axis. One method is the Domb-Sykes plot of f s /f s−1 versus 1/s, where the y-intercept gives 1/R and the slope near the origin gives the exponent α [13,24]. We have applied all the standard methods, recalled above. We found that, when using 50-80 quadruple-precision Taylor coefficients from our L 2 series, such methods allow a determination of the radius of convergence with an error in the range 10%-20%. We also tried more advanced extrapolation techniques, such as convergence acceleration methods (see, e.g., [35,34] and references therein) and the asymptotic extrapolation technique [25,35,34], which have the potential of capturing both leading and subleading asymptotic behaviour. Unfortunately, all these techniques failed to reach the asymptotic regimes and thus gave no improvement over the more standard ones.
We found that for the L 2 series with 50-80 terms a rather straightforward method apparently gives a better accuracy. We again assume that the L 2 series has a single singularity and that (B.4) holds asymptotically for large s. It follows that the logarithms of the Taylor coefficients, F s ≡ ln f s , are given by F s ≈ ln γ + α ln s + βs for s → ∞.

(B.5)
We can determine the coefficients α, β and γ by a least-square fit over the available Taylor coefficients. That is, we minimise S s=1 d 2 s for the discrepancies d s ≡ c + a ln s + bs − F s . (B.6) We now give an example of such determination for the four-mode initial condition (19). Here f s are the L 2 norms of the first 80 time-Taylor coefficients ξ (s) for the displacement, computed at t = 0. Fig. B.13a shows the values F s /s and the fitted curve (c+a ln s+bs)/s. (Note that F s /s ≈ β = − ln R for large s.) The least square fit gives a = −1.94, b = −0.187 and c = 0.116. Thus, the radius of convergence is R = e −b ≈ 1.2 . To estimate the error on the radius R, we show in Fig. B.13b the scaled discrepancy d s /s for large s. We see that, beyond s = 20, the fluctuations do not exceed 1%. Similar determinations of the radius R(t) for the same initial four-mode flow and 0 < t < 5 have been made, using between 80 and 50 Taylor coefficients, yielding errors between 1% and 3%. We have used the same fitting technique to measure the radius δ of the cylinder of spatial analyticity around the real domain (in both Lagrangian and Eulerian coordinates). The analogue of (B.4) is then the following expression for the shell-averaged energy spectrum E(K), (B.7) How can we achieve a much better accuracy on radii of convergence? We need to know a lot more than 80 coefficients, say, a few hundred. Unfortunately, as observed in Section 5, in quadruple precision, rounding noise becomes explosively large just above order 80. We can go a little further by avoiding FFTs and performing the Fourier-space convolutions explicitly, but this is computationally very expensive. Otherwise, we have to resort to higher than quadruple precision. Recent algorithmic developments on FFTs in "medium" precision (up to 400 bits) offer interesting possibilities [26].