Rigorous asymptotics of a KdV soliton gas

We analytically study the long time and large space asymptotics of a new broad class of solutions of the KdV equation introduced by Dyachenko, Zakharov, and Zakharov. These solutions are characterized by a Riemann--Hilbert problem which we show arises as the limit $N\nearrow + \infty$ of a gas of $N$-solitons. We establish an asymptotic description for large times that is valid over the entire spatial domain, in terms of Jacobi elliptic functions.

can be presented as the existence of a simultaneous solution to the pair of equations −ψ xx + uψ = Eψ , (1.2) where E is the spectral parameter and ψ = ψ(x, t). The Lax pair formulation yields a complete solution procedure for the initial value problem for (1.1) via the inverse scattering transform in the case of rapidly decaying or step-like initial data, and has led to large and ever-growing collection of results concerning the analysis of the initial value problem in many different asymptotic regimes, including the behaviour in the small dispersion limit, as well as a complete description of the long-time behaviour for fairly general decaying or step-like initial conditions. In the case of periodic boundary conditions as well, there have been many works that are aimed at understanding the behaviour of solutions as well as the geometry of the space of solutions. These works have all been driven by the physical origins of the KdV equation as a basic model for one-dimensional wave motion of the interface between air and water, and in particular the discovery of the soliton. The soliton is a rapidly decreasing travelling wave solution of the KdV equation, namely a solution of the form u(x, t) = f (x − vt) and takes the form where E = −η 2 with E the energy parameter of Schrödinger equation in the Lax pair (1.2). The periodic travelling wave that can be obtained by direct integration of the KdV equation and takes the form where dn(z m) is the Jacobi elliptic function of modulus m = β2−β3 β1−β3 and β 1 > β 2 > β 3 . In both formulas x 0 is an arbitrary phase. With the above potential (1.5), the Schrödinger equation (1.2) coincides with the Lamé equation and the stability zones (or Bloch spectrum) of the potential are [β 3 , β 2 ] ∪ [β 1 , +∞).
Of the highest importance for applications to the theory of water waves was the discovery of families of explicit more complex solutions, such as N-soliton solutions when the Schrödinger equation in (1.2) has N simple eigenvalues, or a N -gap solution when there are N + 1 disjoint stability zones of the corresponding Schrödinger equation or solutions that connect to Painlevé transcendents.
Since the early days of integrable nonlinear PDEs, researchers have considered the notion of a soliton gas (see [Zak09], and references contained therein). The quest is for an understanding of the properties of an interacting ensemble of many solitons, ultimately in the presence of randomness. However, even in the absence of randomness, the dynamics of a large collection of solitons is only understood with mathematical precision in a few specific settings (the smalldispersion limit of the KdV equation, as considered in the works of Lax and Levermore [LL83a,LL83b,LL83c] could be interpreted as a highly concentrated soliton gas, with a smooth and rapidly decaying function being represented as an infinite accumulation of solitons).
Within integrable turbulence, the interest is in the computation of statistical quantities describing evolving random configurations of solitons. In [DP14] and [SP16] the authors used computational methods to approximate such statistical quantities via the Monte-Carlo method, and presented a formal derivation of evolution equations for the first four statistical moments of the solution. In another direction [Zak71,EK05] the interest is in computing a kinetic equation describing the evolution of the spectral distribution functions. This has been extended to similar formal considerations based on properties of fundamental solutions in the periodic setting, as opposed to solitonic gasses [El16,EKPZ11].
1.1. The soliton gas. Towards the goal of discovering new, broad families of solutions to integrable nonlinear PDEs, the "dressing method" as developed by Zakharov and Manakov [ZM85] has yielded some interesting new results in [DZZ16]. In that paper, the authors show how the dressing method can be used to produce a new family of solutions they refer to as "primitive potentials" which although are not random, can be naturally interpreted as a soliton gas. Cutting to the chase, the authors derive a Riemann-Hilbert problem which is to determine a vector Ξ = Ξ 1 Ξ 2 T satisfying a normalization condition at ∞, and the jump relations where the jump matrix J(λ) is given by J(λ) = 1 1 + r 1 (λ)r 2 (λ) 1 − r 1 (λ)r 2 (λ) 2ir 1 (λ)e −2λx 2ir 2 (λ)e 2λx 1 − r 1 (λ)r 2 (λ) . (1.7) The parameters η 1 and η 2 are taken to be real with 0 < η 1 < η 2 , and the orientation of the jumps is the following: both the intervals (iη 1 , iη 2 ) and (−iη 2 , −iη 1 ) are oriented downwards. The reflection coefficients r 1 (λ) = r 1 (λ; t) and r 2 (λ) = r 2 (λ; t) evolve in time according to r 1 (λ; t) = r 1 (λ; 0)e (8λ 3 −12λ)t , r 2 (λ; t) = r 2 (λ; 0)e −(8λ 3 −12λ)t . (1.8) The authors consider a number of different settings, and use a combination of analytical and computational methods to provide a description of the solutions of the KdV equation determined by this Riemann-Hilbert problem. In the case that r 2 ≡ 0, the potential is exponentially decaying as x ↗ +∞. But the behavior as x grows in the other direction (as well as the the asymptotic behavior for x large in the case that both reflection coefficients are nontrivial) was mentioned as a challenging problem for both analysis and computation. The configuration of solitons considered in [DZZ16] is somewhat different than the solitonic gas configurations considered in [DP14] and [SP16], where they considered a large number of solitons that were spaced quite far apart from each other at t = 0. In other words, they considered a dilute gas of solitons that had enough space between them to evolve as isolated solitons until they interact, usually in a pair-wise fashion. In contrast, the soliton gas considered in [DZZ16] (and considered here as well) is a configuration that cannot be viewed as a collection of isolated solitons. Indeed, as we show, they are overlapping to the extent that, even at t = 0, they effectively behave as a genus-1 dispersive wave, and this poses a conceptual challenge in that finite-soliton interactions may or may not be relevant.
1.2. Statement of the results. In Section 2 we consider a sequence of Riemann-Hilbert problems, indexed by N , for a pure N -soliton solution, with spectrum confined to the intervals (−iη 2 , −iη 1 ) ∪ (iη 1 , iη 2 ) for some η 2 > η 1 > 0 and show that for this sequence, as N ↗ +∞, the solution of the Riemann-Hilbert problem converges to the solution of the Riemann-Hilbert problem studied in [DZZ16], for the case r 2 ≡ 0.
Then in Section 3 (Theorem 3.6) we establish that the potential u(x) determined by this Riemann-Hilbert problem coincides with the Lamé potential as x ↘ −∞: The function dn (z m ) is the Jacobi elliptic function of modulus m = η 1 η 2 . It is periodic with period 2K(m), and satisfies dn (0 m ) = 1 and dn (K(m) m ) = √ 1 − m 2 . Here K(m) is the complete elliptic integral of the first kind of modulus m: K(m) = ∫ π 2 0 dϑ √ 1−m 2 sin 2 ϑ . Therefore the function (1.9) is periodic in x with period 2K(m) η 2 .The minimum amplitude of the oscillations is −η 2 2 − η 2 1 and the maximum amplitude is η 2 1 − η 2 2 so that the amplitude of the oscillations is 2η 2 1 . The average value of u(x) over an oscillation is equal to where E(m) is the complete elliptic integral of the second kind: E(m) = ∫ π 2 0 √ 1 − m 2 sin 2 ϑ dϑ. The phase φ in formula (1.9) depends on the coefficient r 1 (λ) that characterizes the continuum limit of the norming constants of the soliton gas and it is equal to (1.10) A technical issue arises in the control of the error, and when η2 2K(m) (x + φ) = 1 2 + n, n ∈ Z, one requires a different local parametrix in the Riemann-Hilbert analysis or one needs consider, for this particular value of the parameter x, a vector-valued Riemann-Hilbert problem for the error, whereas for all other values of x we essentially use an analysis based on matrix-valued Riemann-Hilbert problems. While the error term is expected to be uniform, in our analysis we explicitly assume that η2 2K(m) (x + φ) ≠ 1 2 + n, n ∈ Z. Finally in Sections 4 -6 we provide a global long-time asymptotic description of the solution u(x, t) to the KdV equation with this initial data u(x). The asymptotic behaviour depends on the quantity ξ = x 4t. There are three main regions, (1) a constant region, (2) a region where the solution is approximated by a periodic traveling wave with constant coefficients specified by the spectral data, (3) a region where the solution is approximated by a periodic travelling wave with modulated coefficients see Figure 1. More precisely: (1) for fixed ξ > η 2 2 , there is a positive constant C = C(ξ) so that u(x, t) = O e −Ct .
The equation (1.14) was used by Gurevich and Pitaevskii [GP73] to describe the modulation of the travelling wave that is formed in the solution of the KdV equation with a step initial data u(x, 0) = −η 2 2 for x < 0 and u(x, 0) = 0 for x > 0. Such a modulated travelling wave is also called a dispersive shock wave. The rigorous analysis of the dispersive shock wave emerging from step-like initial data problem was obtained via inverse scattering in [Hru76] and more recently via Riemann-Hilbert methods in [EGKT13].
We observe that in this paper, if η 1 is set to 0, then m = 0 and lim m→0 dn(z m) = 1 so that our asymptotic description in (1.9) yields a potential that converges to −η 2 2 as x ↘ −∞. The behavior for t large agrees (globally) with the dispersive shock wave emerging from a step like initial data with reflection coefficient equal to zero on the real axis. In other words, in the case that η 1 = 0, the "primitive potential" identified in [DZZ16] is in the class of step-like potentials.

