Global well-posedness of the free-interface incompressible Euler equations with damping

We prove the global well-posedness of the free interface problem for the two-phase incompressible Euler Equations with damping for the small initial data, where the effect of surface tension is included on the free surfaces. Moreover, the solution decays exponentially to the equilibrium.

For each t > 0, the two fluids are described by their velocity and pressure functions, which are given by u ± (t, ·) : Ω ± (t) → R 3 and p ± (t, ·) : Ω ± (t) → R, respectively. For each t > 0, we require that (u ± , p ± , h ± ) satisfy the following free boundary problem for the two-phase incompressible Euler equations with damping, with taking into account of the effect of surface tension on the free surfaces, (1. 3) The third and fourth equations in (1.3) state that the free surfaces Σ ± (t) move with the velocities of the two fluids; the fourth equation implies in particular that the normal components of the velocities are continuous across the free interface Σ − (t).
The first two equations in (1.3) is the incompressible Euler equations with damping (or dissipation), where a ± > 0 is the damping coefficient; the equations is used in geophysical models for large-scale processes in atmosphere and ocean [27,14], where the damping term models the friction due to the bottom of the ocean or the Rayleigh friction (or the Ekman pumping/dissipation) in the planetary boundary layer in the presence of rough boundaries [34,6,16]. We may refer to, for instance, [3,29,10,22,9] for some mathematical results of the incompressible damped Euler equations. p atm is the constant pressure of the atmosphere and σ ± > 0 are the surface tension coefficients. Finally, N ± = (−∂ 1 h ± , −∂ 2 h ± , 1) are the upward (nonunit) normal vectors to Σ ± (t), and H ± are twice the mean curvature of the free surfaces Σ ± (t): To complete the statement of the problem, we must specify the initial conditions. We suppose that the initial two free surfaces are given by the graphs of the functions h ± (0) = h 0,± : T 2 → R, which yield the initial domain Ω ± (0) on which the initial velocities u ± (0) = u 0,± : Ω ± (0) → R 3 are specified. We will assume that + h 0,+ > h 0,− > −b on T 2 .

