A semi-discrete scheme derived from variational principles for global conservative solutions of a Camassa–Holm system

We define a kinetic and a potential energy such that the principle of stationary action from Lagrangian mechanics yields a Camassa–Holm system (2CH) as the governing equations. After discretizing these energies, we use the same variational principle to derive a semi-discrete system of equations as an approximation of the 2CH system. The discretization is only available in Lagrangian coordinates and requires the inversion of a discrete Sturm–Liouville operator with time-varying coefficients. We show the existence of fundamental solutions for this operator at initial time with appropriate decay. By propagating the fundamental solutions in time, we define an equivalent semi-discrete system for which we prove that there exists unique global solutions. Finally, we show how the solutions of the semi-discrete system can be used to construct a sequence of functions converging to the conservative solution of the 2CH system.


Introduction
The Camassa-Holm (CH) equation is first known to have appeared in [25], although written in an alternative form, as a special case in a hierarchy of completely integrable partial differential equations. The equation gained prominence after it was derived in [8] as a limiting case in the shallow water regime of the Green-Naghdi equations from hydrodynamics, see also [18]. Since then, the CH equation has been widely studied due to its rich mathematical structure: it is for instance bi-Hamiltonian, admits a Lax pair and is completely integrable. The solutions may develop singularities in finite time even for smooth initial data, see, e.g., [12,13]. The so-called Camassa-Holm system (2CH) was first introduced in [41]. This is not the only two-component generalization which has been proposed for the CH equation. For instance, in [9,23] the authors showed how similar systems are related to the AKNS hierarchy. However, we will here only consider (1.2), which similarly to (1.1) can be derived as a model for water waves. Indeed, the system was derived in [20] from the Euler equations in the case of constant vorticity, while different derivation based on the Green-Naghdi equations can be found in [14]. The 2CH system shares many properties with the CH equation: the equation is bi-hamiltonian [41], admits a Lax pair and is integrable [14]. Results on the well-posedness, blow-up and global existence of solutions to (1.2) are provided in [21,22,33,34].
Both the CH equation and the 2CH system are geodesic equations, see [15][16][17]21]. The CH equation is a geodesic on the group of diffeomorphisms for the right-invariant norm To clarify this statement, we introduce the notation y : R + × R → R for a path in the group of diffeomorphisms, meaning that y(t, ξ) denotes the path of a particle initially at ξ, and the Eulerian velocity is given by u(t, x) = y t (t, y −1 (x)). The geodesic equation is then obtained as an extremal solution for the action functional The momentum map, as defined in [2], is given by the Helmholtz transform m(u) = u − u xx in Eulerian coordinates. Then we may write the energy as E kin = 1 2 R m(u)u dx. For the 2CH system in [21], the diffeomorphism group is replaced with a semi-direct product which accounts for the variable ρ. Then the 2CH system is a geodesic for the right-invariant norm 1 2 u 2 H 1 + 1 2 ρ 2 L 2 . However, we will not follow this approach here, but rather use the fact that (1.2) can be derived as the governing equation for a different variational problem, where the action functional includes a potential energy term and the variation is performed on the group of diffeomorphisms only. This point of view enables us to derive a discretization which mimics the variational derivation of the continuous case. In this approach, we consider the variable ρ as a density entering the action functional through a potential term (1.4) where ρ ∞ 0 is the asymptotic value of ρ. The mass conservation equation ρ t + (ρu) x = 0 is not a result of the variational derivation, but is instead a given constraint of the problem. Equation (1.4) can be interpreted as an elastic energy: it increases whenever the system deviates from the rest configuration given by ρ ≡ ρ ∞ . In the beginning of section 2, we present the derivation of the 2CH as the critical point for the variation of E kin − E pot . This approach follows the classical framework, see [1], and the potential term E pot depending on the density is added in the same way as in [40], see also [27,39,44] for applications to more complex fluids. In Lagrangian variables, the mass conservation equation simplifies to the expression To derive a discrete approximation of the 2CH system, we propose to follow the same steps of the variational derivation in the continuous case. First, we discretize the path functions y(t, ξ) by piecewise linear functions, y i (t) = y(t, ξ i ) for ξ i = iΔξ, i ∈ Z and Δξ > 0. Then, we approximate the Lagrangian using these discretized variables. Finally, we obtain the governing equation for the discretized system from the principle of stationary action, as in the continuous case. In our opinion, the advantage of using this variational approach as basis for our discretization is that we need only take variations with respect to a single discrete variable, rather than two. This reduction is achieved by the use of the identity (1.5). Note that the group structure of the diffeomorphisms is not carried over to the discrete setting, as the composition rule is not defined at the discrete level. In practice, this means that that our discretized equation will not have a purely Eulerian form and should be solved in Lagrangian variables. We retain two symmetries though, the time and space translation invariance. As a result, we have conservation of discrete counterparts of the integrals R (u 2 + u 2 x ) dx and R u dx, see section 2. We rewrite the 2CH system (1.2) in Lagrangian variables following [31]. We first apply the inverse of the Helmholtz operator Id − ∂ xx to obtain u t + uu x = −P x , P − P xx = u 2 + 1 2 We rewrite the second equation above as a system of first-order equations, for Q = P x and f = u 2 + 1 2 u 2 x + 1 2 ρ 2 . In Lagrangian variables we haveP(ξ) = P(y(ξ)), and the system (1.7) becomes In (1.8), the operator denoted by y ξ corresponds to pointwise multiplication by y ξ . The matrix operator corresponds to the momentum map in Lagrangian coordinates and must be inverted to solve the system. In contrast to its Eulerian counterpart in (1.7), the operator evolves in time. This significantly complicates the analysis, especially in the discrete case. In section 4, we introduce the operators G and K which define the fundamental solutions of the momentum operator, Note that the operator becomes singular when y ξ vanishes. In the discrete case, the momentum operator and its fundamental solution are given by where D ± denotes forward and backward difference operators, see section 2. This is a form of Jacobi difference equation, cf [43]. To establish solutions of (1.10), we shall invoke results from [24,42] which generalize the Poincaré-Perron theory on difference equations. Section 3 is completely devoted to this analysis. The CH equation and 2CH system can blow up in finite time, even for smooth initial data. The blow-up scenario for CH has been described in [11,12,19] and consists of a singularity where lim t→t c u x (t, x c ) = −∞ for some critical time t c and location x c . However, since the H 1 -norm of the solution is preserved, the solution remains continuous. In fact, the solution can be prolongated in two consistent ways: conservative solutions will recover the total energy after the singularity, while dissipative solutions remove the energy that has been trapped in the singularity, see [5, 6, 30-32, 36, 38]. If ρ > 0 initially, no blow-up occurs and the 2CH system preserves the regularity of the initial data, see [31]. We can interpret this property as a regularization effect of the elastic energy: the particles cannot accumulate at a given location because of a repulsive elastic force. The peakon-antipeakon collision is a good illustration of the dynamics of the blow-up. We present this scenario in figures 1 and 2. In the peakon-antipeakon solution, which corresponds to ρ 0 ≡ 0, we observe the breakdown of the solution and the concentration of the energy distribution into a singular measure. At collision time, u 2 + u 2 x = 0 and the energy reduces to a pure singular Dirac measure, which naturally cannot be plotted. For the same u 0 , but ρ 0 ≡ 1, the potential energy prevents the peaks from colliding, which is clear from the plot of the characteristics in figure 1. The potential energy grows as the characteristics converge and results in an apparent force which diverts them.
The global conservative solutions of the 2CH system are based on the following conservation law for the energy, (1.6). This equation enables us to compute the evolution of the cumulative energy defined from the energy distribution as in Lagrangian coordinates, for which we obtain dH dt = −(uR) • y. This evolution equation is essential to keep track of the energy when the solution breaks down. To handle the blowup of the solution, we need also to have a framework which allows the flow map ξ → y(t, ξ) to become singular, that is where y ξ can vanish and the momentum operator in Lagrangian coordinates become ill-posed. In [31], explicit expressions for P and Q are given. Here, we adopt a different approach where we propagate the fundamental solutions K and G from (1.9) in time. Introducing U = u • y, the equivalent system for (1.2) is given by with the evolution equations for K and G given by Here [U, K] denotes the commutator of U and K, see section 4. In the case where ρ ∞ = 0, R and Q in (1.12a) are given by (1.12c) The derivation of (1.12) can be carried over to the discrete system, and this is done in section 4. The short-time existence for the solution of the semi-discrete system relies on Lipschitz estimates. At this stage, one of the main ingredients in the proofs is the Young-type estimate for discrete operators presented in proposition 5.1. For the global existence, we have to adapt the argument of the continuous case and complement it with a priori estimates of the fundamental solutions (g, k, γ, κ). These estimates follow from monotonicity properties of these operators, see lemma 4.1. Section 5 is devoted to establishing existence and uniqueness for global solutions of the discrete 2CH system. In section 6, we explain how the solution of the semi-discrete system can be used to construct a sequence of functions that converge to the solution of the Solutions for peakon-antipeakon initial data. For ρ 0 ≡ 0 we plot u in (a) and u 2 + u 2 x in (b). We observe the blow-up of u x at t c ≈ 1.5 and the concentration of energy. For the same initial u 0 , but ρ 0 ≡ 1, we plot the corresponding solution in (c)-(e) and observe that u x does not blow up. In (e) we plot the distribution of the potential energy given by (ρ(t, x) − ρ ∞ ) 2 , and observe how it grows when the peaks get closer to each other.
2CH system (1.2). Finally, in section 7 we present how to construct appropriate initial data for the semi-discrete system in order to achieve the convergence in section 6.

