Structural stability for the splash singularities of the water waves problem

In this paper we show a structural stability result for water waves. The main motivation for this result is that we would like to exhibit a water wave whose interface starts as a graph and ends in a splash. Numerical simulations lead to an approximate solution with the desired behaviour. The stability result will conclude that near the approximate solution to water waves there is an exact solution.

The function c(α, t) is arbitrary since the boundary is convected by the normal component of the velocity of the fluid. Also, we notice that, in order to get an explicit equation for ∂ t ω, we need to invert the operator I + T = I + 2 BR(z, ·), z α and we have taken the acceleration due to gravity and the density ρ equal to one. Once one has solved this system for (z, ω) the velocity of the fluid and the pressure in the domain Ω(t) can be recovered by using Biot-Savart and Bernoulli laws. For details see [3].
In the last two decades these equations have been intensively studied. For an extensive survey about analytical results on water waves see the monograph [9].
In this paper we are concerned with the problem of the existence of water waves which start as a graph and become a splash curve in finite time. Roughly speaking, a splash curve is a smooth curve that collapses with itself in a single point such as the curve of fig.  1. A rigorous definition can be found in [3] where the existence of splash singularities has been shown. Coutand and Shkoller [5] have proven the existence of splash singularities in presence of vorticity. Fefferman, Ionescu and Lie [6] have proven the non existence of splash singularities for internal waves, i.e. for an interface between two incompressible fluids. We are interested in the following statement: Conjecture 1.1 There exist initial data z 0 (α), ω 0 (α) of solutions of the water wave equations such that at time 0 the curve z 0 (α) can be parameterized as a graph, the interface then turns over at a finite time T 1 > 0, and finally produces a splash at a finite time T 2 > T 1 .
We should remark that this conjecture is a combination of the scenarios in theorems [3,Theorem I.1] and [4,Theorem 7.1] and is supported by numerical evidence that we can see in Fig. 2. This numerical simulation was carried out using the method of Beale, Hou and Lowengrub [1].
(From now on we will omit the superscript tilde in the notation). We start computing a numerical approximation of a solution to the water waves equation 5 that starts as a splash, turns over and finally is a graph. Such a candidate is depicted in Fig. 2. With this aproximation we can construct explicit functions (x, γ) that solve the system where f and g are errors that we hope are small. By using the computer we are able to give rigorous bounds for these errors. The question we want to answer is if there exists an exact solution (z, ω) of the water waves equation close to these functions (x, γ). That means we need to prove the following theorem: where (x, γ, ψ) are the solutions of x β |x β | 2 dβ be γ t +2BR t (x, γ) · x α = −(Q 2 (x)) α |BR(x, γ)| 2 + 2bBR α (x, γ) · x α + (bγ) α − Q 2 (x)γ 2 4|xα| 2 α − 2(P −1 2 (x)) α + g ψ(α, t) = Q 2 x (α,t)γ(α,t) 2|xα(α,t)| − b s (α, t)|x α (α, t)|, (7) where (z, ω) are the solutions of (7) with f ≡ g ≡ 0, ϕ is the function ϕ = Q 2 z (α, t)ω(α, t) 2|z α (α, t)| − b(α, t)|z α (α, t)|, and E is the following norm for the difference Then we have that d dt E(t) ≤ C(t)(E(t) + E k (t)) + cδ(t) where ) and 2 (t)) 2 , k big enough depends on the norms of f and g, and E(t) is given by where the L ∞ norm of the function is the Rayleigh-Taylor function, and finally m(q l )(t) ≡ min α∈T |z(α, t) − q l | for l = 0, ..., 4, with which are the singular points of the transformation P .