1.2.
Background. The local well-posedness of the two-phase incompressible Euler equations with surface tension, without damping, has been proved by Cheng, Coutand and Shkoller [7]; see also Shatah and Zeng [30,31], Cheng, Coutand and Shkoller [8] and Pusateri [28]. We may also refer to Ambrose [1] and Ambrose and Masmoudi [2] for the local well-posedness of the irrotational flows. These results of the local well-posedness were established for the large initial data, which in particular implies the short-time structural stability of vortex sheets with surface tension. The presence of surface tension is necessary for the stability of vortex sheets in the two-phase incompressible Euler equations; without surface tension, we have the well-known Kelvin-Helmholtz instability, see Caflisch and Orellana [5] and Ebin [15]. For the short-time structural stability of vortex sheets in the two-phase compressible Euler equations, we refer to Coulombel and Secchi [11,12] and Stevens [33].
Up to our best knowledge, the global stability of vortex sheets remains open. With including damping, there are no trivial vortex-sheet solutions to the Euler equations, and hence there is no meaning to consider their stability for the Euler 1.3. Reformulation in flattening coordinates. In order to transform the free boundary problem (1.3) to be one in the fixed domain, we will use a flattening transformation as [36,24] rather than the Lagrangian coordinate transformation. To this end, we define the fixed domain for which we have written the coordinates as x ∈ Ω ± . We shall write Σ + := {x 3 = } for the upper boundary, Σ − := {x 3 = 0} for the internal interface and Σ b := {x 3 = −b} for the lower boundary. Throughout the paper we will write Ω = Ω + ∪ Ω − and Σ = Σ + ∪ Σ − . We think of h ± as functions on Σ ± according to h + : R + × (T 2 × { }) → R and h − : R + × (T 2 × {0}) → R, respectively. We will transform the free boundary problem in Ω ± (t) to one in the fixed domains Ω ± by using the unknown free surface functions h ± . Indeed, we will flatten the coordinate domains via the following special coordinate transformation: Here η =b 1 P + h + +b 2 P − h − , where P ± are the Poisson extensions defined by (A.3) and (A.8), respectively, andb Note that if h ± are sufficiently small in an appropriate Sobolev space, then ∂ 3 ϕ = 1 + ∂ 3 η > 0 and hence the mapping Φ is a diffeomorphism. This allows us to transform the free boundary problem to one on the fixed domains Ω ± . Since the domains Ω ± and the boundaries Σ ± are now fixed, we henceforth consolidate notation by writing f to refer to f ± . When we write an equation for f we assume that the equation holds with the subscripts added on the domains Ω ± or Σ ± , etc.; when we write an equation involving f on Σ − , then it means that the equation holds for both f = f + and f = f − ; when necessary to distinguish the two, we will write out the subscripts explicitly. To write the jump conditions on Σ − , for a quantity f = f ± , we define the interfacial jump as We assume that the initial surface functions satisfy the zero-average conditions Indeed, by the third, sixth and second equations in (1.11), we have The conditions (1.13) will allow us to apply Poincaré's inequality for h ± on Σ ± for all t ≥ 0. If it happens that the initial surface functions do not satisfy the zero average conditions (1.12), then we can shift the data and the coordinate system so that (1.12) is satisfied, provided that the extra condition + h 0,+ > h 0,− > −b is satisfied; see [36] for instance. The problem (1.11) possesses a natural physical energy. For sufficiently regular solutions, we have an energy evolution equation that expresses how the change in physical energy is related to the dissipation induced by the damping: . (1.18) The local well-posedness of the problem (1.11) with surface tension in our energy functional E n (n ≥ 3) can follow exactly in the same way as that of the problem without damping; indeed, the local well-posedness of the two-phase incompressible Euler equations with surface tension in the energy functional E 2 (i.e., n = 2) was proved by Cheng, Coutand and Shkoller [7], and one could then check step by step that there would not be any essential difficulties to generalize the local wellposedness result in [7] to the case n ≥ 3 or with damping. We may then record the local well-posedness of (1.11) as the following. Here P is a generic polynomial.
Our main results of this paper are stated as follows.
and the zero-average conditions (1.12). There exists a universal constant ε 0 > 0 such that if E n (0) ≤ ε 0 , then there exists a global unique solution (v, q, h) solving (1.11) on [0, ∞). Moreover, there exists universal constants C, γ > 0 such that for all t ≥ 0, and (1.21) Since h is such that the mapping Φ(t, ·), defined in (1.6), is a diffeomorphism for each t ≥ 0. As such, one may change coordinates to y ∈ Ω(t) to produce a global-in-time, exponentially decaying solution to (1.3).
Remark 1.2. We can also consider the two-phase Euler equations with different densities under the influence of the gravitational force. In this case, if the upper fluid is heavier than the lower fluid, then the fluids are susceptible to the well-known Rayleigh-Taylor gravitational instability. We may prove that there exists a critical surface tension value σ c such that if σ − > σ c , then the problem is nonlinearly stable; if 0 < σ − < σ c , the the problem is nonlinearly unstable. The critical value σ c is same as the viscous problem [35,36,18,24,25]. We will report this in a forthcoming paper.
We now present a sketch of main ideas in the proof of Theorem 1.2, and by Theorem 1.1 it suffices to derive the a priori estimates as recorded in Theorem 6.1.
The first step is to utilize the geometric structure of (1.11) and the energydissipation structure (1.16) to derive the following temporal energy evolution estimate:Ē The derivation of the estimate (1.22), which is in spirit of the local well-posedness theory in [7,34], relies heavily on the definition of E n (cf. (1.18)), especially that ∂ n−1 For the one-phase incompressible Euler equations [13], one can use either ∂ n−1 to close the highest order temporal energy estimate. However, for the two-phase incompressible Euler equations, especially in our flattening coordinates, ∂ n t v ∈ H −1/2 (Σ) is necessary to be included in E n . Indeed, compared to the one-phase incompressible Euler equations, the new and most technical difficulty is the control of the following boundary integral on Σ − : by (1.24) Here + · · · means a sum of terms that can be controlled easily by C(E n ) 3/2 . The first integral in the right hand side of (1.24) is controlled by the To estimate the second integral we use the third equation with v − on Σ − to write, upon an integration by parts in x * , It is now crucial to convert the first integral in the right hand side of (1.25) to be an integral over the interior domain Ω − , by using the divergence theorem, where ξ is a nonnegative cutoff function so that supp ξ ⊂ Ω and ξ = 1 on Σ − and we extend v + to be defined on Ω. The second integral in the right hand side of (1.26) is controlled by using the second equation of (1.11) in Ω − . For the first integral, we use the first equation of (1.11) in Ω − to write . The second integral in the right hand side of (1.27) is controlled by integrating by parts in x * . Now the situation is that it remains to estimate the first integral in the right hand side of (1.27) and the second integral in the right hand side of (1.25). However, there are both out of control and can not be estimated by using further the equations. To overcome this difficulty, we will introduce the differential operator The key point lies in that when considering D − t ∂ n−1 t instead of ∂ n t , those two troublesome integrals are canceled out; these allow us to conclude the estimate (1.22).
The next goal is to replaceĒ n andD n in the left hand side of (1.22) by the full one E n . To derive the estimates for the normal derivatives of v, the natural way is to estimate instead the vorticity curl ϕ v = ∇ ϕ × v to get rid of the pressure term ∇ ϕ q, and the conclusion is that We also use the equations of curl ϕ v to derive the bound for the (H 1 (Ω)) -norm of curl ∂ n t v. The similar estimates for div v follow easily by using the incompressible condition. The bound of the (H 1 (Ω)) -norm of curl ∂ n t v and div ∂ n t v yield in particular that ∂ n t v ∈ H −1/2 (Σ). Note that to utilize the vorticity and divergence estimates so as to employ the Hodge-type elliptic estimates of ∂ j t v for j = 0, . . . , n − 1, we need to show that This will follow from the third equation in (1.11) by showing that ∂ j t h ∈ H 3 2 (n−j)+1 (Σ) for j = 1, . . . , n. By chaining the Hodge-type elliptic estimates of v, the regularizing elliptic estimates of h by the presence of surface tension and the derivation of the estimates of q, we can show an iteration- )+1 (Σ), see Proposition 5.1 for more details. On the other hand, by the third equations in (1.11), we have (1.29) So to conclude the estimates, it remains to show that ∂ n t h ∈ H 1 (Σ). Since |∂ n t h| 2 1 ≤ E n , this implies that for the estimates of E n in the energy, we have However, for the estimates of E n in the dissipation, the difficulty is that |∂ n t h| 2 1 is not controlled byD n . The most crucial point of our global analysis is the observation of the following: JIALI LIAN See (5.24) for the proof of (1.31). Hence, we have the following time-integrated dissipation estimate We remark that the L 2 norm of h is controlled by using Poincaré's inequality since Consequently, combining (1.22), (1.28), (1.30) and (1.32), we can close the a priori estimates for E n is small and then the decay estimate can be derived also, as recorded in Theorem 6.1. Hence, the proof of Theorem 1.2 is completed.
1.5. Notation. We use the Einstein convention of summing over repeated indices. Throughout the paper C > 0 denotes a generic constant that does not depend on the data and time, but can depend on the the parameters of the problem, e.g., a, b, σ and n. We refer to such constants as "universal". Such constants are allowed to change from line to line. We employ A 1 A 2 to mean that A 1 ≤ CA 2 for a universal constant C > 0.
We write N = {0, 1, 2, . . . } for the collection of non-negative integers. When using space-time differential multi-indices, we write N 1+d = {α = (α 0 , α 1 , . . . , α d )} to emphasize that the 0−index term is related to temporal derivatives. For α ∈ N 1+d , For any differential operator ∂ α , we use the commutators We denote x * = (x 1 , x 2 ) for the horizontal coordinates and v * = (v 1 , v 2 ) for the horizontal components of v. We denote ∇ * for the horizontal gradient, div * for the horizontal divergence and ∆ * for the horizontal Laplacian. For simplify the notations, we still use ∇f to mean ∇ * f , etc., when without confusion, for function f defined on Σ. We omit the differential elements dx and dx * of the integrals over Ω and Σ, and also sometimes the time differential elements.