Derivation of the semi-discrete CH system using a variational approach
As a motivation for our discretization, we will here outline how the system (1.2) can be derived from a variational problem involving a potential term, as indicated in the previous section. In a standard way, see, e.g., [1], the Lagrangian L is defined as the difference of a kinetic and potential energy where the energies are given by (1.3) and (1.4). The governing equation is then derived by the least action principle, also called principle of stationary action, on the group of diffeomorphisms. A first step in this direction is to introduce the particle path, denoted by y(t, ξ). Then, we rewrite E kin and E pot in Lagrangian variables. For the kinetic energy, we obtain From the principle of mass conservation for a control volume, the density satisfies 3) The definition of ρ given by (2.3) is equivalent to the conservation law (1.2b). We can check this statement directly: We introduce the Lagrangian density r defined as r(t, ξ) = ρ(t, y(t, ξ))y ξ (t, ξ), and by requiring it to be preserved in time, we impose the definition of the density ρ in the system. The identity (2.3) allows us to reduce the number of variables in our Lagrangian. Indeed, we rewrite the potential energy (1.4) in terms of the particle path y only and obtain which must satisfy the Euler-Lagrange equation Inserting the above identities in the Euler-Lagrange equation (2.6) we get For the remainder of this section we give details of how we discretize the variational derivation outlined above. Let us start by discretizing the kinetic and potential energies given by (2.2) and (2.4), respectively. First we divide the line into a uniform grid by defining ξ j = jΔξ for some discretization step Δξ > 0 and j ∈ Z. We approximate y(t, ξ j ) with y j (t) for j ∈ Z, and the spatial derivatives y ξ (t, ξ j ) with the finite difference D + y j . The finite difference operators D + and D − are defined as and they satisfy the discrete product rule When we later encounter operators in the form of grid functions with two indices, such as g i, j for i, j ∈ Z, we will indicate partial differences by including the index in the difference operator, for instance D j+ g i, j = (g i, j+1 − g i, j )/Δξ. We use the standard notation p and ∞ for the Banach spaces with norms Turning back to the energy functionals, we discretize the kinetic energy (2.2) using finite differences and set The Lagrangian velocity is as usual defined as U i =ẏ i and, using this notation, (2.10) becomes The discrete counterpart of (2.4) is similarly defined as (2.11) where y 0, j = y 0 (ξ j ) and ρ 0, j := ρ 0 (y 0, j ). Now we define the Lagrangian as the difference between the kinetic and potential energy, Now we compute the Fréchet derivatives of L dis with respect to y andẏ. This derivative is given in 2 , the space of square summable sequence using the duality pairing defined by the scalar product, v, w 2 := Δξ Formally, we have where in the final identity we have used the summation by parts formula This leads to the Fréchet derivatives where the rightmost equality in (2.13) is a consequence of E pot Δξ being independent of U. For the potential term we find which gives the Fréchet derivative The Euler-Lagrange equation is then see, e.g., [1]. From (2.14) we then have d dt which leads to the following system of governing equationṡ for j ∈ Z. Note that we have omitted ρ 2 ∞ on the right-hand side in (2.15) as D − maps constants to zero.
We can use the Legendre transform to define the Hamiltonian Writing out the above Hamiltonian explicitly we have We observe that the Lagrangian L dis does not depend explicitly on time. Then it is a classical result of mechanics, which follows from Noether's theorem, that H dis is time-invariant, The Lagrangian L dis is also invariant with respect to translation so that an other time invariant can be obtained. We denote by ψ : 2 × R → 2 the transformation given by the uniform translation (ψ(y, ε)) j = y j + ε. For simplicity, we write y ε (t) = ψ(y(t), ε). From the definition of ψ, we haveẏ ε (t) =ẏ(t) and D + y ε (t) = D + y(t).
Hence, the Lagrangian L dis is invariant with respect to the transformation ψ. Then Noether's theorem gives us that the quantity δL dis δU , δy ε δε 2 is preserved by the flow. In our case, δy ε δε j = 1 and δL dis δU j (2.13). Thus, we obtain that the quantity is preserved. Note that I corresponds to a discretization of in Eulerian coordinates, which is preserved by the 2CH system. Let us return to (2.15), and in particular to the left-hand side which containsU j , but not in an explicit form. For a given sequence a = {a j } j∈Z ∈ ∞ and an arbitrary sequence w = {w j } j∈Z ∈ 2 , we define the operator A[a] : 2 → 2 by When a = D + y, (2.19) corresponds to the momentum operator m in Lagrangian coordinates, and takes the form of a discrete Sturm-Liouville operator. This operator is symmetric and positive definite for sequences a such that a j > 0, as we can see from where we once more have used (2.12). When A[D + y] is positive definite, it is invertible and we may formally write (2.15) as a system of first order ordinary differential equations, (2.20) When solving the above system, we obtain approximations of the fluid velocity and density in Lagrangian variables, U j (t) ≈ u(t, y(t, ξ j )) and ρ 0, j /(D + y j (t)) ≈ ρ(t, y(t, ξ j )). We conclude this section with some comments on the Hamiltonian form of the equations. Hamiltonian equations in generalized position and momentum variables follow from the Lagrangian approach in classical mechanics, see, e.g., [1]. The generalized momentum is defined as p = δL dis δU (y, U). When we express the Hamiltonian H dis given in (2.17) in term of y and p, the Hamiltonian equations are then given as usual bẏ If we introduce the fundamental solution g i, j of the operator A[D + y], see section 3, the Hamiltonian can be rewritten as In the case ρ = 0 (that is E pot Δξ = 0), we recognized the similarity of this expression with given in [8]. The Hamiltonian H mp defines the multipeakon solutions, which can be seen as another form of discretization for the CH equation, see [37] for the global conservative case. Then, the two discretization appear as the results of two different choices of discretization for the inverse momentum operator: g i, j in the case of this paper andĝ i, j = e −|y i −y j | in [8]. We note that a numerical study of discretizations of the periodic CH equation considering both multipeakons and the variational method presented in this paper can be found in [26].