Remark 1.3
We can absorb the terms in E(t) by E(t) raised to an appropriate power and terms in (x, γ) by performing the splitting z = z − x + x (or the analogous one for a different variable) for any norm or any quantity that appears in E(t). Theorem 1.2 was announced in [2].
If we knew C(t), f (t), g(t), k or bounds on them, a priori, then we could provide bounds on E(t) at any time T . We point out here that E(t) controls the norm ∂ α z 1 (α) − ∂ α x 1 (α) L ∞ . Let T g be a time in which the approximate solution is a graph, i.e. ∂ α x 1 (α, T g ) > 0 ∀α.
and this shows that z is a graph. In other words, the possible set of solutions of the water waves equation is a ball centered at (x, γ, ζ) with the topology given by E. All of the elements of this ball are graphs, therefore the solution is necessarily a graph. Thus, the problem is reduced to study and find bounds for C(t), f (t), g(t), k.
The recent developments of computer architecture have boosted their use in mathematics, giving birth to a full set of new results only achievable by this enormous power. However, it has the drawback that floating-point operations can not be performed exactly, resulting in numerical errors. In order to overcome this difficulty and be able to prove rigorous results, we use the so-called interval arithmetics, in which instead of working with arbitrary real numbers, we perform computations over intervals which have representable numbers as endpoints. On these objects, an arithmetic is defined in such a way that we are guaranteed that for every We can also define the interval version of a function f (X) as an interval I that satisfies that for every x ∈ X we have f (x) ∈ I.
The article is organized as follows: in sections 2 and 3 we give some details about how to control the errors f , g and the constants that arise in Theorem 1.2 by using the computer. Finally, in section 4 we give a complete proof of Theorem 1.2.
2 Bounds for f (t) and g(t)