2.
Preliminary. In this section we record some preliminary results that will be used in the derivation of the a priori estimates for the solutions to (1.11). We will assume throughout the rest of the paper that the solutions are given on the interval [0, T ] and obey the a priori assumption for an integer n ≥ 3 and a sufficiently small constant δ > 0. This implies in particular that (2.1) and (2.2) will be used frequently, without mentioning explicitly.
In order to derive the estimates for the time derivatives of the solutions to (1.11) (essentially for the highest order), as for the local well-posedness, it is natural to utilize the geometric structure of the equations given in (1.11); see also Guo and Tice [19,20,21] of utilizing this structure to handle with the pressure term. We first apply ∂ j t for j = 0, . . . , n − 1 to (1.11) to find that where Furthermore, we may write that for j ≥ 1, and (2.8) It thus follows that As explained in the introduction we can not simply use the equations (2.3) with taking j = n to derive the estimates for the n-th order time derivative of the solutions, and we introduce the following differential operator where V − is the extension of v − , defined in Ω − , onto Ω, with taking the value of v + near Σ + . Note that D − t is the material derivative operator in Ω − but not in Ω + ; the important advantage of introducing V − not rather using v is that it is continuous across the interface Σ − . Note also that In particular, D − t is a tangential derivative on ∂Ω ± and is continuous across Σ − . We now restate the two key points, mentioned in the introduction, of using∂ n t instead of ∂ n t . First, it follows from the first equation of (2.3) in Ω − with j = n − 1 that∂ (2.13) Second, it follows from the third equation of (2.3) with j = n − 1 that (2.16) Now applying D − t to the equations (2.3) with j = n − 1, we find that, by (2.16), where Moreover, applying D − t to the equation (2.9) with j = n − 1, we havẽ ∂ n t H = ∇ · ∇∂ n t h 1 + |∇h| 2 − ∇h · ∇∂ n t h 1 + |∇h| 2 3 ∇h + F 4,n , We present the estimates of these nonlinear terms F i,n .
Lemma 2.1. It holds that Proof. We expand these commutators in F i,j into a sum of products and then control each product with the highest order derivative term in H r (r = 0, 1/2, 1, accordingly) and the lower order derivative terms in H m for m depending on r, using Lemmas A.2 and A.1, the trace theory along with the definition of E n (cf. (1.18)); all of them are bounded by E n . We remark that it is needed to have included the term ∂ n+1 t h 2 −1/2 in the definition of E n so that when estimating F 1,n , by Lemma A.1, The estimate (2.23) follows by summing.
In order to utilize the linear structure of (1.11), it is convenient to write the equations in the linear perturbed form where Here we have used the fact that, by (1.10), ∂ ϕ i − ∂ i = −∂ i η∂ ϕ 3 for i = t, 1, 2, 3. We present the estimates of these nonlinear terms G i .