Soliton gas as limit of N solitons as N ↗ +∞
The Riemann-Hilbert problem for a pure N -soliton solution (see for example [GT09]) is described as follows: find a 2-dimensional row vector M such that where c j ∈ iR; where M 1 (λ) is the first entry of the vector M (λ). In particular, for a one soliton potential, namely N = 1, one recovers the expression (1.4) where the shift x 0 is given by We are interested in the limit as N ↗ +∞, under the additional assumptions are sampled from a density function (λ) so that ∫ −iλj η1 (η)dη = j N , for j = 1, . . . , N . (ii) The coefficients {c j } N j=1 are purely imaginary (in fact c j ∈ iR + ) and are assumed to be discretizations of a given function: where r 1 (λ) is an analytic function for λ near the intervals (iη 1 , iη 2 ) and (−iη 2 , −iη 1 ), with the symmetry r 1 (λ) = r 1 (λ), and is further assumed to be a real valued and non-vanishing function of λ for λ ∈ [iη 1 , iη 2 ]. In the regime x ↗ +∞, it is easy to notice that all residue conditions contain only exponentially small terms and therefore, by a small norm argument, the potential is exponentially small.
On the other hand, for x ↘ −∞ all of those terms are exponentially large. To show that the solution is also exponentially small in this latter case, one may reverse the triangularity of the residue conditions, by defining (2.4) Now the residue conditions are while the potential u(x) is still extracted from A via the same calculation: The quantity e −2iλj x now decays exponentially as x ↘ −∞, and this implies (again by a standard small-norm argument) exponential decay of the potential u(x) for x ↘ −∞. On the other hand, the product term is exponentially large in N : for each j = 1, . . . , N (2.7) for some C ∈ R + . Therefore this exponential decay does not set in until x is rather large. Indeed, in order for the residue conditions to all be exponentially small, it must be that x ≪ −CN . In other words, the N-soliton solution that we are considering has very broad support, and in the large-N limit, it is not exponentially decaying for x ↘ −∞.
We will show here how to derive the Riemann-Hilbert problem for a soliton gas with one reflection coefficient r 1 (as described in [DZZ16]) from a meromorphic Riemann-Hilbert problem for N solitons in the limit as N ↗ +∞. First, we remove the poles by defining within a closed curve γ + encircling the poles counterclockwise in the upper half plane C + , and within a closed curve γ − surrounding the poles clockwise in the lower half plane C − . Outside these two sets, we take Z(λ) = M (λ).
Then the jumps are where, for λ ∈ γ + or γ − , the boundary values Z + (λ) are taken from the left side of the contour as one traverses it according to its orientation, and the boundary values Z − (λ) are taken from the left. The quantity Z(λ) is normalized so that Z(λ) = 1 1 + O λ −1 as λ → ∞.
We assume now that in the limit as the number of poles goes to infinity, the poles are distributed according to some distribution (λ) with density compactly supported in (iη 1 , iη 2 ) (and extended by symmetry on the corresponding interval in the lower half plane).
For the sake of simplicity, we can assume that the N poles are equally spaced along (iη 1 , iη 2 ) with distance between two poles equal to ∆λ = η2−η1 N and with (atomic) density: for some normalization constant Z N .
Remark 2.1. In the case where the poles are distributed according to a more general measure (λ), the steps to follow are very similar. The entries of the jump matrices will carry the density function along, which can be eventually incorporated in the definition of the reflection coefficient r 1 .
As the number of poles increases within the support of the measure, the following result holds.
Proof. Using (2.3), the expressions in the jumps can be rewritten as (2.14) The convergence follows from the convergence of the Riemann sum to the Riemann-Stieltjes integral for x ∈ K any compact subset of R.
Thanks to the proposition above and a small norm argument, we arrive at a limiting Riemann-Hilbert problem (which we still call Z with abuse of notation) At this point it is important to point out to the reader that the contour (iη 1 , iη 2 ) and (−iη 2 , −iη 1 ) are both oriented upwards. Next, we define within the loop γ + , and within the loop γ − . Outside these two curves, we define X(λ) = Z(λ).
The jumps across the curves are no longer present, but there are jumps across (iη 1 , iη 2 ) and (−iη 2 , −iη 1 ) because the integrals have jumps across those intervals. Using Sokhotski-Plemelj formula, we arrive at a Riemann-Hilbert problem for X This Riemann-Hilbert problem is equivalent to the one described in [DZZ16] with r 2 = 0, up to a transposition (X is a row vector here, while the solution of the Riemann-Hilbert problem in [DZZ16] is a column vector), and using the symmetry that r 1 (λ) = r 1 (λ) for λ ∈ (−iη 1 , −iη 2 ). We note that there is a sign discrepancy between this Riemann-Hilbert problem and the one appearing in [DZZ16], which is resolved by a careful interpretation of the sign conventions therein.
3. Behaviour of the potential u(x, 0) as x ↘ −∞ We consider a soliton gas Riemann-Hilbert problem as in (2.19) -(2.20) with 0 < η 1 < η 2 , and reflection coefficient r 1 (λ) defined on (iη 1 , iη 2 ) such that it has an analytic extension to a neighbourhood of this interval. Furthermore we assume that r 1 (−λ) = r 1 (λ) on the imaginary axis. We set Σ 1 = (η 1 , η 2 ) and Σ 2 = (−η 2 , −η 1 ). The vector valued function X that will determine the KdV potential u(x) is the solution to the following Riemann-Hilbert problem: As explained in [DZZ16], we can recover the potential u(x) of the Schrödinger operator via the formula where X 1 (λ; x) is the first component of the solution vector X. We first perform a rotation of the problem in order to place the jumps on the real line, by setting the Riemann-Hilbert problem for Y reads as follows The contour setting is shown in Figure 2. We can recover u(x) from 3.1. Large x asymptotic. Introduce the following new vector function where g(λ) and f (λ) are scalar functions to be determined below. We require that where Ω is a constant independent of x and needs to be determined and (3.10) In order to solve the scalar Riemann-Hilbert problem (3.7) -(3.9) for g we observe that From the above, we can write g ′ (λ) as where is real and positive on (η 2 , +∞) with branch cuts on the contours Σ 1 and Σ 2 and κ is a constant to be determined. By integration we obtain The condition (3.7) implies that η1 −η1 ζ 2 + κ R(ζ) dζ = 0 and the condition (3.8) implies that This gives is the complete elliptic integral of the first kind with modulus m = η 1 η 2 and (3.20) In order to solve the Riemann-Hilbert problem for T (λ) we wish to obtain a constant jump matrix J T . For this purpose we make the following ansatz on the function f It is easy to check that the function f (λ) is given by (3.26) The constraint (3.25) gives ∆ equal to where in the last equality in (3.27) we use the fact that r(−λ) = r(λ). We remind the reader that we are assuming the function r to be real positive and non-vanishing on Σ 1 and Σ 2 .
3.2. Opening lenses. We start by defining the analytic continuationr(λ) of the function r(λ) off the interval (−η 2 , −η 1 ) ∪ (η 1 , η 2 ) with the requirement that We can factor the jump matrix J T on Σ 1 as follows We can now proceed with "opening lenses". We define a new vector function S as follows where the Riemann-Hilbert problem satisfied by S is depicted in Figure 3. In order to proceed we need the following lemma Lemma 3.1. The following inequalities are satisfied where C 1 and C 2 are the contours defining the lenses as shown in Figure 3.
Lemma 3.1 guarantees that the off diagonal entries of the jump matrices along the upper and lower lenses are exponentially small in the regime as x ↘ −∞, therefore those jump matrices are asymptotically close to the identity outside a small neighbourhoods of ±η 1 and ±η 2 . We are left with the model problem The Riemann-Hilbert problem for S (∞) appeared already in the long time asymptotic for KdV with step-like initial data [EGKT13]. Below we follow the lines in [EGKT13] to obtain the solution.
3.3. The global parametrix S (∞) . For solving the Riemann-Hilbert problem (3.32) and (3.33) we introduce a two-sheeted Riemann surface X of genus 1 associated to the multivalued function R(λ), namely The first sheet of the surface is identified with the sheet where R(λ) is real and positive for λ ∈ (η 2 , +∞). We introduce a canonical homology basis with the B cycle encircling Σ 1 clockwise on the first sheet and the A cycle going from Σ 2 to Σ 1 on the first sheet and coming back to Σ 2 on the second sheet. The points at infinity on the surface are denoted by ∞ ± where ∞ + is on the first sheet and ∞ − on the second sheet of X. See Figure 4. We introduce the holomorphic differential We also have Next we introduce the Jacobi elliptic function which is an even function of z and satisfies the periodicity conditions We also recall that the Jacobi elliptic function with modulo τ vanishes on the half period τ 2 + 1 2 . Finally, we define the integral w(λ) = λ η2 ω and we observe that We introduce the following functions , and we observe that Figure 4. Construction of the genus-1 Riemann surface X and its basis of cycles.
It follows that the functions ψ 1 and ψ 2 have simple poles at λ = ±η 1 . Furthermore the following jump relations are satisfied: while for λ ∈ (−η 1 , η 1 ) (3.43) Next we introduce the quantity We are now ready to construct the solution of the Riemann-Hilbert problem (3.32) -(3.33).
Proof. We observe that S (∞) (λ) has at most fourth root singularities at the branch points and it is regular every where else on the complex plane. Because of (3.36) and (3.37) we have S (∞) (∞) = 1 1 , namely the condition (3.33) is satisfied. Combining (3.42), (3.43) and (3.44), we conclude that the jump conditions (3.32) are satisfied.
This vector solution provides the asymptotic behavior of the solution S to Riemann-Hilbert problem depicted in Figure 3, for all z bounded away from the endpoints. However, in order to prove this, we need to construct a matrix solution to this Riemann-Hilbert problem, which we call P (∞) . This will be accomplished in the next two subsections, by creating a second, independent, vector solution.
3.4. A second vector solution. In order to construct a second vector solution, note that for any complex numbers a and b, we have that the following relations are satisfied: and In addition, for any constants D, we have that the vector (3.50) We choose D so that the quantity ϑ 3 (w(λ) + D; τ ) has a zero at the point ∞ − on the second sheet of X namely In this way the quantity ϑ 3 (−w(λ) + D; τ ) has a simple zero at ∞ + on the first sheet of the surface X. Next we choose a = b = 1 so that the equation has a double zero at ∞ − and the equation has a double zero at ∞ + . Therefore the vector H(λ) = H 1 H 2 : , has no poles and at most fourth root singularity at the points λ = ±η 1 and λ = ±η 2 and satisfies the boundary value relations (3.54) Computing the behavior for λ → ∞, we find Such solution is well defined and linearly independent from the solution S (∞) defined in (3.45) when Remark 3.3. When xΩ + ∆ = πi + 2πi mod n, n ∈ Z, a technical / analytical issue arises concerning control of the error in the Riemann-Hilbert analysis. First, the solution S (∞) vanishes at λ = 0. Indeed, at this special value, we have and, taking λ = 0, we see that both entries of S (∞) vanish. Moreover, the second vector solution H(λ) computed above vanishes identically as λ → ∞. Indeed, at the special value xΩ + ∆ = πi, we have