Construction of the fundamental solutions of the discrete momentum operator
In this section we construct a Green's function, or fundamental solution, for the operator defined in (2.19). Note that when a = D + y coincides with the constant sequence 1 = {1} j∈Z we have from (2.19) that A[1] = Id − D − D + , which corresponds to the operator used in the difference schemes studied in [10,35]. As the coefficients are constant, the authors are able to find an explicit Green's function g which can be written as and fulfills (Id − D − D + )g = δ 0 . Here δ 0 = {δ 0, j } j∈Z for the Kronecker delta δ i, j , equal to one when the indices coincide and zero otherwise. In our case, the coefficients appearing in the definition of A[D + y] are varying with the grid index j, which significantly complicates the construction of the Green's function. Let us consider the operator A[a] from (2.19) and the equation (A[a]g) j = f j . We want to prove that there exists a solution which decreases exponentially as j → ±∞. To this end, we want to find a Green's function for the operator A[a], and the first step is to realize that the homogeneous operator equation (A[a]g) j = 0 can be written as This can again be restated as a Jacobi difference equation, see [43, equation (1.19)], or equivalently in matrix form Observe thatÃ j is not symmetric and always contains positive, negative and zero entries under the assumption a j > 0. Moreover,Ã j is ill-defined when a j−1 = 0, which corresponds to the occurrence of a singularity in the system. We want to allow for this in our discretization in order to obtain solutions globally in time. If we go back to the first restatement of the operator equation and introduce the variable we get the following characterization of the homogeneous problem or equivalently Here A j is a symmetric matrix with positive entries whenever a j > 0, and it reduces to the identity matrix when a j = 0. We will use (3.5) rather than (3.2) to construct our Green's function, and it will also significantly simplify the analysis of the asymptotic behavior of the solutions.

Lemma 3.1 (Properties of matrix A j ).
Consider A j from (3.5) and assume a j = 1 + uniformly with respect to j when a j > 0. Moreover there is the obvious identity λ ± j = 1 when a j = 0. Asymptotically we have lim j→±∞ A j = A, where A is given by A j after setting a j = 1, and the eigenvalues λ ± of A satisfy for M > m > 0 depending only on Δξ. Moreover, as the eigenvalues are strictly positive it follows that the spectral radius of A j , spr(A j ) := max{|λ + j |, |λ − j |} satisfies A j = spr(A j ) = λ + j , A = spr(A) = λ + , and both matrices can be diagonalized: A j = R j Λ j R j , A = RΛR . Proof of Lemma 3.1. To see that det A j = 1 one can compute it directly, or see it from the eigenvalues which shows that A j is invertible irrespective of the value of a j . As A j is real and symmetric, it can be diagonalized with orthonormal eigenvectors r ± j as follows Since D + b ∈ 2 , for any j ∈ Z we have the bound meaning a j is bounded from above and below. Then it follows that corresponding to (3.6). Furthermore, since D + b j ∈ 2 , we have lim j→±∞ a j Δξ = Δξ. We denote by A, Λ, R, and λ ± the matrices and eigenvalues given by A j , Λ j , R j , and λ ± j after replacing a j by 1. From the preceding limit, (3.8) and (3.9) we obtain (3.10) Bounds for λ ± are obtained similarly to the bounds for λ ± j . As A j , A are symmetric and hence normal, their norms coincide with the spectral radius spr(·) which here coincides with the largest eigenvalue.
Note that (3.5) corresponds to a transition from (g j , γ j−1 ) to (g j+1 , γ j ), so that A j can be considered as a transfer matrix between these two states. Thus, solving the homogeneous operator equation (A[a]g) j = 0 bears clear resemblance to propagating a discrete dynamical system, and this is also the idea employed in the analysis of Jacobi difference equations in [43, equation (1.28)]. However, in making the change of variable to obtain (3.5) we lose the symmetry of the difference equation, and so the results in [43] are no longer directly applicable. On the other hand, our system can be regarded as a more general Poincaré difference system, and our idea is then to apply the results [24, theorem 1.1] and [42, theorem 1] to the matrix product which is the transition matrix from (g j , γ j−1 ) to (g k , γ k−1 ). Note that in the lemma below, the norms can be taken to be the standard Euclidean norm, but one could use any vector norm.

Lemma 3.2 (Existence of exponentially decaying solutions). Consider the matrix equation
That is, there exist solutions v n with exponential decay in either direction, owing to the Lyapunov exponent λ − < 1. Moreover, the initial vectors are unique up to a constant factor.

Proof of Lemma 3.2.
We begin with the case of increasing n, and we want to apply [24, theorem 1.2] which states that for sequences of positive matrices {A n } satisfying lim n→+∞ A n = A for some positive matrix A we have for some vectors v and w with positive entries such that Av = spr(A)v. As mentioned in [3, remark 4], there is in general no easy way of determining the vector w explicitly. We recall that our A n has positive entries, unless a n = 0 in which case we have A n = I. Because of (3.10), there can only be finitely many n 0 for which A n reduces to the identity. If we instead consider the sequence of positive matrices consisting of our {A n } where we have omitted the finitely many identity matrices, they clearly still satisfy (3.10) and so (3.14) holds with spr(A) = λ + and v = r + from (3.8) and (3.9). However, as the matrices we omitted were identities, it is clear that the limit in (3.14) for both sequences coincide. Hence, [24, theorem 1.1] holds for our nonnegative sequence as well. Now, as A n I entrywise it follows that the entries of Φ n,0 are nondecreasing for n 0, which means that Φ n,0 is also nondecreasing for such n. Therefore, by (3.14) we have that any initial vector v 0 such that w v 0 = 0 leads to a solution v n with nondecreasing norm, and which then by [42, theorem 1] must satisfy , an asymptotically exponentially increasing solution. Indeed, the nondecreasing norm rules out the possibility of v n = 0 for n large enough. It follows that (3.15) holds for equal to either λ + or λ − , but if it were λ − < 1, then v n could not be nondecreasing. However, choosing instead a nonzero v 0 such that w v 0 = 0, we obtain an asymptotically exponentially decreasing solution v n satisfying (3.15) with = λ − < 1. This follows by once more excluding the scenario of v n = 0 for large enough n, since v 0 is nonzero and each A n has full rank. Then the only remaining possibility is v n satisfying (3.15) with = λ − . An obvious For decreasing n, we will be able to reuse the arguments from above. From (3.11) we find that Φ n,0 is a product of inverses of A n for n < 0, and by (3.5) we have Since (A n ) −1 contains entries of opposite sign, it would appear that we may not be able to use our previous argument. However, a change of variables will do the trick for us. First recall (3.3) which shows that γ j corresponds to a rescaled forward difference for g j , hence its sign indicates whether g is increasing or decreasing at index j. For an increasing solution in the direction of increasing n it is then necessary for g n and γ n−1 to share the same sign as n → +∞. On the other hand, for an increasing solution in the direction of decreasing n, the forward difference for γ n−1 should have the opposite sign of g n as n → −∞. Therefore, a change of variables allows us to rewrite the previous equations as , n < 0 and for this system we may use the positive matrix technique from before. The eigenvalues of B j in (3.16) are the same as those of A j , but they switch positions in the corresponding eigenvectorsr ± j compared to r ± j of A j : The same argument as in the case of increasing n then proves the existence of v 0 giving exponentially decreasing solutions as n → −∞. The uniqueness follows from the uniqueness of limits in (3.14), which for a given eigenvector v of A means that w is unique up to a constant factor. But then again, since we are in R 2 , the vector orthogonal to w is unique up to a constant factor.

Remark 3.3 (Signs of the initial vectors).
Here we underline that the form of Φ n,0 implies that the entries of v ± 0 in lemma 3.2 must be nonzero, with opposite signs for v − 0 and same sign for v + 0 . Indeed, by (3.12) and (3.13) we have Let us then assume v − 0 = 0 with nonnegative entries of the same sign, for n 0, and so it is impossible for the norm to tend to zero. Hence, the entries of v − 0 must be nonzero and of opposite sign. For n → −∞, we can use (3.16) and the same argument to arrive at the same conclusion for [g 0 , −γ −1 ] , implying that v + 0 = [g 0 , γ −1 ] has nonzero entries of equal sign.