Representation of the functions and Interpolation
The first thing one has to decide is how to represent the data and how to pass from the cloud of points in space-time obtained by non-rigorous simulation to a function defined everywhere in [−π, π] × [0, T ]. We need to interpolate in some way.
In our case, we chose to represent the functions x and γ by piecewise polynomials (splines) of high degree (10) in space, and low degree (3) in time. To do so, we first interpolate in space for every node in the time mesh. The interpolation is made via B-Splines. Since the interpolation is reduced to solve a linear (interval) system Ac = y, where A is constant in time and space and y depends on the values of the function at time t since the mesh in space is constant, we precondition by multiplying by the non-rigorous inverse of the midpoints of the entries of A. We remark that the system is interval-based because we need to produce a curve that is a splash (i.e. there have to be two points α 1 , α 2 such that we can guarantee x 0 (α 1 ) = x 0 (α 2 ). Finally, the system is solved using a rigorous Gauss-Seidel iterative method. We also remark that the need for interval-based calculations is only strictly necessary at time t = 0 since it is the only point in which we have to guarantee some equality. By working with multiprecision (1024 bits) we can get widths in the coefficients of the order of 10 −300 . In order to perform interpolation in time, we fix the values of the function and its time derivative at the mesh points. This gives us lots of systems of 4 equations (the values of the function and its derivative at both endpoints) and 4 unknowns (the 4 coefficients of the degree 3 polynomial) but with an explicit formula for each of them. With this method, our spline will be C 1 in time but it might not be C 2 .

Rigorous bounds for Singular integrals
In this section we will discuss the computational details of the rigorous calculation of some singular integrals. In particular we will focus on the Hilbert transform, but the methods apply to any integral kernel whose main singularity is homogeneous of degree -1. Parts of the computation (the N part) are slightly related to the Taylor models with relative remainder presented in M. Joldeş' thesis [8].
Let us suppose that we have a function f given explicitly by a spline (piecewise polynomial) which is C k−1 everywhere and C k except at finitely many points (the points in which the different pieces of the spline are glued together). We need to calculate rigorously the Hilbert Transform of f , that is and we want to approximate it by a piecewise polynomial function with less regularity, plus an error that can be bounded in H q , 0 ≤ q ≤ c < k and in L ∞ . Let us assume that the knots of the spline are α i , i = 0, . . . , N − 1 and that we fix x ∈ [α i , α i+1 ] where the indices are taken modulo N and the distance between the indices is taken over Z N . We can split our integral in Now, if we want to express Hf F (x) as a polynomial, it is easy since the integrand does not have a singularity. Hence where E accounts for the error and is a polynomial with interval coefficients. Typically, we will use as the points for the Taylor expansions x * (i) = α i since we will compare the resulting polynomial with another one of the form j b j (x − x i ) j and we will also choose . This choice is useful for two reasons: first, we will only have to integrate half of the terms since the rest will integrate to zero; and second, the error estimates will be better for this choice of y * (j) in the sense that the coefficients will be smaller. All the computations will be carried out using automatic differentiation. We should remark that we can get estimates for the error E in any of the above mentioned norms without having to recompute it since the relation holds for every q < k. Now, we move on to the term Hf N (x). In this case, we perform a Taylor expansion in both the denominator and the numerator where η belongs to an intermediate point between x and y, which we can enclose in the convex hull of [α i , α i+1 ] and [α j , α j+1 ] where the convex hull is understood in the torus.
Since typically K will be very small (compared to N ) there is no ambiguity in the definition. Finally, we can factor out (x − y) and divide both in the numerator and the denominator.
Since we know f (y) explicitly, we can perform the explicit integration and get a piecewise polynomial as a result.

Estimates of the norm of the Operator I + T
In this subsection we will outline how to compute the norm of the operator I + T = I + 2 BR(z, ·), z α . Since the operator T behaves like a Hilbert Transform plus smoothing terms, we will describe how to calculate rigorously with the help of a computer an estimate for the norm of its inverse. The procedure is more general and can be applied to a bigger family of kernels. Let T = R/2πZ, and let A(x), B(x) be real-valued functions on T. Also, let E(x, y) be a real-valued function on T × T. We assume A, B and E are given by explicit formulas such as as perhaps piecewise trigonometric polynomials or splines, and E(x, y) is a trigonometric polynomial on each rectangle I × J of some partition of T × T. We suppose A, B, E are smooth enough.
Let H be the Hilbert transform acting on functions on T, i.e.
Assume that A and B have no common zeros on T. Let Thus, S is a singular integral operator. We hope that S −1 exists and has a not-so-big norm on L 2 , but we don't know this yet. Our goal here is to find approximate solutions F of the equation SF = f for suitable given f ∈ L 2 (T), and to check that SF − f L 2 (T) < δ for suitable δ. Our computation of F will be based on heuristic ideas, but the computation of an upper bound for SF − f L 2 (T) will be rigorous. In our case, A(x) = 1, B(x) = 1.
To carry this out, let H 0 ⊂ H 1 ⊂ L 2 (T) be finite-dimensional subspaces, e.g. with H i consisting of the span of wavelets (from a wavelet bases) having lengthscale ≥ 2 −N i . Here N 1 ≥ N 0 + 3 (say). Let π i be the orthogonal projection from L 2 (T) to H i , and let us solve the equation If f is given explicitly in a wavelet bases, then (10) is a linear algebra problem, since π 1 Sπ 1 is of finite rank, and its matrix (in terms of some given basis for H 1 ) can be computed explicitly.
• If π 0 f ∈ Range(π 1 Sπ 1 ), then we find F ∈ H 1 such that π 1 Sπ 1 F = π 0 f , i.e. π 1 SF = π 0 f . We then have and both norms on the right-hand side may be estimated explicitly. Now, our goal is to make a heuristic computation of an operator of the form such that SS − I has small norm on L 2 (T).
Here, we will make a heuristic computation ofS; later we will give a rigorous upper bound for the norm of SS − I on L 2 (T). By a heuristic computation ofS we mean a heuristic computation ofÃ,B andẼ.
We first findÃ andB by setting Then, this means that So, from now on, we suppose thatÃ andB are known. For the operator I + T , this means A = 1/2,B = −1/2. We want to computeẼ. Now, let {φ ν } be some orthonormal basis for L 2 (T), for example a wavelet basis. By the previous methods, we can try to find functions ψ ν ∈ L 2 (T) such that Sψ ν − φ ν has small norm. We carry this for ν = 1, . . . , N for a large N . We now try to makeẼ satisfỹ Thus, we want Note that ψ # ν can be computed explicitly. Since the φ ν (all ν) form an orthonormal basis for L 2 (T), it is natural to definẽ This can be computed explicitly, and it satisfies (12). Thus, we can compute We claim that all terms enclosed in curly brackets are integral operators of the form for an E # that we can calculate. Let us go term by term Note that ifÃ is a piecewise trigonometric polynomial and C k , then E # can easily be computed modulo a small error in C k−1 .
• BHẼ has the form S # , with • EBH has the form S # , with This proves the claim. Letting dy be the operator in curly brackets in (13), we see that and that the function E # (x, y) can be computed modulo a small error in C 0 (T×T). Therefore, we obtain an upper bound for the norm of SS − I, namely Defining S err := SS − I, we obtain an explicit upper bound δ for the norm of S err on L 2 (T). We hope that δ < 1. If not, then we fail.
Suppose δ < 1. Then so we obtain a right inverse for S, namelyS(I + S err ) −1 , which has norm at most where S denotes the norm ofS as an operator on L 2 (T). Recall Therefore, Plugging that bound into (14), we obtain an explicit upper bound for the norm on L 2 of a right inverse for S. Similarly (by looking atSS instead of SS), we obtain an upper bound for the norm on L 2 of a left inverse for S.

Remark 2.1
To estimate e.g. max x T |E # (x, y)|dy it may be enough just to use the trivial estimate where (for each t),A(·, t), B(·, t), E(·, ·, t) are as assumed above. If A, B, E depend in a reasonable way on t, then one shows easilly that We can make η small by taking t 1 close enough to t 0 . Suppose we prove that S −1 t 0 ≤ C 0 by the previous methods. Then, of course we obtain an upper bound for 3 Bounds for C(t) and k 3.1 Writing the differential inequality as a differential system of equations The calculation of a bound for C(t) requires more effort than the previous one since one needs to calculate the terms one by one and add all their contributions to C(t). For example, in order to calculate the evolution of the norm D H k (t) a systematic approach is to take k derivatives (k ranging from 0 to 4) in the equation for the evolution of z (7 with f = g = 0), take another k derivatives in the equation for x (7 with arbitrary f, g) and subtract them. Let us focus from now on in the term Q(z) 2 BR(z, ω) − Q(x) 2 BR(x, γ) and its derivatives. One notices that in order to write a term in the variables (z, ω, ϕ) composed of a factors minus its counterpart in the variables (x, γ, ψ) in a suitable way (i.e. as a sum of terms that only have factors x, γ, ψ, D, d, D) then the number of terms is 2 a − 1. The way of writing it is the classical way of adding and subtracting the same term with the purpose of creating differences of terms and eliminate all the occurrences of the variables (z, ω, ϕ). An example for the Birkhoff-Rott operator (with Q = 1) is given next. We should remark that the computation and bounding of the Birkhoff-Rott is the most expensive one, the rest of the terms being easier.
After having seen this, it is clear that a tool that can perform symbolic calculations (derivation and basic arithmetic at least) and the correct grouping of the factors is required since the performance at this task by a human is not satisfactory. We developed a tool in 900 lines of C++ code that could do all this and output the collection of terms in Tex. We show an excerpt of the terms concerning the fourth derivative of BR(z, ω) − BR(x, γ). The total number of terms in that case is 2841.
However, there is a significant way to reduce the number of terms in the estimates: writing the equation in complex form instead of vector form. Thus, we can write the evolution for z in the following way: In this formulation, the fourth derivative accounts for only 140 terms. We present the first 10 below.
The final observation is that if we consider E(t) as a scalar, we might not get suitable estimates. In order to get better estimates, we will modify the energy into a "vectorized" version E v (t), which we will also denote by E(t) by abuse of notation. This new vectorized energy will be as follows where the homogeneous spacesḢ k have their norm defined by f Ḣk = ∂ k α f L 2 . With this vectorized system, we avoid both the bounding of any given norm by the full energy and any constant factor arising from interpolation between two Sobolev spaces. Thus, our constant C(t) will roughly be of a size comparable to the largest eigenvalue of the linearized system.

Estimates for the linear terms with Q = 1
Since we expect E(t) to be small, the terms that affect more to the evolution of E(t) are the linear ones. We now report on the non-rigorous experiments over the linear terms to obtain an approximate bound of the behavior of the full system (i.e. an approximation to the largest eigenvalue of the linearized system). We remark that a multiplication of the estimates by a constant, even a small factor 2 for example, has a big impact on the system, rendering the estimates useless and the estimations not tight enough, because the type of estimates we are going to get are exponential in the product of the time elapsed between the splash and the graph and the constant. Therefore, we should be very careful and fine estimates have to be developed.
First of all, we will work with Q = 1 and later move on to the case Q = 1. We will adopt the following convention to denote the different Kernels (integral operators) that appear: The operators for which b 2 = −1 will act on D or its derivatives whereas the operators for which b 2 = −1 will act on d or its derivatives. We now describe how to split the Kernels in such a way that they can be computed. For the case where b 2 = −1 we illustrate this by splitting Θ 0,0,0,0 2,0 , but the technique can be applied to any Kernel. where .
We can think of c 1 (α) and c 2 (α) as the Taylor coefficients of Θ(α, β) around β = α. We can bound the terms in (15) in the following way: We have then the estimates We now move on to T 1 . We will estimate it in the following way: To estimate the kernel T 2 we will use the Generalized Young's inequality [7]: we have that We finally show how to estimate the Kernels with b 2 = −1. We will do this by showing how to estimate Θ 0,0,0,0 1,−1 but the technique can be applied to any Kernel.
We can easily estimate these two terms applying to T 1 the same estimates (Young's inequality) as for T 2 in the previous case and by noting that T 2 is 1 2 c 1 (α)H(d).

Estimates for the linear terms with Q = 1
To perform the real estimates, where Q = 1 we will use the estimates from the previous sections. We will explain how to pass from the former ones to the latter ones. We will illustrate this by computing the linear terms of the Birkhoff-Rott operator.
First of all, the total number of terms will increase by a factor 2, since we will have In order to calculate the old terms with Q = 1, the only thing we have to do is to incorporate a factor of ∂ k α Q 2 (x)(α) in the estimates. The new terms can easily be calculated using that, up to linear order

Proof of Theorem 1.2
In this section, we will prove the stability Theorem 1.2. The equations are: f will be the error for z and g will be the error for ω.

4.1
Computing the difference z − x and ω − γ We define now: The energy and the Rayleigh-Taylor condition where ) and 2 (t)) 2 , k big enough depend on the norms of f and g.
Remark 4.1 From now on, we will denote E(t) + E(t) k by P (E(t)).
is left to the reader. We compute The first integral is easy to bound by CP (E(t)), we proceed as in the local existence Theorem I.7 in [3]. We split We have: Thus, we are done with I 3 . We now split I 1 = l.o.t + I 1,1 + I 1,2 + I 1,3 + I 1,4 where l.o.t stands for low order terms, nice terms easier to deal with.
In I 1 1,3,1 we find a commutator, which can be handled as before. It is also easy to estimate I 2 1,3,1 .
We show how to deal with I 4,2 1,4,1 . We compute we consider the most singular terms x α |x α | J 5 will be given later. In J 1 and J 2 we find 4th order terms in derivatives in z and x so they are fine. In J 3 we find inside the integrals This implies that we find "Hilbert" transforms applied to four derivatives of x and z. We are done with J 3 .
In J 4 we also find them inside the integrals (16) and (17) so it is easy to check that we have kernels whose main singularity is homogenous of degree 0 applied to four derivatives of ∂ 4 α ω and ∂ 4 α γ. This implies that we have a Hilbert transform applied to ∂ 3 α ω and ∂ 3 α γ so we are done with J 4 . The most dangerous term is J 5 which is given by We split further In J 5,1 we find a Hilbert transform applied to ∂ 4 α z ⊥ and ∂ 4 α x ⊥ so it is fine. We split further: J 5,2,1 can be estimated as before (there are more derivatives: 5 in total, but they are in x). In J 5,2,2 we find a commutator. Finally: We use that H(Λ) = −∂ α and z α · ∂ 4 α D ⊥ = −z ⊥ α · ∂ 4 α D to obtain: Then we are done with I 4,2 1,4,1 , I 4 1,4,1 , I 1,4,1 , I 1,4 and I 1 . To finish with I it remains to control I 2 . We split it as: The low order terms are easier to deal with. We further split I 2,1 .