57)
and clearly both quantities converge to 0 as λ → ∞. We see that for this special value of the parameters, the model Riemann-Hilbert problem defining S (∞) is somewhat special in that there is a well-behaved vector-valued solution that vanishes identically as λ → ∞, and the question of the existence of a suitably well-behaved matrix-valued solution is nontrivial. We can construct a matrix-valued solution in this case, as follows. We takê It is straighforward to verify that detĤ = 1. (3.59) While this represents a reasonable solution to the outer model Riemann-Hilbert problem for xΩ+∆ 2πi = 2n+1 2 , n ∈ Z , it raises the issue of the construction of a local parametrix near λ = η 2 which is new, and as far as we know, open, and will not be addressed here.
On the other hand, the global parametrix P (∞) is a good approximation of the solution S to the Riemann-Hilbert problem away from the endpoints λ = ±η 2 , ±η 1 , where P (∞) exhibits a fourth-root singularity. So, we need to introduce four local parametrices P (±ηj ) (j = 1, 2) in a suitable neighbourhood of each endpoint.
We show here the construction of a (matrix) local parametrix P (η2) around λ = η 2 . The constructions of the other local parametrices near λ = ±η 1 , −η 2 follow a similar argument.
In the above formulae I 0 (ζ), K 0 (ζ) are the modified Bessel functions of first and second kind, respectively, and H (j) (ζ) the Hankel functions.
In conclusion the local parametrix around the endpoint λ = η 2 is where A is a prefactor that is determined by imposing that (3.74) Therefore, we set (3.75) By construction, A is well-defined and analytic in a neighbourhood of η 2 , minus the cut (−∞, η 2 ]; additionally, it is easy to see that A is invertible (det A(λ) = 1).
Lemma 3.5. A(λ) is analytic everywhere in the neighbourhood B ρ of η 2 .
Proof. To prove the statement, one needs to check that A has no jumps across the interval [η 1 , η 2 ] ∩ B ρ and that it has at most a removable singularity at λ = η 2 . Careful computation (using the definition (3.75), with tears) reveals that A + (λ) = A − (λ) on (−∞, η 2 ) ∩ B ρ . Next, we notice that ζ(λ) has a simple zero at η 2 by construction, thus ζ(λ) 1 4 σ3 has at most a fourth-root singularity at the point λ = η 2 . Also the global parametrix P (∞) (λ) has at most a fourthroot singularity near η 2 and consequently all the entries of A(λ) have at most a square root singularity at λ = η 2 .
On the other hand A(λ) is analytic in B ρ {η 2 }, therefore the point λ = η 2 is a removable singularity and A(λ) is indeed analytic everywhere in B ρ . 3.7. Small norm argument and determination of u(x, 0) for large negative x. Consider the following "remainder" Riemann-Hilbert problem: where U is the (matrix) ensemble of the global parametrix P (∞) and the four local parametrices P (±ηj ) , j = 1, 2; thanks to the previous arguments, the vector-valued function E satisfies a Riemann-Hilbert problem with jumps that are asymptotically close to the identity matrix: on the upper and lower lenses on the circles around the endpoints (3.77) and E(λ) = 1 1 + O λ −1 as λ → ∞. Therefore, by a standard small norm argument (see, for example [Its11, Section 5.1.3]) the solution E behaves as follows: E(λ) = 1 1 + O x −1 as x ↘ −∞. We note in passing that the construction of a matrix-valued global approximation is very useful, in that we arrive directly at a small-norm Riemann-Hilbert problem. In addition, we note that the error term has not been established for x in a vicinity of those values where xΩ+∆ 2πi = 2n+1 2 , n ∈ Z . Keeping into account all the transformations we performed, we are now able to explicitly solve the original Riemann-Hilbert problem Y in the large negative x regime: (3.78) and in place of U (λ) we use the global or local parametrix. We recall that the potential u(x) can be calculated from the solution Y (λ) as where Y 1 (λ; x) is the first entry of the vector Y .
Theorem 3.6. In the regime x ↘ −∞, with xΩ+∆ 2πi ≠ 2n+1 2 , n ∈ Z , the potential u(x) has the following asymptotic behaviour where E(m) and K(m) are the complete elliptic integrals of the first and second kind respectively with modulus m = η 1 η 2 , φ is given by The formula (3.80) can be written in the equivalent form where dn (z m ) is the Jacobi elliptic function of modulus m.
Proof. We are interested in the first entry of the vector Y (λ) (for λ large), and we have, from (3.64), From the expression of the g function (3.16), we have From the formula of f (λ) in (3.26) we have where we have used the relations and the parity of the function ϑ 3 (z + τ 4 ; τ 2 ). From the above expansions, using the explicit expression of ∆ in (3.27), and the parity of ϑ 3 (z; 2τ ) we obtain the expression of u(x) in (3.79). In order to obtain the expression (3.82) we need the following identity ([Law89, pg.59]) where dn (z m ) is the Jacobi elliptic function of modulus m and period 2K(m) and we recall that 2τ = iK(m ′ ) K(m). Then we can write so that the expression for u(x) in (3.80) can be written in the form (3.82).
Remark 3.7. While the above theorem requires xΩ+∆ 2πi ≠ 2n+1 2 , n ∈ Z , the reason is because we have not constructed the requisite new parametrix near η 2 in this case. It is expected that the leading order behavior should not change for x near such values, whereas the error term may or may not be altered.