Theorem 3.4 (Existence of a discrete Green's function).
Assume {a j } j∈Z to be a nonnegative sequence such that a j = 1 + D + b j with D + b ∈ 2 . Then, for any given index i, there exists a unique sequence Proof. Our strategy follows the standard approach for constructing Green's functions: we first construct solutions of the homogeneous version of (3.17) with exponential decay, and then we combine them in order to obtain a delta function at a given point. We start by constructing g 0, j centered at i = 0.
Choosing v ± 0 from lemma 3.2 we set and define the sequences where by construction g ± , γ ± have exponential decay for n → ∓∞. Then, applying the operator A[a] to g ± we find by construction of g ± and γ ± . Let us then define for some hitherto unspecified constant C, and observe from the homogeneous equation that a j g 0, j − D − γ 0, j = 0 for j = 0. Moreover, we have D + g 0, j = a j γ 0, j for all j by construction. Now we would like to show that the constant C can be chosen to obtain A[a]g 0,0 = 1/Δξ. From (2.12), we get  22) and the last equality follows from the identity g ± n+1 = g ± n + Δξa n γ ± n . Since the left-hand side of (3.21) vanishes by definition of g ± , we have W n (g − , g + ) = W m−1 (g − , g + ) for any n, m ∈ Z. That is, the Wronskian W n (g − , g + ) is a constant W(g − , g + ) for the constructed sequences g + and g − .
Next, we want to show that the Wronskian is nonzero. Considering and the definition (3.18), we use the sign properties stated in remark 3.3 to conclude that the two terms in the final sum are always nonzero and of the same sign, implying W(g − , g + ) = 0. Finally, we will determine the constant C by considering the backward difference which leads to Consequently, setting C −1 = W(g − , g + ) in (3.20) gives the desired Green's function.
Note that there is nothing special about the index i = 0 where we centered the Green's function. We can simply use the sequences (3.19) from before and define to obtain a Green's function g i, j centered at an arbitrary i. The uniqueness of g i, j follows from the vectors v ± 0 in lemma 3.2 being uniquely defined up to constant factors. Indeed, when constructing the Green's function in (3.23) these factors disappear since we are dividing by the Wronskian W(g − , g + ), and so we have no degrees of freedom left in our construction of g i, j , hence it is unique.

Note that A[a]
is not the only way to discretize the operator with first order differences, we may also consider In fact, we will need the Green's function for this operator as well to close our upcoming system of differential equations. Fortunately, the existence of Green's function for (3.24) follows from the considerations already made in theorem 3.4.
Proof of Corollary 3.5. Manipulating the homogeneous version of (3.25) we find it to be equivalent to Introducing the previous equation can be written as where we recognize the matrix A j from (3.5). Going backward we find (3.16). Hence, we get the solution for free from (3.4). Indeed, choosing we know that these sequences have the correct decay at infinity. Defining are exactly those found in (3.1) for the operator Id − D − D + . In fact, for a j ≡ 1 the sequences g and k coincide since D − D + = D + D − , and their explicit expression (3.1) can be recovered from the columns of Λ n R −1 in the diagonalization A n = RΛ n R −1 . Observe that by (3.3) and (3.26) we can rewrite (3.17) and (3.25) in the compact form (3.28)

Lemma 3.7 (Sign properties of the discrete Green's functions).
Assume a j 0 for j ∈ Z, and let g, γ, k, and κ be solutions of (3.28) which decay to zero for | j − i| → +∞. Then the following sign properties hold, In particular, this leads to the monotonicity properties where the arrows denote monotone decrease.
In figure 3 we have included a sketch of g i,n , γ i,n , k i,n , and κ i,n for Δξ = 0.2, i = 0, 4 and a n = a(nΔξ) given by We say sketch, as they have been computed on a finite grid n ∈ {−20, . . . , 20} with boundary conditions γ i,−21 = g i, 21 = 0 and k i,−21 = κ i,21 = 0, and consequently we find that neither of −20 or k i,20 are exactly zero. However, the exponential decay makes them very small and the qualitative behavior indicated in lemma 3.7 is still the same. Note how a(ξ) being zero on the interval (0.5, 1] leads to constant kernel values in that neighborhood, even at the peaks for the kernels centered at ξ 4 = 0.8.

Proof of Lemma 3.7.
We prove this only for g and γ as the proof for k and κ is similar. The proof relies on the reasoning in remark 3.3.
As a first step we want to show that the properties (a) and (b) hold for g i,i , g i,i+1 , γ i,i−1 , and γ i,i . To this end, we recall from the proof of theorem 3.4 that since g i, j and γ i, j satisfy (3.28), they must also satisfy with A k and B k as defined in (3.5) and (3.16). By our assumptions, the Green's functions must tend to zero asymptotically, and we recall from remark 3.3 that a necessary condition for this is for the vectors where we stress the importance of a j 0 for this argument to hold. Using only (3.28) we calculate Since a j 0, it follows that g i,i 0. Recalling that g i,i must be nonzero according to the sign requirements, we necessarily have g i,i > 0, and then γ i,i−1 > 0 follows. Moreover, multiplying the identity g i,i+1 − Δξa i γ i,i = g i,i by g i,i+1 and using a i 0, g i,i > 0, and g i,i+1 γ i,i < 0, we must have g i,i+1 > 0, which then implies γ i,i < 0.
Next we must prove that (a) and (b) hold for the remaining values of j, and this will be achieved with a contradiction argument. We define the vectors v + j := If we can prove that they retain the sign property under the above propagation, then we are done. Let us consider Assume that v + j does not retain the sign property, then there is some k i + 1 which is the first index such that v + k+1 does not have a positive first component and negative second component. We consider the two possible cases.
The first case is v + k+1 0 (v + k+1 0) considered entrywise. First of all, v + k+1 cannot be the zero vector as A k has full rank, since then v + k would also have to be zero which contradicts k + 1 being the first problematic index. Otherwise, the entrywise inequality A k+1 I leads to ). This is however impossible, as it contradicts the assumed decay of the Green's functions.
The remaining case is that the entries interchange sign from v + k to v + k+1 . However, then we would have Since a k 0, v + k would also have negative first component and positive second component, which contradicts k + 1 being the first problematic index. Hence, v + j always has positive first component and negative second component for j > i, thus for j i it follows that g i, j is always positive, while γ i, j is always negative which shows that g i, j is decreasing in this direction.
A similar argument holds in the other direction when considering v − j and B j . Then −γ i, j is always negative for j < i, which means that g i, j is increasing with j for these indices. Thus, (a) and (b) hold for {g i, j } j∈Z and {γ i, j } j∈Z .

An equivalent semi-discrete system for global solutions in time
We now return to the initial value problem (1.2). We use the Lagrangian formulation introduced in earlier works, see [31], but reformulate the governing equations by propagating the fundamental solutions of the momentum operator.

Reformulation of the continuous problem using operator propagation
The 2CH system can be written as for P implicitly defined by Let us introduceρ := ρ − ρ ∞ to ease notation. Note that most expressions simplify when we consider ρ ∞ = 0. We have chosen to cover the case of arbitrary ρ ∞ , to allow for the initial condition ρ(0, x) = ε, for any ε > 0. Such initial data lead to solutions without blow-up, see [29]. In the case of the 2CH system, the conservation law for the energy is given by where we have used P from (4.1) to define We can check that the first order system is equivalent to (4.1). Hence, and (4.4) is yet another form of the 2CH system. We introduce as before the Lagrangian position y(t, ξ) and velocity U(t, ξ). Moreover, we define the Lagrangian density r(t, ξ): = ρ(t, y(t, ξ))y ξ (t, ξ), and the cumulative energy H given by as well as the Lagrangian variablesQ = Q • y andR = R • y. From (4.5), we get U t = −Q and r t = 0, while the conservation of energy (4.2) yields H t = −UR. Finally, we rewrite the system (4.4) in terms of the Lagrangian variables. To simplify the notation, we replaceQ by Q, and similarly forR. The equivalent system in Lagrangian variables is then given by . (4.8) In (4.8), we use the same notation for the variable y ξ and the operator for pointwise multiplication by y ξ . We will use this convention for the rest of the paper. The equivalence between (4.4) and (4.8) holds only assuming the that y ξ 0 and all the functions are smooth enough to do the manipulation. Note that we need to decompose the variables y and r in (4.7) to give them a decay which enables us to define them in a proper functional setting. We define ζ andr as y(t, ξ) = ζ(t, ξ) + ξ and r(t, ξ) =r(t, ξ) + ρ ∞ y ξ (t, ξ).
The Banach space which contains ζ and H is the subspace of bounded and continuous functions with derivative in L 2 , endowed with the norm f V := f L ∞ + f ξ L 2 . Then we let be a Banach space tailored for the tuple X = (ζ, U, H,r) with norm The unique solution of (4.7), as studied in [31], is then completely described by this tuple.
An alternative viewpoint of the equivalent Lagrangian system is the following. Let us define the operators G and K as the fundamental solutions to the operator in (4.8), meaning that they satisfy As we mentioned in the introduction, the operators K and G can be computed explicitly, using the fundamental solution of the Helmholtz operators in Eulerian coordinates. If we define and then we can check that the operators defined as G( f )= R g(η, ξ) f (η)dη and K( f )= R κ(η, ξ) f (η)dη are solutions to (4.12), again assuming y is monotone increasing in ξ. This means that we can obtain explicit expressions for R and Q given by In [31,36], the authors prove that the right-hand side of their respective versions of (4.7) is locally Lipschitz, and consecutive contraction arguments yield the existence of a unique shorttime solution. In the same manner, we would like to prove that there exists a unique shorttime solution for our semi-discrete system, but the explicit forms for g and κ in (4.13) are not available in the discrete setting. As a remedy, we propagate the kernel operators corresponding to K and G by incorporating them in the governing equations. Given the evolution of y, that is, y t = U, we can derive evolution equations for G and K. Let us see how this can be done in the continuous case before dealing with the discrete case. Formally we have Here we assume again that we know a priori that y remains a monotone function with respect to ξ. Then, we can rewrite the last equality as For a function U, we can associate a pointwise multiplication operator which we denote by U.
That is, we may write U( f )(ξ) = U(ξ) f (ξ) for any function f and any point ξ. The integral kernel of U would be singular and equal to U(ξ)δ(ξ − η). Using this notation, we can rewrite (4.15) as This can equivalently be stated as where the evolution equation for K is derived analogously. An equivalent system of equations for the 2CH system is then given by with R and Q given as For all initial conditions we will consider, the new system (4.17) and (4.18) gives rise to the same solutions as the one given by (4.7), (4.13) and (4.14). It can be shown that the evolution equation (4.17) for G and K can be obtained directly from the product identity (4.12) by differentiating it and using integration by parts.