Lemma 2.2. It holds that
Proof. We first apply the spatial-time mixed differential operators to G i to expand them as sums of products and then proceed in the same way as Lemma 2.1, expect when estimating ∂ n t G 4 we need to use first the structure of G 4 : The estimate (2.30) follows by summing.
3. Temporal estimates. In this section we will derive the energy evolution estimates for the temporal derivatives of the solutions to (1.11). For a generic integer n ≥ 3, we define the temporal energy bȳ and the corresponding dissipation bȳ To derive the estimates for the highest order time derivative of the solutions, we shall use the equations (2.17). A key point is that the various norms of the error between∂ n t and ∂ n t of the solutions are bounded by CE n . We record the following time-integrated temporal energy evolution estimate.
Proof. We only present the estimates of the highest order time derivative of the solutions, and the lower orders can be estimated in a much simpler way. Taking the inner product of the first equation in (2.17) with∂ n t v and then integrating by parts over Ω ± by using the second, third and sixth equations in (1.11), we obtain By (2.23), we have By integrating by parts over Ω ± , using the sixth and second equations in (2.17), we obtain By (2.23), we obtain (3.10) Integrating by parts in both x * and t, we find that (3.12) Upon an integration by parts in x * , we find that Hence, (3.7)-(3.14) implies Similarly, we have Now, we estimate I 3 . By integrating by parts in t and x * , we obtain Here we have used the H 1/2 (Σ − )-H −1/2 (Σ − ) dual paring, the trace theory, Lemma A.3 and the fact that |∂ n t v| 2 −1/2 E n . Hereafter, for the conciseness of presentations, we omit to write out explicitly the estimates of many terms which are bounded by C(E n ) 3/2 ; one can check them easily from the context. Since by the third equation of (1.11), v · N = 0 on Σ − , we have Then we obtain, by a further integration by parts in x * , Hence, we obtain, by (2.23), Now we convert the remaining integral in the right hand side of (3.22) to be one over the interior domain Ω − . For this, let ξ be a nonnegative cutoff function so that supp ξ ⊂ Ω and ξ = 1 on Σ − . Let V + be an extension of v + , defined in Ω + , onto Ω − . By using the divergence theorem, we have By using the second equation of (2.17) and (2.23), we obtain By integrating by parts in x * and using the first equation of (2.17) and (2.23), we have By integrating by parts over x * , we obtain Finally, we turn to estimate I 4 . By integrating by parts in t, we obtain To estimate the second integral in (3.28), since ∇ ϕ · v = div * v * + 1 ∂3ϕ ∂ 3 v · N = 0, we may use the following decomposition We easily have Integrating by parts in x 3 , we obtain, using |∂ n t v| Consequently, in light of the estimates (3.5), (3.15), (3.16), (3.27) and (3.39), we deduce from (3.4) that, by integrating in time from 0 to t, Recalling the definition (2.11) of∂ n t , sõ Then as Lemma 2.1, we may bound On the other hand, by the trace theory and using again that v · N = 0 on Σ − , we obtain Also, by Lemma 2.1, we have Therefore, by the estimates (3.43), (3.44) and (3.45), we deduce from (3.40) that where we have used Poincaré's inequality since Σ h = 0. Deriving the same estimates for the lower order time derivatives in a much simpler way, the estimate (3.3) then follows.
4. Vorticity and divergence estimates. In this section we will derive the estimates for the normal derivatives of v. The natural way is to estimate instead the vorticity curl ϕ v = ∇ ϕ × v to get rid of the pressure term ∇ ϕ q; applying curl ϕ to the first equation in (1.11), we have Note that the third and sixth equations in (1.11) implies that the equations (4.1) are characteristic in each of Ω ± , and the estimates follow in the same way as the one-phase incompressible Euler equations. For the completeness, we still repeat the procedure. We derive the following time-integrated energy estimates of the vorticity.
We now record an estimate for the divergence.
Proposition 4.2. It holds that Proof. It follows directly by the second equation of (2.25) and (2.30).
Proposition 4.3. It holds that Proof. It follows from (4.1) that where It is easy to see that F n 0 E n + ∂ n−1 t curl v 0 . We now take ψ ∈ H 1 (Ω), then we have By integrating by parts over Ω − , we have and Hence, we obtain which, together with the similar arguments on Ω + , implies that Similarly, we have ∂ n t div v (H 1 (Ω)) E n . By the trace estimates (A.14), (4.18) and (4.19), we conclude (4.10).
5. Elliptic regularity. In this section we will chain the Hodge-type elliptic estimates of v, the regularizing elliptic estimates of h by the presence of surface tension and the derivation of the estimates of q by using directly the fourth, fifth and first equations in (2.25) to derive the full energy estimates. We first present the result in the "energy".
Proposition 5.1. It holds that Proof. We first derive some preliminary estimates resulting from the control ofD n (weaker thanĒ n ). By the normal trace estimate (A.13), the second equation in Now integrating (6.6) in s from t/2 to t implies that, by (6.9), This together with (6.2) yields (6.3).
Appendix A. Analytic tools.
A.1. Poisson extensions. We will now define the appropriate Poisson sums that allow us to extend h ± , defined on the surfaces Σ ± , to functions defined on Ω, with "good" boundedness. Suppose that Σ + = T 2 × { }. We define the Poisson sum in T 2 × (−∞, ) by Here "−" stands for extending downward and " " stands for extending at x 3 = , etc. We can then extend h + to be defined on Ω by This allows us to extend h − to be defined on Ω − . However, we do not extend h − to the upper domain Ω + by the reflection since this will result in the discontinuity of the partial derivatives in x 3 of the extension. We instead to do the extension through the following. Let 0 < λ 0 < λ 1 < · · · < λ m < ∞ for m ∈ N and define the (m + 1) × (m + 1) Vandermonde matrix V (λ 0 , λ 1 , . . . , λ m ) by V (λ 0 , λ 1 , . . . , λ m ) ij = (−λ j ) i for i, j = 0, . . . , m. It is well-known that the Vandermonde matrices are invertible, so we are free to let α = (α 0 , α 1 , . . . , α m ) T be the solution to It is easy to check that, due to (A.5), ∂ l 3 P +,0 f (x , 0) = ∂ l 3 P −,0 f (x , 0) for all 0 ≤ l ≤ m and hence ∂ α P +,0 f (x , 0) = ∂ α P −,0 f (x , 0), ∀ α ∈ N 3 with 0 ≤ |α| ≤ m. (A.7) These facts allow us to extend h − to be defined on Ω by We may assume that m is sufficiently large for our use.
We have the following boundedness of P = P ± .