Behaviour of the potential u(x, t) as t ↗ +∞
Letting the potential u(x, t) evolve in time according to the KdV equation, the reflection coefficient evolves as r 1 (λ; t) = r 1 (λ)e −8λ 3 t . This will lead to the study of a Riemann-Hilbert problem Y for the soliton gas described as follows We are interested in the asymptotic behaviour of Y (λ) in the long-time regime (t ↗ +∞).
The phase appearing in the exponents in the jumps shows different sign depending on the value of the quantity It is clear that in the case ξ > η 2 2 , the phases in the jumps are exponentially decaying in the regime t ↗ +∞, therefore by a straightforward small norm argument we conclude as t ↗ +∞ , (4.4) and the potential u(x, t) becomes trivial. The more interesting case ξ ≤ η 2 2 will be studied below. It will be clear that we will observe the presence of critical value ξ crit at which a phase transition will occur when passing from ξ > ξ crit (the "super-critical" case) to ξ ≤ ξ crit (the "sub-critical" case). In the first case the asymptotic description gives an asymptotic solution that is a modulated travelling wave (the wave parameters are changing slowly in time), while in the sub-critical case, the asymptotic solution is a travelling wave.

Super-critical case: the α-dependency
We first consider the case where the value of ξ crit ∈ R will be defined in (5.18).
Before continuing our analysis we want to comment on equation (5.16). We can rewrite it in the form This relation describes the modulation of the parameter α as a function of ξ. The quantity η 2 2 W (m α ) was derived by Whitham in his modulation theory of the traveling wave solution of the KdV equation [Whi74]. In the general theory there are three parameters involved, while in our case, two parameters are fixed, one being zero and the other one η 2 . This specific case gives a self-similar solution to the Whitham equations. This solution was derived and used by Gurevich-Pitaevskii [GP73] to describe the modulation of the travelling wave that is formed in the solution of the KdV equation with step initial data u(x) = −η 2 2 for x < 0 and u(x) = 0 for x > 0 and was called a dispersive shock wave in analogy with the shock wave that is formed in the solution of the Hopf equation u t + 6uu x = 0 for step initial data.
Using the expansion of the elliptic functions one has , as m α → 1 , The Whitham equations are strictly hyperbolic ( [Lev88]), so that ∂ ∂α W (m α ) > 0 for 0 < α < η 2 . Hence by the implicit function theorem, the equation (5.17) defines α as a monotone increasing function of ξ for ξ ∈ [ξ crit , η 2 2 ] where ξ crit is given by Then, clearly ξ crit > − 3η 2 2 2 . From g ′ (λ), we also have a representation of g(λ): This, together with (5.5), yields the formulã For future use we will need the x derivatives of tg(λ) and tΩ. Before calculating them, let us observe thatΩ which gives, using the Riemann bilinear relations, (5.21) Lemma 5.1. The following identities are satisfied Proof. We observe that g ′ (λ)dλ defined in (5.11) is a meromorphic one-form on the Riemann surface X α defined as X α = (η, λ) ∈ C 2 η 2 = R 2 α (λ) = (λ 2 − α 2 )(λ 2 − η 2 2 ) . We define the homology basis on X α in the following way: the B cycle encircles the cut [α, η 2 ] clockwise and the A cycle starts on the cut [−η 2 , −α] on the upper semi-plane, goes to the cut [α, η 2 ] and then goes back to [−η 2 , −α] on the second sheet of X α . Then we have (5.24) Regarding the first relation in (5.22) we have dλ vanishes because it is a holomorphic one-form, (no singularity at ±α or infinity) and it is normalized to zero on the A cycle because of (5.24) and therefore it is identically zero [Kri88] (see also [Gra02]). Another alternative proof is to calculate the derivative and use the explicit formulae of the constants c 1 and c 2 in (5.15). We conclude that Regarding the relation (5.23), by (5.22) and (5.24) we have As we did in Section 3, we choose the function f to simplify the jumps on Σ 1,α and Σ 2,α via It is easy to check that the function f (λ) is given by (5.32) The constraint (5.31) yields a formula for∆: where in the last relation in (5.33) we use the fact that r(−λ) = r(λ).
As a consequence, T satisfies the following Riemann-Hilbert problem: (5.34) (5.35) 5.1. Opening lenses. It is useful to provide representations of the entries appearing in the jump matrix for T (λ) in either Σ 1,α or Σ 2,α , that clearly demonstrate their analytic continuation off of these intervals, as was done in Section 3. The following formulae are valid on both intervals: The following formulae are valid on Σ 1,α : . (5.38) And the following ones are valid on Σ 2,α : As was done in Section 3, we can factor the jump matrix on Σ 1,α as follows These factorizations permit us to open lenses as shown in Figure 5.
Opening lenses: the grey entries in the jumps represents exponentially small quantities in the limit t ↗ +∞. The contours C 1 and C 2 are the lens boundaries.
Because of lemma 5.2, letting t ↗ +∞, the jump matrices (as depicted in the Figure 5) will converge to constant jumps exponentially fast outside neighbourhoods of ±α and ±η 2 . We then obtain the following model Riemann-Hilbert problem forS (∞) : The global parametrixP (∞) . Along the same lines as we did in Section 3, we construct a (matrix) model problem whose solution will yield a solution of the above (vector) Riemann-Hilbert problem. Since the solution of this model problem will be invertible, one is able to arrive at a small-norm Riemann-Hilbert problem for the error in the large-time regime, more directly than if one considers only vector Riemann-Hilbert problems. And, as with Section 3, we build this local parametrix under the assumption that tΩ+∆ 2πi ≠ 2n+1 2 , n ∈ Z . We want to determine the matrix valued functionP (∞) that is analytic in C (−η 2 , η 2 ) and satisfies the following Riemann-Hilbert problem The solution is obtained from the solution of P (∞) in (3.64) by replacing η 1 with α.
5.3. The local parametrix P (±α) . We will construct now a (matrix) local parametrix around the point λ = −α. The construction of the other local parametrix near the endpoint λ = α is analogous, while the construction of the local parametrices near λ = ±η 2 is the same one as in the Section 3.6 (again, under the assumption tΩ+∆ 2πi ≠ 2n+1 2 , n ∈ Z ).
In conclusion, our local parametrix is then defined as where A is an analytic prefactor whose expression is determined by imposing that In light of this asymptotic behaviour we set (5.62) By construction, A is defined and analytic in a neighbourhood of −α, minus the cuts (−∞, −α] ∪ [−α, +∞); moreover, A is invertible (det A(λ) ≡ 1).
Lemma 5.3. A(λ) is analytic in a whole neighbourhood of −α.
Proof. The proof entails verifying that A has no jumps across the interval (−α− , −α+ ) and that it has at most a removable singularity at λ = −α. We leave the verification that A + (λ) = A − (λ) across the interval (−α − , −α + ) to the reader, using the jump relations satisfied byP (∞) and the above definitions. The conformal map ζ(λ) has a simple zero at λ = −α (by construction), therefore ζ(λ) − 1 4 σ3 has at most a fourth-root singularity at −α. Similarly,P (∞) (λ) has a fourth-root singularity at −α, as well; therefore, all the entries of A(λ) have at most a square root singularity at λ = −α, and A(λ) is analytic in B ρ {−α}. The point λ = −α is a removable singularity. This implies that A(λ) is indeed analytic everywhere in B ρ .
Proof. Expanding each term of Y 1 (λ) in (5.4) in a neighbourhood of infinity gives the following. Regarding f (λ) defined in (5.32) we have .
Regarding e −tg(λ) we are interesting in the x derivative of this expression. Using (5.22) we have ∂ ∂x e −tg(λ) = − 1 λ where ′ stands for the derivative with respect to the argument of the theta-function. By (5.23) we have ∂ ∂xS where ′′ stands for second derivative with respect to the argument. Taking into account that for any smooth function F (α(ξ), η 2 ), ∂ ∂x F (α(ξ), η 2 ) = O(t −1 ) by (5.16), we can write the above expression in the form ∂ ∂xS Gathering the above expansions and using the explicit expression ofΩ and∆ in (5.21) and (5.33) respectively we obtain (5.68). Also in this case, using theorem 3.6 we can reduce the expression of u(x, t) to the form (5.69).
In order to show that the usual contour deformations can be carried out, as they were in Sections 3 and 4, we need to verify that the quantity Re 2g(λ) + 8λ 3 − 8ξ 2 λ is positive on the contour C 1 , and negative on the contour C 2 , where these contours are as shown in Figure 3.

Conclusions
In this paper we have considered the Riemann-Hilbert problem of [DZZ16] in the case of one non-trivial reflection coefficient. We have shown how this Riemann-Hilbert problem describes a soliton gas as the limit of a finite N -soliton configuration as N tends to +∞. Then we established rigorous asymptotics of the KdV potential in several different regimes. First, for the initial configuration, we studied the challenging behaviour as x ↘ −∞, and obtained a universal asymptotic description in terms of the periodic travelling wave solution of KdV. Then, we provided a complete analysis of the long-time behavior of the solution of the KdV equation determined by the Riemann-Hilbert problem of [DZZ16]. For large t, there are three fundamental spatial domains, in which the solution u(x, t) displays different asymptotic behaviour, either exponentially small, or in terms of the periodic travelling wave solution of KdV with slowly varying parameters, or the periodic travelling wave solution of KdV with fixed parameters.
Several challenges remain, like the asymptotic analysis when there are two nontrivial reflection coefficients or when the case where the spectral parameters of the soliton gas accumulates in disconnected components of the imaginary axis. Beyond these, it is enticing to consider the interaction of one large soliton with this gas like in [CDE16] or the interaction between two such soliton gases.