Reformulation of the semi-discrete system
Turning back to the formal expression (2.20), we use the the Green's functions from theorem 3.4 and corollary 3.5 to write out the right-hand side explicitly. Considering (3.28) where we now have a j = D + y j , we observe that they correspond to the discrete versions of (4.12). Indeed, we have the following identity which has to be compared with (4.12) in the continuous case. Thus, the second equation in (2.20) can be rewritten aṡ where we have defined and From the expressions in (4.20) and (4.22), it seems that, if D + y i goes to zero for some index i and time t, thenU j and h i blow up. However, it turns out that these quantities remain bounded, which allows us to extend the solution globally in time. To obtain a well-defined system, we are going to remove the explicit dependence on 1/D + y i by adding h to the set of variables of the system. With the discrete kernels g, k, γ, and κ from section 3 we are able to express A[D + y] −1 in (2.20) to obtain (4.20). However, since we do not know their explicit form as functions of D + y j , we derive a system analogous to (4.17) by introducing g, k, γ, and κ as variables. To compute the evolution of g, k, γ, and κ, we repeat the procedure from the continuous case. By differentiating (4.19) and using the fact thatẏ i = U i , we get Here we denote by (g * f ) j the action of the operator g i, j as a summation kernel on a sequence f i , defined as Moreover, we introduce the following norms for the operators, We establish in the next lemma some important properties for the fundamental solutions.

Lemma 4.1 (Preservation of identities).
Let T > 0, and assume that, for t ∈ [0, T], (D + y j (t)) t = D + U j (t) for j ∈ Z, and that g, k, γ, κ and D + U are bounded in 2 -norm in [0, T]. Then, for t ∈ [0, T] the sequences g i,j (t), k i,j (t), γ i,j (t), κ i,j (t) satisfy the following identities: (a) We introduce the four operators z l for l = 1, 2, 3, 4 defined as Using (D + y j (t)) t = D + U j (t) and (4.23) we find that Similarly, one shows that Integrating the first of these, taking absolute values, applying Hölder's inequality and taking supremum over i we obtain Treating the three other relations similarly and defining we may add the four inequalities to obtain an inequality of the form is bounded by assumption. Since Z(0) = 0, Grönwall's inequality yields Z(t) = 0 for t ∈ [0, T], which proves the result. (b) We prove the symmetry of g.
Δξ δ i,m , such that summation by parts shows Then, we use the identity D m+ g j,m = (D + y m )γ j,m from (4.19) twice, first for j and then for i, to obtain After summation by parts and using (4.19) once more, we end up with g j,i = Δξ m∈Z g i,m (D + y m )g j,m + D m− γ j,m = g i, j , and the symmetry of g is proved. A similar procedure shows the symmetry of k i, j . For the antisymmetry we also use (4.19) and summation by parts to compute Returning to (4.20), the second term in the right-hand side of the governing equations can be simplified as follows, where we have used (4.19) and (4.27). We define Then, the evolution of U is given bẏ The form of Q in (4.28) motivates the definition Indeed, with these definitions we have We recognize this as the discrete version of (4.8).
The relationU j = −Q j shows that we have a differential equation for U in the variables y, U, H,r, g, and κ. From (4.21) we obtaiṅ (4.32) Next, we introduce the cumulative energy H j as after using (4.29) and (4.32). Then, we use the relation between Q and R given in (4.31) to obtain Simplifying further, we obtaiṅ This leads toḢ where in the last equality we have used the decay at infinity together with (2.12). Collecting all the equations and applying the relations (4.26) and (4.27) we obtain the closed systemζ where y j = jΔξ + ζ j , and we recall

Existence and uniqueness of the solution to the semi-discrete 2CH system
In this section, we show that the semi-discrete system (4.35) has a unique, globally defined solution. Let us first introduce the functional setting for the analysis. We define the discrete analogue of the H 1 (R)-inner product, which induces a norm in the usual manner. The discrete Sobolev-type inequality can be proven in a very similar way as its continuous version, see, e.g., [7]. We introduce the subspace V Δξ of ∞ defined as We define the discrete version of the space used in the continuous setting, namely Since we have included the operator kernels as solution variables in (4.35), we have to introduce a space for them as well. To account for that the kernels are well-behaved, we choose their space to be * := 1 ∩ ∞ with norm · * = · 1 + · ∞ , with p -norms defined in (4.25). We note that * ⊂ 2 , since we have the inequality Thus, we will consider solution tuples of the form where E ker Δξ denotes the space E Δξ augmented with the space for the kernel operators * . Moreover, for the kernel operator g we have that the transpose g of g is given by (g ) i, j = g j,i . Then, the following result, reminiscent of Young's convolution inequality, will prove useful.