Error term
We deal with I 2,2,1 more carefully. We use that We can integrate by parts in all of the above terms to get low order terms. We are finally done with I.

Computing the difference ϕ − ψ
From the local existence proof we find the equation for ϕ t : We will show how to find the equation for ψ t . We start from The equation for γ t reads: Then We should remark that we have used that Computing, we find that are error terms. We consider It is easy to check that With this formula it is easy to find that In order to deal with II we take a derivative in α in the equation for ω and ψ to reorganize the most dangerous terms. If we find a term of low order, we will denote it by NICE. Since the equations for ϕ t and ψ t are analogous except for the E 1 term, the NICE terms are going to be easier to estimate in terms of CP (E(t)) + cδ(t).
Expanding (3): We use that The term (|x α |B x (t)) t depends only on t so it is not going to appear in computing II.
2|x α | is a transparent term which is NICE (even if we have to deal with Λ 1/2 ) The first term is at the level of ∂ α x so it is NICE. The second term is at the level of ∂ α x or ψ α so it is NICE. We write the last one as The first term is at the level of x α or ψ so it is NICE. For the second term we have used that Finally: The first term is at the level of x α , x t , BR ∼ x α so it is NICE. We use that Using that That yields NICE (at the level of xα,xt,BR) Finally: which means We use (19) to deal with (Qx)αt |xα| . We find that The fact that the last two terms are NICE, allows us to find that Finally: We gather all the formulas from (4) to (12) absorbing the error terms byẼ 1 α whenever we encounter them.
We add and subtract terms in order to find the R-T condition. We remember here that In σ x there are error terms but they are not dangerous. Then, we find Line (19) can be written as We expand x αt to find We claim that NICE, we use that |xα| 2 only depends on time Hilbert transform applied to γα This means that We write whereÊ is an error term. To simplify we write Setting the above formula in the expression of (20)+(21) allows us to find being E 3 α a new error term. We now complete the formula for σ x in (20) to find x αα · x ⊥ α |x α | 3 + errors Thus (22) + (23) + (24) = Q 3 1 |x α | 3 D x (α) + errors = NICE + errors Finally, we obtain ψ αt = NICE(x, γ, ψ) − Q 2 x σ x x αα · x ⊥ α |x α | 3 + E 4 α II 2,3,3 ≤ CP (E(t)), l.o.t In II 2,3,2 we split further: II 2,3,2 = II 1 2,3,2 + II 2 2,3,2 Then For II 2 2,3,2 : It remains Then II 1 2,3,1 ≤ CP (E(t)) + cδ(t) At this point we remember that we had to deal with so in II 1 2,3,1 we find one derivative less (or 1/2 derivatives less) and this shows that we can bound II 1 2,3,1 ≤ CP (E(t)) + cδ(t) by brute force. It remains to show claim (22). We remember α this term is also in H 3 We write

α) + errors = NICE + errors
Finally, the most singular terms in Q 2 x σ x are We take 3 derivatives and consider the most dangerous characters: In M 2 we find