Remarks on geometric properties of SQG sharp fronts and $\alpha$-patches

Guided by numerical simulations, we present the proof of two results concerning the behaviour of SQG sharp fronts and $\alpha$-patches. We establish that ellipses are not rotational solutions and we prove that initially convex interfaces may lose this property in finite time.


Introduction
Our goal in this article is to investigate some simple questions about the evolution of the surface quasi-geostrophic (SQG) front, some of which are well known for the vortex patch problem. The surface quasi-geostrophic equation was introduced in the mathematical community by Constantin, Majda and Tabak in [5]. This equation is derived considering small Rossby and Ekman numbers and constant potential vorticity. It provides a mathematical description of the evolution of the temperature from a general quasi-geostrophic system for atmospheric and oceanic flows (see [20] for more details).
Leaving aside its interest from a geophysical point of view, the SQG equation is a two dimensional model of the principal equations of inviscid fluid dynamics in 3D, namely, the incompressible Euler equations. The SQG system is the following set of equations: where Λ = (−∆) 1 2 and Ψ is the stream function. The first equation simply represents the fact that the temperature, θ, is advected by the velocity u. The second and third equations relate the temperature, θ, with the velocity, u, through a non local operator. Those two equations can be rewritten as u = (−R 2 θ, R 1 θ) where R i are the Riesz Transforms, given 1 arXiv:1401.5376v1 [math.AP] 21 Jan 2014 Motivated by the articles [5] and [12], a lot of effort has been devoted to understanding these equations. In particular, the problem of whether the SQG system presents finite time singularities or not is open.
The existence of global weak solutions in L 2 was proven by Resnick in [21] using an extra cancellation due to the oddness of the Riesz transform. A particular kind of weak solution for an active scalar are the so called patches, i.e., solutions for which θ is a step function: where Ω(0) is given by the initial distribution of θ, and Ω(t) is the evolution of Ω(0) under the velocity field u given by u( The evolution of such distribution is completely determined by the evolution of the boundary, allowing the problem to be treated as a non-local one dimensional equation for the contour of Ω(t). In this setting, local existence of smooth solutions was first obtained for analytic curves by Rodrigo in [22]. The question of local existence of simply connected patches with Sobolev regularity of its boundary was addressed by Gancedo in [9].
To get a better understanding of the behaviour of solutions of these interface problems, several numerical experiments have been performed. In [6] the problem of the evolution of two patches is studied. Their simulations suggest an asymptotically self-similar singular scenario in which the distance between both patches goes to zero in finite time while simultaneously the curvature of the boundaries blows up. Recently, Gancedo and Strain [10] proved that in fact, no splash singularity can be formed, i.e., two interfaces can not collapse in a point, if the interfaces remain smooth.
Also in [6], a new family of patch problems interpolating between the vortex patch and the SQG patch is introduced: the α-patch model. For the vortex patch problem, the velocity field is obtained as v = ∇ ⊥ ∆ −1 θ, where θ is the vorticity, while for SQG the velocity is given by . The local existence proof for α ∈ (1, 2) can be found in [3].
The evolution equation for the interface of an α− patch, which we parametrize as a 2π periodic curve z(x), can be written as if 0 < α < 2, and as for the vortex patch, due to the fact that the evolution of the boundary of the patch is unchanged if the velocity is modified by adding a tangential term to the boundary, since that amounts just to a change of parametrization of the interface. For more details see [9], [22]. We make use of this property in our numerical simulation, adding a tangential term in order to keep the modulus of the tangent vector constant in space.
In [23], based on numerical simulations, it is suggested that an elliptical patch with a big aspect ratio between its axes may develop a self-similar singularity with an explosive growth of the curvature. In [24], it was already pointed out that small perturbations of thin strips may lead to a self similar cascade of instabilities, leading to a possible arc chord blow up.
Beyond that, very little is known about the qualitative behaviour of patch like solutions of the SQG equation.
The analogous problem for the vorticity formulation of 2D Euler (α = 0) is better understood. The global existence and uniqueness of weak solutions of the 2D Euler in vorticity formulation is due to Yudovich [27]. Regularity preservation for C 1,α patches was obtained by Chemin using techniques from paradifferential calculus in [4]. Another proof of that result, which highlights the extra cancellation on semi spheres of even kernels, can be found in [1] by Bertozzi and Constantin. In recent years, Denisov has studied the process of merging for the vortex patch problem. This is the scenario showed by the numerics of [6] for the alphapatch. However, for the vortex patch problem, the collapse in a point can not happen in finite time, the distance between the two patches can decay at most as fast as a double exponential. Denisov proves, in [8], that this bound is sharp if one is allowed to modify slightly the velocity by superimposing a smooth background incompressible flow.
In addition, it is known that there exist several solutions which evolve by rotating with constant angular velocity around its center of mass. The ellipses are the most well known example of solutions exhibiting this behaviour, and in fact, the only simply connected ones for which there are explicit formulae. The existence of that kind of solutions can be proven by a bifurcation analysis from the circular solution. For more details, see [13] and references therein.
Nowadays, the bigger computation capacity of computers has lead to their use as a mathematical tool. However, floating-point operations can result in numerical errors. To deal with this matter 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 x ∈ X, y ∈ Y x y ∈ X Y, for any operation . For example, We can also define the interval version of a function f (X) as an interval I that satisfies that for every x ∈ X, f (x) ∈ I. Rigorous computation of integrals has been theoretically developed since the seminal work of Moore and many others [2,16,18,19,25], and has had applications in physics and fluid mechanics [15,11]. In order to perform them we used the C-XSC library [14].
Guided by new numeric simulations of the evolution of elliptical patches, we establish the following two results for the general α-patch: • The ellipses are not rotating patches for the equation (3) with 0 < α < 2.
• For every 0 < α < 2, there exists a solution of the α-patch equation that begins convex, and after a finite time it is no longer convex. Moreover, there exists a solution of the vortex patch problem that has the same property.
The proof of those two results relies on establishing a sign for certain integrals. For the first theorem, we can do it by performing certain manipulations in the integral. However, for the second theorem we recur to a computer assisted proof.
The article is organized as follows: In the second section we discuss our numerical experiments. Section 3 presents a simple proof of the fact that ellipses are not rotating solutions of the α-patch problem for α > 0. This is followed by a computer assisted proof of the loss of convexity in Section 4.

Numerical Simulations
In order to obtain insight on how the ellipses evolve under SQG dynamics we performed simulations for elliptical SQG patches with varying aspect ratio. From these simulations three main differences with the behaviour under vortex patch evolution are easily appreciated: • The ellipses are not rotating solutions.
• Initially convex patches may lose that property.
For ellipses with an aspect ratio close to 1, the evolution is similar to the one of the vortex patch problem, since all points at the boundary behave in similar way. However, as the aspect ratio is increased, the fact that the ellipses do not rotate and lose convexity becomes clearer. Figures 1-3 represent four timestamps for the evolution of the ellipses with principal semiaxes 1-3, 1-5 and 1-6. The patch with initial data z(x) = (cos(x), 3 sin(x)) displays a combination of a rotating motion with a smaller scale oscillation which leads to loss of convexity.
When the bigger semiaxis is increased to 5, z(x) = (cos(x), 5 sin(x)), the patch exhibits an almost periodic behaviour in which the arc chord, i.e., the distance between different points at the boundary, decreases but eventually increases again. In addition, the loss of convexity is evident. See Figure 2.
For the ellipses with a bigger difference between its principal semiaxes, the simulations suggest that the arc chord goes to zero in finite time, developing a singularity. This question has already been addressed in [23] and [24], where more precise numerical simulations are discussed which suggest a self similar behaviour. We remark that our simulations are here used as a tool to gain intuition, not as a final purpose themselves and are by no means of a comparable accuracy to those in [23,24].
The numerical simulations were performed using a contour dynamics algorithm. The interface is discretized by a 512 point grid and the interface evolves in time with a tangential component such that it maintains |∂ x z| constant in space. In addition, derivatives are computed by a Fourier method with an exponential filter. Evolution in time is carried over by an adaptive step Runge-Kutta 45 pair.

Rotating solutions
The purpose of this section is to show some differences between the rotating solutions of the vortex patch and the α-patch problem.
It is well known that the circular patch is a steady solution of all those problems, as can be easily seen by noticing that the operator (−∆) 1− α 2 is a radial multiplier, and hence, it maps radial functions into radial functions. Since the velocity field is the skew gradient of (−∆) 1− α 2 θ, where θ denotes the step function of the patch, the velocity field is tangent to the level sets of θ, leaving them invariant.
In addition to the circular patch, for the vortex patch problem, there is a family of rotating solutions, the so called V-states. These are patches which only rotate around their center of mass with constant angular velocity and which have n-fold symmetry. The other remarkable property of V-states is that they have n-axes of symmetry. Some numerical results about the shape of V-states can be found in [7] and [26] and references therein. The existence of V-states can be proven by a bifurcation analysis from the circle. However, there is no explicit formula for the V-states with more than two axes of symmetry, which correspond to ellipses. The fact that the ellipses are rotating solutions for the vortex patch problem is a classical result [17]. However, we will show that ellipses are not rotating α-patches, for α > 0. Proof: Assuming that the center of mass of the patch is the origin, the condition that the patch rotates with constant angular velocity is equivalent to the existence of an Ω such that: for all x ∈ [0, 2π], where z(x, t) is a complex parametrization of the boundary. Since the evolution of the patch is completely determined by the normal component of the velocity at the boundary, it is enough to check that the previous identity is not satisfied in the normal direction.
Plugging z(x, t) = (R 1 cos(x), R 2 sin(x)) into the integral formula for the velocity and taking the inner product of both sides of the previous equation with ∂ x z(x) ⊥ , we are led, after some lengthy manipulations, to check that , 0 < R < 1, and are not proportional.
We will do so by showing that H(x) = F (x) G(x) has different limits at x = 0 and x = π 2 . F vanishes both at x = 0 and x = π 2 due to the fact that the integrand is odd with respect to those points, and hence, H(x) is continuous at x = 0 and x = π 2 with limits We will check that H(0) − H( π 2 ) = 1 2 (F (0) + F ( π 2 )) > 0.
Taking derivatives with respect to x and evaluating at x = 0 and x = π 2 , we obtain where the constant C α,R might change from line to line. Finally, the last integral is strictly positive since cos(x) > sin(x) > 0 for 0 ≤ x ≤ π 4 .

Loss of convexity
This section is devoted to prove the following theorem: We start with some technical lemmas: Then, the following inequalities are true: where = 0.0078125 = 1 128 .
Proof: The proof is computer-assisted and the codes can be found in the supplementary material. To compute the positiveness or negativeness, the program checks recursively halving the interval that is validating, until it either achieves a sign condition or the interval is less than a minimum size (in our case 2·10 −10 ). The program passed all the checks without problems.
From the previous lemma follow the next corollaries: be defined on [−π, π] and let C = 0.45 or C = 0.15. Then the curvature is only zero at x = π.
Proof: The proof follows from the fact that the denominator of K(x) is never zero and the numerator is equal to k C (x)e Corollary 4.4 Let z 1 (x), z 2 (x) be defined as in the previous corollary and let = 1 128 . Then, for every x ∈ [−π, −π + ], ∂ k x z 1 (x) belongs to the convex hull of ∂ k x z 1 (−π) and ∂ k x z 1 (−π + ) and for every x ∈ [π − , π], ∂ k x z 1 (x) belongs to the convex hull of ∂ k x z 1 (π − ) and ∂ k x z 1 (π) when k = 0, 1, 2, 3, 4, 5.
Proof: The proof is immediate after noticing that Proof of the Theorem: The proof of the theorem is computer-assisted and can be found in the supplementary material. Let us explain the necessary reductions of the problem to a computationally tractable one.
The general strategy is the following: we will start with an initial condition such that the curvature is zero at exactly one point x and positive otherwise. Then, we will rigorously prove that the derivative of the curvature at that point is different than zero. Since the α-patch equations are time reversible, we can go a sufficiently small time forward and backward so that we are able to show that there exists a solution that starts being convex and loses the convexity after a finite time by choosing the initial jump of temperatures across the boundary to be either +2π or −2π.
Remark 4.5 By the time reversibility of the equations, this also shows that there are initial conditions that are non convex and in finite time they become convex.
More precisely, we will show the following: if z 1 and z 2 are defined as in Corollary 4.3, then: [0.09, 2.00).
We now detail the formulas needed for the computation of the time derivative of the curvature.
where the last term is zero since the curvature is zero at the point at which we are evaluating the expression. For the vortex patch equation, taking derivatives from the contour equation (4) and integrating by parts yield c α < ∞ and that I(0) = 0, it is enough to show that DI(α) = ∂ α I(α) < 0. This is the sign condition that we will validate in the second region. In the third and fourth regions, we will validate I(α) < 0 or I(α) > 0 depending on the range of α. We chose α cr = 0.04 and α br = 1.95. We remark that both I(α) and DI(α) have the same sign if multiplied by |z x | 3 . Therefore, in order to avoid the division, we will compute the values of I(α)|z x | 3 and DI(α)|z x | 3 . The integral that we need to calculate in the small α region is • Two intervals, Left and Right, which set the limits for the bounded and singularity regions (i.e. singularity = [Left, Right], bounded = [−π, π]\ singularity). In our proof, Left = − 1 128 , Right = 1 128 . See below for the precise definition and the integration procedures in these regions.
• Two doubles, AbsTol and RelTol, which limit the precision up to which the integrals are computed. In our proof, AbsTol = RelTol = 10 −6 .
• One interval, α, which is the interval in the parameter space we are calculating.
• One interval, C, which corresponds to the different initial conditions.
We will use a queue (implemented using the Standard Template Library (STL) Queue), in which we store all the ParameterSets to be computed. While the queue is not empty, we take the top element, pop it and give an enclosure of integral we are validating for this region. Three different things can happen: • The enclosure is positive.
• The enclosure is negative.
• We can not say anything about its positivity.
In the first two cases, the result is output to its corresponding file (one for the regions for which the aforementioned integral is positive, another for the ones for which it is negative). In the third case, the ParameterSet is split into other narrower ParameterSets which are pushed in the queue. This splitting is only done if the diameter of α is bigger than a given threshold, which in our case was set to 5 · 10 −6 .
We now describe how to perform the integral over the bounded region. The bounded part is calculated using a Gauss-Legendre quadrature of order 2, given by Whenever the result does not satisfy some tolerance requirements in the form of having absolute or relative (with respect to the volume of the integration region) width smaller than the two constants AbsTol and RelTol the integration domain is split by the midpoint and we call the integrator again with the new two subdomains recursively. Otherwise we will save it and add it to the total. We also limit the depth of the levels of splitting in order to prevent infinite loops or stack overflows because of too stringent tolerances since the uncertainty of the parameters might yield wide enclosures of the integral even with infinite precision. In our case, the maximum number of subdivisions was 13, totalling a maximum number of subintervals equal to 2 13 .
Regarding the singularity region, we remark that if we try to evaluate the (singular) integrand directly, we get expressions of the type 0 0 which can only be enclosed in unbounded intervals. In order to avoid this phenomenon we will bound the absolute value of the integrand and consider the contribution of the singularity region as a residual error. By virtue of the mean value theorem, we will substitute any expressions of the form |(∂ k x z(a) − ∂ k x z(b))| by the interval |(a−b)||∂ k+1 x z([a, b])|, and these derivatives can be easily enclosed only by having to evaluate at the endpoints because of the monotonicity proved in Corollary 4.4. After estimating the integrals in this way, we are left with the task of having to bound quantities of the form C|x| 0 dx, C|x| 1−α dx, C 1 |x| 1−α dx + C 2 | log(C 3 x)||x| 1−α dx, or C|x| 2−α dx depending if we are in the first, second, third or fourth regime of α respectively, for some constants (intervals) C, C 1 , C 2 , C 3 . These integrals can be computed explicitly in closed form. In practice, due to the exponential nature of z 1 , the contribution of these integrals is extremely small compared to the true result. We ran the computation in parallel (every core was allocated a different initial region) over 16 cores by splitting the regions of α on an Intel i5 processor with 4 GB of RAM. The total combined computation time was roughly 120 hours.