Proposition 5.1 (Young's inequality for general operators)
Above, we use the convention q/∞ = 0 for q < ∞, and ∞/∞ = 1. Note that the standard Young's inequality is usually given for a translation invariant kernel where g takes the form g i, j =ĝ i− j for some sequenceĝ. For an operator of this form, we can check that g = τ • g • τ , where the operator τ inverts the indexing, that is τ ( f ) j = f −j . Since the operator τ is an isometry in all q -spaces, the expression (5.6) simplifies to Proof of Young's inequality. Below we will use the following discrete version of the generalized Hölder inequality, where the jth component of a product of sequences is interpreted as n k=1 a k j = n k=1 (a k ) j . We note that the proof of (5.7) follows that of the continuous case, see, e.g., [7, exercise 4.4]. Let us denote h = g * f . Note that r < ∞ =⇒ p, q < ∞, which shows that some configurations are impossible and can be excluded. We deal with the three remaining cases: (a) r < ∞: from the generalized Hölder inequality we obtain where we have used Fubini's theorem in the second inequality. Taking rth roots we obtain the result.
To prove the short-time existence of (4.35), we consider an auxiliary system which corresponds to (4.35), except that we have decoupled ζ, U and H from their discrete derivatives D + ζ, D + U and D + H by introducing the sequences α, β and h. The reason for this is that we cannot take for granted that the kernels satisfy (4.19) for t > 0, and then we cannot use (4.31) when estimating the right-hand side of (4.35b) in h 1 -norm. Once the short-time existence of solutions to the auxiliary system is established, we will prove that the coupling between y, U, H and their discrete derivatives is indeed preserved if it holds initially. The auxiliary system readsζ where we have momentarily redefined R and Q as The evolution equations (5.8c) and (5.8d), and the second equation of (5.8b) have been obtained formally by applying D + to (4.35a), (4.35b) and (4.35c), in combination with (4.31). We collect all the variables in a tuple Δξ and introduce the corresponding norm Note how we require U ∈ ∞ to account for the fact that the decoupling of U and D + U deprives us of the continuous inclusion h 1 ⊂ ∞ .

Lemma 5.2 (Short-time solution for (5.8)). Let Y 0 ∈ E aux
Δξ be such that 1 + α j 0 for all j, and with initial auxiliary variables g 0 , k 0 , γ 0 , κ 0 constructed according to theorem 3.4 and corollary 3.5 with a j = 1 + α j . Then, there exists a time T > 0 depending only on Y 0 E aux Δξ such that (5.8) has a unique solution Y ∈ C 1 ([0, T], E aux Δξ ) with initial data Y 0 .
Proof of Lemma 5.2. We are going to use the symmetry and anti-symmetry identities (4.26) and (4.27) in our estimates and we explain now why it can be done. First, we note that these identities hold initially by the construction of (3.23) and (3.27). Then, from the evolution equations (5.8e)-(5.8h) one can check that the symmetry identities are preserved by the Picard fixed-point operator which we will use here to prove the short-time existence of (5.8). Then, by establishing local Lipschitz regularity of the right-hand side, we can prove the existence of a short-time solution in the closed subset of E aux Δξ where (4.26) and (4.27) hold. Let us consider two functions in E aux Δξ , U, H, r, α, β, h, g, k, γ, κ) andỸ = ζ ,Ũ,H,r,α,β,h,g,k,γ,κ .
For the Lipschitz estimates, we first treat the right-hand sides of (5.8e)-(5.8h). We only provide details for (5.8e) as (5.8f)-(5.8h) can be treated similarly. We start by considering the ∞ -norm using the following splitting, We estimate the first term as follows Δξ m∈Z β m g j,m g i,m − Δξ m∈Zβ mg j,mgi,m g ∞ g 2 β −β 2 + g ∞ β 2 g −g 2 + β 2 g 2 g −g ∞ and the second term has a similar estimate. For the 1 -norm, use the same splitting and consider again only the first term. We make use of the symmetry properties of the kernel operators, as given in lemma (4.1), to switch between indices and obtain Δξ i∈Z Δξ m∈Z β m g j,m g i,m − Δξ m∈Zβ mg j,mgi,m g 1 g 2 β −β 2 + g 1 β 2 g −g 2 + β 2 g 2 g −g 1 , From (5.5) we can then conclude that the right-hand side in (5.8e) is locally Lipschitzcontinuous with respect to the E aux Δξ -norm. Let us consider Lipschitz properties of R and Q. We decompose Q in Q 1 + Q 2 where Starting with Q 2 , we have For the first term above, applying the Young's inequality (5.6) with r = p = 2 and q = 1, we get Using the antisymmetry property (4.27) of κ andκ, namely κ = −γ andκ = −γ, we get Hence, we obtain the following estimate in 2 -norm, For the ∞ -norm, we use the same splitting Applying (5.6) for r = ∞ and p = q = 2, and the symmetry property of κ, we obtain in a similar way as before that In a similar fashion as for Q 2 we find Furthermore, analogous applications of (5.6) and (4.27) produce For the 1 -norms above we then apply the Cauchy-Schwarz inequality to obtain which contain the relevant norms. From the preceding estimates on Q 1 and Q 2 the local Lipschitz property of the right-hand side of the second equation in (5.8a) in the 2 ∩ ∞ -norm is clear. Furthermore, since U ∈ ∞ , the previous ∞ -estimates on R and Q also show that the right-hand sides of (5.8c) and (5.8d) are locally Lipschitz in the 2 -norm. For the last equation in (5.8a), we introduce the right-shift operator (τ R) j = R j−1 and we have The remaining right-hand sides of (5.8a) and (5.8b) are linear in the solution variables, and thus Lipschitz in their respective norms. Hence, for (5.8) written asẎ =F(Y) we have which is what we set out to prove.
The final step in obtaining short-time existence for (4.35) from the auxiliary system, is to show that if the initial data for (5.8) satisfy then these identities are preserved in time by the solution. The result for (5.11a) has been proved in lemma 4.1, as it only depends on the identity (D + y) t = D + U, which is replaced here byα = β. Using (5.11a), we infer from (4.31) that From the definition of (5.8) we get while the expression for D + Q j from (5.12) yields and from the expression for D − R j we obtain Hence, the equation (5.13) give us that (5.11b) holds for all time if it holds initially. Then we have proved the following theorem. -time solution for (4.35)). Given X 0 ∈ E ker Δξ such that 1 + D + ζ j 0 and g 0 , k 0 , γ 0 , and κ 0 are constructed according to theorem 3.4 and corollary 3.5 with a j = 1 + D + ζ j . Then, there exists a time T depending only on X 0 E ker Δξ such that (4.35) has a unique

Theorem 5.3 (Short
The next step is to prove that there exists a subset, denoted by B, of E ker Δξ which is preserved by the evolution equation. For this subset, the solution exists globally in time. The subset B is defined as follows. Definition 5.4. The set B is composed of all (ζ, U, H,r, g, k, γ, κ) ∈ E ker Δξ such that (a) g, k, γ, κ satisfy the properties listed in lemma 4.1 for a = D + y, j for all j, (d) D + y j 0, D + H j 0, D + y j + D + H j > 0 for all j.

Lemma 5.5 (Properties preserved by the flow). Given initial datum X
Δξ ) be the corresponding short-time solution given by theorem 5.3. Then Proof of Lemma 5.5. Property (a) follows from lemma 4.1, since the solution variables in X(t) satisfy D +ẏ j = D + U j and D + U ∈ 2 , where we as usual have D + y j = 1 + D + ζ j . The proof of property (b) essentially follows [31, lemma 3.3], and so we omit it. The proof of (c) is similar to the proof of [31, lemma 3.5], while the proof of (d) is analogous to that of [36, lemma 2.7], and they are also omitted here.
For the rest of the paper we will only consider X ∈ B ∩ E ker Δξ , as solutions in this set contains all the relevant solutions to the original 2CH system (1.2). Lemma 5.5, and in particular the preservation of the identity 2(D + y j )h j = U 2 j (D + y j ) 2 + (D + U j ) 2 +r 2 j (5.14) allows us to prove useful estimates for the solutions in B. We have where H ∞ (t) = lim n→+∞ H n is the total energy of the discrete system. This quantity corresponds to H dis in (2.16). Indeed, the Hamiltonian (2.17) is We denote the preserved total energy H ∞ (t) by H ∞ . Turning back to the inequality (5.15), it can be proved as follows, where in the first inequality we have used (5.14), and in the second inequality we have used D + y j 0 together with the Cauchy-Schwarz inequality. An immediate consequence of (5.15) is that U ∞ can be uniformly bounded by a constant depending only on H ∞ . To show this, we add and subtract in (2.8) to find the identity Taking advantage of the decay of U at infinity, we may then write follows. From (5.16) and (4.35a), we then obtain the estimate Another useful estimate coming from (5.14) is Now that lemma 5.5 has established D + y j (t) 0 in the short-time solution for t ∈ [0, T], we can apply lemma 3.7 with a j = D + y j . Indeed, the sequences g, γ, k, and κ solve (4.19) and belong to * for t ∈ [0, T], and so they correspond to the unique decaying solution. These properties contained in lemmas 3.7 and 4.1 are essential to establish the a priori estimates contained in the next lemma.

Lemma 5.6 (A priori relations and inequalities for the kernels). As a consequence of establishing the preservation of the summation kernels and their sign properties over time, we have the identities
as well as and the bounds Proof of Lemma 5.6. To prove (5.19a) we use D + y j 0 and (4.19) for the leftmost equalities, while for the rightmost equalities we use the monotonicity properties of (3.29) to write We obtain (5.19a) in the same way. To obtain (5.20), we use the definitions of the operators A in (2.19) and B in (3.24), and apply telescopic cancellation to the differences D j+ γ i, j and D j− κ i, j in the identities (4.19). In the same manner, telescopic cancellation applied to (4.19) yields Using the fact that D + y j 0 and g i, j 0, the triangle inequality and (5.20) yield (5.21) for γ. We proceed similarly for κ. For g, observe that, using (4.19), we can rewrite them as Using the decay at infinity we can then write where we have used (3.29) for the first inequality, and (5.23) for the final identity. The bound g i,i 1 follows, and we use 0 g i, j g i,i from (3.29) to conclude. A similar procedure can be applied to prove that k i, j 1. Furthermore, we have and the result on the 1 bound of g follows from (5.19) and (5.21). A similar procedure proves the bound on k 1 . For the bound on γ 1 we find where in the second equality we use lemma 3.7, the third equality uses summation by parts (2.12), and the fourth is due to the kernel definition property (4.19). Then the result follows from (5.20), (5.21), and 0 γ i,i−1 1. A similar procedure proves the bound on κ 1 .
A direct consequence of (5.21) is that the ∞ -norms of the kernels remain bounded by 1 for all time. Moreover, combining (5.22) with (5.17) we find that the 1 -norms remain bounded for any finite t, namely Furthermore, lemma 5.6 allows us to find a bound similar to (5.16) for R ∞ and Q ∞ . Indeed, for Q we find Using (5.18) and the Cauchy-Schwarz inequality, we have which by (5.19) and (5.21) simplifies to Using (5.15), we get UD + U 1 H ∞ . Hence, from (5.25), we get An analogous estimate for R can be obtained so that we can conclude with the bounds Now we are set to prove global existence for solutions of (4.35).

Theorem 5.7 (Global existence).
Given initial datum X 0 in the set B from definition 5.4, the system (4.35) admits a unique global solution X ∈ C 1 ([0, ∞), E Δξ ), such that X ∈ B for all times. In particular, for t > 0, the norm X(t) E Δξ is bounded by C X(0) E Δξ for a constant C depending only on t, the total energy H ∞ , the asymptotic density ρ ∞ , and ζ(0) ∞ .

Proof. The solution has a finite time of existence T only if
blows up as t approaches T. Otherwise the solution can be prolonged by a small time interval by theorem 5.3. Let X be the short-time solution given by (5.3) for initial datum X 0 . We will prove that sup 0 t T X E Δξ < ∞. From the definition of the h 1 -norm and (5.2) we find that the right-hand side of (4.35a) is bounded in the V Δξ -norm by 2+ √ 2 2 U h 1 , while the right-hand side of (4.35d) is bounded in 2 -norm by ρ ∞ U h 1 . Next, we estimate the right-hand side of (4.35b), where we have used the definition of the h 1 -norm, (4.31) and the decomposition D + y j = 1 + D + ζ j . Then, recalling the definitions (5.9a) and (5.9b) and applying the Young inequality (5.6) to the final expression above we see that it is bounded by Then, applying (5.15), (5.16), (5.24) and (5.26) and the definitions of the V Δξ -and h 1 -norms we obtain that Q h 1 is bounded by Finally, the V Δξ -norm of the right-hand side of (4.35c) can be estimated as where we again use the notation (τR) j = R j−1 . In the first identity above we have employed (4.31), while in the final line we have used the definitions of the V Δξ -and h 1 -norms together with (5.2), (5.16) and (5.26). Gathering all the above estimates of the right-hand sides, writing (4.35) in integral form, and taking norms we obtain the following inequality for X(t) = (ζ, U, H,r)(t), for t ∈ [0, T] and some constant C(H ∞ , ζ(0) ∞ , ρ ∞ ) depending only on H ∞ , ζ(0) ∞ and ρ ∞ . Grönwall's inequality then yields which shows that X(T) E Δξ is bounded, and we may according to theorem 5.3 extend our solution indefinitely. In retrospect, with the estimates (5.16) and (5.26) in hand, we can apply a Grönwall estimate to the evolution equations for α, β, h, andr in (5.8). From this we find that the ∞ -norms of D + y, D + U, h, andr at time t ∈ [0, T] are bounded by their ∞ -norm at time t = 0 times a factor exp{C(H ∞ , ρ ∞ )t}, where the constant C(H ∞ , ρ ∞ ) depends only on H ∞ and ρ ∞ .
As mentioned in the introduction, if ρ > 0 initially for the 2CH system (1.2), then the smoothness of the initial data is preserved, see [31]. This is because the characteristics do not collide in this case, and y ξ remains positive for all time. In the discrete case, this property takes the form of a lower bound on D + y. For any given time T, there exists a constant C > 0 depending on max t∈[0,T] X(t) E Δξ , ρ ∞ , and T such that for all j and t ∈ [0, T]. This follows from property (c) in definition 5.4. Thus, if ρ 0, j > 0, we will have y j (t) < y j+1 (t) for all time.

Convergence of the scheme
In this section we interpolate the solutions of the semi-discrete scheme analyzed in section 5, with initial data constructed in section 7. We shall then show that these interpolated functions converge to the solution of the 2CH system as written in (4.7) and (4.8). Let us in this section use Y Δξ to denote the tuple of grid functions obtained in theorem 5.7 for t ∈ [0, T], Since these functions also belong to the set B in definition 5.4, we will augment the E Δξ -norm (5.4) as follows: In order to ease notation below, we will write Y Δξ for sup 0 t T Y Δξ (t) B . We define the interpolated functions as follows where V is a placeholder for ζ, U, H, and Q, while χ j (ξ) denotes the indicator function for the interval [ξ j , ξ j+1 ). We also introduce the functions Observe that the interpolated functions above are piecewise linear and continuous, except for r Δ ,r Δ which are piecewise constant. In particular we note the identity which shows R Δ (t, ξ j ) = R j−1 . Let us also recall the definition of the space E in (4.10). A consequence of theorem 5.7 is that tuple of interpolated functions satisfies X Δ (t) ∈ C 1 ([0, T], E) for any fixed T > 0 and Δξ > 0. Let us now consider a given initial datum X 0 = (ζ 0 , U 0 , H 0 ,r 0 ) ∈ E for the equivalent 2CH system (4.7). We assume there exists a sequence of discrete initial data Y Δξ,0 ∈ E Δξ such that the interpolation of Y Δξ,0 , denoted X Δ,0 , converges to X 0 , i.e., lim Δξ→0 X Δ,0 − X 0 E = 0. (6.5) We will explain how to construct such sequence in the next section. For T > 0 and each Y Δξ,0 , let Y Δξ be the corresponding solution given by theorem 5.7. Furthermore, we denote by X ∈ C([0, T], E) the solution to (4.35) with initial data X 0 , while X Δ ∈ C([0, T], E) is the function interpolated from Y Δξ using (6.2). Then, we have the following convergence result.
due to (4.31), (4.35a), (4.35b), and (4.35d). Thus, the three linear equations in (4.7) are satisfied exactly by our interpolants. The next step is to control the evolution of the error for the variable H Δ . We find This identity then implies almost everywhere. Combining the above identities we can estimate the error in the V-norm as follows, Now, for the relations (4.8), we measure the error in L 2 -norm. From (4.31), we obtain the relation and find Finally, using (4.31) once more, we have which can be estimated as The estimate (6.6) is exactly as we want it, (4.7c) is satisfied in the appropriate norm up to some small remainder. However, the estimates (6.7) and (6.8) require some more work, as we shall see next.
Let us estimate the E-norm of the difference between X Δ (T ) and the exact solution X(T) := (ζ, U, H,r)(T). From the above estimates and (4.7) we find where we have used that the final expression in (6.6) can be bounded by ΔξC H ( Y Δξ ) for some constant C H depending only on Y Δξ . From (6.9), it is clear that we need estimates of and by definition of the H 1 -norm and the Sobolev inequality f L ∞ 1 √ 2 f H 1 , it will be sufficient to bound Q Δ − Q H 1 and R Δ − R H 1 . To this end, we note that by the estimates (6.7) and (6.8), it follows that for some functions v Δ , w Δ ∈ L 2 which are bounded by a constant depending only on the norm Y Δξ of (6.1). Recalling (4.14) and the operators defined in (4.12) we know that R(t, ξ) and Q(t, ξ) can be written as Due to the obvious similarities between (6.10) and (4.18) we would like to generalize the operator identity (4.12) by replacing y(t, ξ) with any function b(t, ξ) such that b(t, ·) − Id ∈ V and b ξ (t, ξ) 0, in particular this holds for our y Δ (t, ξ) in (6.3) by virtue of lemma 5.5. This is can be done, and the unique H 1 -solution of for v(t, ·), w(t, ·) ∈ L 2 is then Consequently, we can generalize G and K from (4.12) to be operators from V × L 2 to H 1 as follows, Using these operators, we may write the general solutions φ(t, ξ), ψ(t, ξ) as An argument analogous to [ and Finally turning back to the functions we are interested in, we note that, since our interpolants R Δ and Q Δ are solutions of (6.10), they can be written as These should then be compared to R and Q for the exact solution, which now can be written as Then, we write and it follows from the Lipschitz property that for constants C Q,1 , C Q,2 , and an analogous estimate holds for R Δ (t, ·) − R(t, ·) H 1 .
From the above estimates, the obvious inequality f ξ L 2 f H 1 , and f V we may add the equations in (6.9) to obtain where we have used X := sup 0 t T X(t) E and sup 0 t T X Δ E C( Y Δξ ) are bounded by constants depending on T and the E-norm of their initial data. In particular, by theorem 5.7 we know Y Δξ is bounded by a constant depending only on T, H ∞ , ζ(0) ∞ , and ρ ∞ . Grönwall's inequality then yields the estimate Combining this estimate with (6.5), we obtain the desired result.
Since convergence in Lagrangian coordinates implies convergence in the corresponding Eulerian coordinates, see [28] for details, this shows that interpolated solutions of the discrete two-component Camassa-Holm system can be used to obtain conservative solutions of the 2CH system (1.2). In particular, as conservative solutions of (1.1) are unique according to [4], our discretization of the CH equation corresponds to the unique conservative solution of the CH equation.

Construction of the initial data
In this section we consider initial data for the continuous system given by u 0 ∈ H 1 ,ρ 0 = ρ 0 − ρ ∞ ∈ L 2 , and a measure μ 0 which corresponds to the energy distribution, see [31]. To ease notation we omit the subscript 0 and the dependence on t for the rest of this section, as we are always considering t = 0. The absolutely continuous part of the measure μ satisfies and may in general contain singular parts. Here we will restrict ourselves to the case where the singular part is purely atomic, and construct corresponding initial data for the discrete scheme. The ability to handle singular initial data was one of the motivations for the effort put into section 3 to allow for D + y i = 0. From [31,Thm. 4.9], we know the functions (y, U, H, r) defined as give us the initial data for the equivalent system (4.7) which provides us the global conservative solutions of (1.2) with initial data (u, μ). We define the discrete initial data y j = y(ξ j ) and U j = U(y j ). For the L 2 -functionr we definer The discrete identity (5.14) is essential to obtain global existence of solution to the semidiscrete system. The identity reflects the strong connection between the energy variable H j and the other variables. To fulfill (5.14), we set h j as follows: if D + y j > 0, we define h j 0 such that it satisfies (5.14), that is 2h j = U 2 j D + y j + (D + U j ) 2 D + y j +r 2 j D + y j (7.2) and if D + y j = 0 we set h j = 1 2 . Then we define H j = Δξ j−1 m=−∞ h m to ensure D + H j = h j . Note that in the norms below we will use U to denote both the continuous-case function U(ξ) and the discrete function {U j } j∈Z . However, the norm used will indicate which of them we are considering: p and V Δξ are used for discrete functions, and L p and V are used for continuouscase functions. Theorem 7.1. We consider the initial data of the two-component Camassa-Holm system (1.2) given by u 0 ∈ H 1 , ρ 0 such that ρ 0 − ρ ∞ =:ρ 0 ∈ L 2 for some ρ ∞ 0, and a positive finite Radon measure μ 0 whose absolutely continuous part satisfies μ 0,ac = (u 2 0 + u 2 0,x +ρ 2 0 )dx, while its singular part may be an atomic measure (the singular continuous part is zero). By definition, the global conservative solution of 2CH is obtained by solving (4.7) for the initial datum X 0 ∈ E constructed from (u 0 , ρ 0 , μ 0 ), where X 0 is given in (7.1). For this X 0 , we can construct sequences of initial data for the semi-discrete scheme, X 0,n = (ζ 0,n , U 0,n , H 0,n ,r 0,n ) ∈ E Δξ such that each element of the sequence belongs to the set B defined in definition 5.4 and the interpolation sequence defined in (6.2) converges to X 0 in E.
Proof. Let us start by verifying that these initial data satisfy properties (a)-(d) in definition 5.4. Clearly, from (7.1a) it follows that D + y j 0, and so the construction in Section 3 gives fundamental solutions satisfying property (a). Properties (c) and (d) have already been satisfied through our definition of h j . To verify (b), we need to show that the discrete initial data are uniformly bounded. Following [31,36] we have |y(ξ) − ξ| μ(R), and since the total energy μ(R) is bounded, we have y − Id L ∞ μ(R). Since y j = y(ξ j ), this carries directly over to our setting, |y j − ξ j | μ(R), meaning ζ ∞ μ(R). Moreover, in the aforementioned works, the authors prove that ξ → y(ξ) is Lipschitz with Lipschitz constant 1, which yields |y(ξ j+1 ) − y(ξ j )| |ξ j+1 − ξ j | = Δξ =⇒ |D + y j | 1.
Applying once more the continuous-case inequality (ρ • y)y ξ 1, the above estimate yields r j D + y j . From the preceding estimates, D + y j 1, and (7.2) we find Hence, h ∈ ∞ . We have ξ jr 2 (η)dη, and multiplying with Δξ and summing over j we obtain r 2 2 r 2 L 2 . Then it remains to check that H(0) ∈ V Δξ , and from (5.14) we estimate 2h j = U 2 j D + y j + (D + U j ) 2 +r 2 j − 2h j D + ζ j U 2 j + (D + U j ) 2 +r 2 j + h j + h j |D + ζ j | 2 . (7.5) Now, summing over j we find h 1 U 2 h 1 + r 2 + h ∞ D + ζ 2 2 , where the right-hand side is bounded by our previous estimates. Since h j > 0, it follows from our definition of H j that H j < H j+1 and H j < h 1 , which yields H ∞ = h 1 . Finally, we have h 2 h ∞ h 1 , so H ∈ V Δξ . Thus, we have proved that X j belongs to B.
Let us now prove that the interpolants for these initial data defined by (6.2) converge to the continuous initial data in E-norm. We start with ζ in L ∞ -norm, Next, we consider the L 2 -norm of U, Then, for the L 2 -norm of U ξ , we have where in the final limit we use [7, lemma 4.3]. A completely analogous estimate holds for the convergence of ζ Δ,ξ in L 2 . Consideringr Δ we find (r(η) −r(ξ)) 2 dη dξ, and following the proof for U ξ we find that this also converges. It remains to prove H Δ → H in V. We shall first prove that h Δ converges to h in L 1 , and we do it as follows. For a given n ∈ {1, 2, . . .}, we consider the partition of R defined by the points ξ i,n = i2 −n , which corresponds to Δξ = 2 −n . In this way, each partition is a subdivision of a coarser partition. We denote h Δ by h Δ n and similarly for all the other variables. We consider the sets B = {ξ ∈ R s.t. y ξ (ξ) = 0}, B n = {ξ ∈ R : there exists i ∈ Z, s.t. ξ ∈ (ξ i,n , ξ i+1,n ) and y(ξ i+1,n ) = y(ξ i,n )}.
Let us also define B o as the union B o = ∪ n 0 B n . Since y is increasing, we have B n ⊂ B. Moreover, as partitions for larger n are obtained by further subdivision, we have B n ⊂ B n+1 . Let L be the set of Lebesgue points for y ξ . We know that the set of Lebesgue points have full measure, that is m(L c ) = 0. We have, for some i, |y Δ n ,ξ (ξ) − y ξ (ξ)| 1 Δξ ξ i+1,n ξ i,n |y ξ (ξ) − y ξ (η)| dη 1 Δξ ξ+Δξ ξ−Δξ |y ξ (ξ) − y ξ (η)| dη, which tends to zero for any Lebesgue point ξ ∈ L. We consider a measure μ such that the singular part does not contain any singular continuous part, that is of the form for a sequence a j 0 such that a j 1 < ∞. When μ s takes this form, the set B can be written as for some values γ i ∈ R, for which we do not need explicit expressions in this proof. In this case, we have Indeed, this is a consequence of m((γ i , γ i + a i ) ∩ (B o ) c ) = 0, which can be proved as follows.