Rigorous Asymptotics of a KdV Soliton Gas

Abstract

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\rightarrow + \infty$ of a gas of N-solitons. We show that this gas of solitons in the limit $N\rightarrow \infty$ is slowly approaching a cnoidal wave solution for $x \rightarrow - \infty$ up to terms of order $\mathcal {O} (1/x)$, while approaching zero exponentially fast for $x\rightarrow +\infty$. We establish an asymptotic description of the gas of solitons for large times that is valid over the entire spatial domain, in terms of Jacobi elliptic functions.

1 Introduction

This paper concerns the concept of a gas of solitons for the Korteweg–de Vries (KdV) equation,

$$\begin{aligned} u_{t} - 6 u u_{x} + u_{xxx} = 0 \ . \end{aligned}$$

(1.1)

It is well known that this nonlinear partial differential equation is integrable, arising as the compatibility condition of a Lax pair of linear differential operators. The compatibility condition can be presented as the existence of a simultaneous solution to the pair of equations

$$\begin{aligned}&-\psi _{xx} + u \psi = E \psi \, , \end{aligned}$$

(1.2)

$$\begin{aligned}&\psi _{t} -4\psi _{xxx} + 6 u \psi _{x} + 3 u_{x} \psi = 0 \ , \end{aligned}$$

(1.3)

where E is the spectral parameter and $\psi =\psi (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 a 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, 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

$$\begin{aligned} u(x,t)=-2\eta ^2{\text {sech}}^2\left( 2\eta (x-4\eta ^2t-x_0) \right) \end{aligned}$$

(1.4)

where $E=-\eta ^2$ is the energy parameter of Schrödinger equation in the Lax pair (1.21.3). The periodic travelling wave that can be obtained by direct integration of the KdV equation takes the form

$$\begin{aligned} u(x,t)=\beta _1+\beta _2-\beta _3-2(\beta _1-\beta _3)\text{ dn}^2(\sqrt{\beta _1-\beta _3}(x+2(\beta _1+\beta _2+\beta _3)t+x_0)|m) \end{aligned}$$

(1.5)

where $\text{ dn }(z|m)$ is the Jacobi elliptic function of modulus $m^2=\frac{\beta _2-\beta _3}{\beta _1-\beta _3}$ and $\beta _1>\beta _2>\beta _3$. In both formulas $x_0$ is an arbitrary phase. Let us introduce the $\vartheta$ function

$$\begin{aligned} \vartheta _3 (z;\tau )= \sum _{n\in \mathbb {Z}} e^{2\pi i \, nz + \pi n^2 i \tau } \ , \qquad z \in \mathbb {C} \ ,\quad {\text {Im}}\tau >0. \end{aligned}$$

Using the standard relation between Jacobi elliptic functions and $\vartheta$-function (see eg. [Law89] pg. 45 exercise 16 and 3.5.5) we re-write (1.5) as

$$\begin{aligned} \begin{aligned} u(x,t)&= \bar{u}-2\frac{\partial ^2}{\partial x^2}\log \vartheta _3\left( \dfrac{\sqrt{\beta _1-\beta _3}}{2 K(m)}[x+2 t(\beta _1+\beta _2+\beta _3) +x_0]-\frac{1}{2};\hat{\tau }\right) \,\\ \bar{u}&=\beta _1+\beta _2-\beta _3-2(\beta _1-\beta _3)\dfrac{E(m)}{K(m)}\, \end{aligned} \end{aligned}$$

(1.6)

with $K(m)=\int _0^{\pi /2}\frac{\mathrm{d}\vartheta }{\sqrt{1-m^2\sin ^2\vartheta }}$ and $E(m)=\int _0^{\pi /2}\mathrm{d}\vartheta \sqrt{1-m^2\sin ^2\vartheta }$, the complete elliptic integrals of the first and second kind respectively, $\hat{\tau }=iK'(m)/K(m)$ and $K'(m)=K(\sqrt{1-m^{2}})$. We observe that $\bar{u}$ is the average value of u(x, t) over an oscillation. The above formula coincides with the genus-one case of the more general Its-Matveev and Dubrovin-Novikov formula [Its75, Dub74] for finite-gap solutions of KdV.

With the above potential (1.5), the Schrödinger Eq. (1.2) coincides with the Lamé equation and the stability zones (or Bloch spectrum) of the potential are $[\beta _3,\beta _2]\cup [\beta _1,+\infty )$.

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.21.3) 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 small-dispersion 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 the evolution of random configurations of solitons. In [DP14, 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 [ET20, 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 not random, can be naturally interpreted as a soliton gas. Cutting to the chase, the authors derive a Riemann–Hilbert problem which seeks a vector $\Xi = \begin{bmatrix} \Xi _{1}&\Xi _{2} \end{bmatrix}^{T}$ satisfying a normalization condition at $\infty$, and the jump relations

$$\begin{aligned} \Xi _{+}(i \lambda ) = J(\lambda ) \Xi _{-}( i \lambda ) \, , \qquad \Xi _{+}(-i \lambda ) = J^{T}(\lambda ) \Xi _{-}( -i \lambda ) \, , \qquad \lambda \in (\eta _1,\eta _2)\, \end{aligned}$$

(1.7)

where the jump matrix $J(\lambda )$ is given by

$$\begin{aligned} J(\lambda ) = \frac{1}{ 1 +r_1(\lambda )r_2(\lambda )}\begin{bmatrix} \displaystyle 1 -r_1(\lambda )r_2(\lambda ) &{} \displaystyle 2 ir_1(\lambda ) e^{ - 2 \lambda x} \\ \displaystyle 2 ir_2(\lambda ) e^{2 \lambda x} &{} \displaystyle 1 -r_1(\lambda ) r_2(\lambda ) \end{bmatrix} \ . \end{aligned}$$

(1.8)

The parameters $\eta _{1}$ and $\eta _{2}$ are taken to be real with $0< \eta _{1} < \eta _{2}$, and the intervals $(i\eta _1,i\eta _2)$ and $(-i\eta _2,-i\eta _1)$ are oriented downwards.

The reflection coefficients $r_1(\lambda )=r_1(\lambda ; t)$ and $r_2(\lambda )=r_2(\lambda ; t)$ evolve in time according to

$$\begin{aligned} r_1(\lambda ; t)=r_1(\lambda ; 0) e^{ ( 8 \lambda ^{3} - 12 \lambda ) t} \, , \qquad r_2(\lambda ; t) =r_2(\lambda ; 0)e^{-( 8 \lambda ^{3} - 12 \lambda ) t} \ . \end{aligned}$$

(1.9)

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} \equiv 0$, the potential is exponentially decaying as $x \rightarrow +\infty$. 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, 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, at $t=0$ the potential approaches zero exponentially fast as $x \rightarrow +\infty$, while for $x\rightarrow -\infty$ the potential approaches the cnoidal wave solution of KdV very slowly—the error decays with a rate of $O(\frac{1}{x})$. Because of these different behaviors, these potentials represent a new large class of potentials which have not been previously considered in the literature. This model is substantially different from the model of infinite solitons considered in [Boy84, Zai83] where an infinite number of equally spaced and identical solitons can be identified with the cnoidal wave solution of KdV.

1.2 Statement of the results

In Sect. 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\eta _2,-i\eta _1)\cup (i\eta _1,i\eta _2)$ for some $\eta _2>\eta _1>0$ and show that for this sequence, as $N \rightarrow + \infty$, the solution of the Riemann–Hilbert problem converges to the solution of the Riemann–Hilbert problem studied in [DZZ16], for the case $r_2(\lambda ) \equiv 0$.

Remark. Since the Riemann–Hilbert problem emerges in a limit, the existence and uniqueness of a solution is not a-priori known. For completeness, we provide a proof of existence which is valid for all x and t in the Appendix.

In Sect. 3 (Theorem 3.6) we establish that the potential u(x, 0) determined by this Riemann–Hilbert problem coincides with the periodic travelling wave as $x \rightarrow - \infty$:

$$\begin{aligned} u(x,0)=\eta _2^2-\eta _1^2-2\eta _2^2{\text {dn}}^2\left( \eta _2(x+\phi ) + K(m)\left| \, m\right. \right)+ \mathcal {O}\left(x^{-1}\right) \,. \end{aligned}$$

(1.10)

The function ${\text {dn}}\left( z\left| \, m\right. \right)$ is the Jacobi elliptic function of modulus $m=\eta _1/\eta _2$. It is periodic with period 2K(m), and satisfies ${\text {dn}}\left( 0 \left| \, m\right. \right)=1$ and ${\text {dn}}\left( K(m) \left| \, m\right. \right)=\sqrt{1-m^2}$. The expression (1.10) for the elliptic solution of KdV coincides with the travelling wave solution (1.5) in the introduction by identifying $\beta _1=0$, $\beta _2=-\eta _1^2$ and $\beta _3=-\eta _2^2$.

The function (1.10) is periodic in x with period $2K(m)/\eta _2$.The minimum amplitude of the oscillations is $-\eta _2^2-\eta _1^2$ and the maximum amplitude is $\eta _1^2-\eta _2^2$ so that the amplitude of the oscillations is $2\eta _1^2$. The average value of u(x) over an oscillation can be obtained from (1.6).

The phase $\phi$ in formula (1.10) depends on the coefficient $r_1(\lambda )$ that characterizes the continuum limit of the norming constants of the soliton gas and it is equal to

$$\begin{aligned} \phi =\int _{\eta _1}^{\eta _2}\dfrac{\log 2r_1(i\zeta )}{\sqrt{(\zeta ^2-\eta _1^2)(\zeta ^2-\eta _2^2)}}\dfrac{\mathrm{d}\zeta }{\pi i}\in \mathbb {R}\, . \end{aligned}$$

(1.11)

Remark 1.1

The potential u(x, 0) is a step-like finite gap potential. The slow decay rate as $x \rightarrow -\infty$ implies that such potential does not fall in the class considered in [BdMET08]. When $\eta _1=0$ the potential $u(x,0)=-\eta _2^2+ \mathcal {O}\left(x^{-1}\right)$ as $x \rightarrow - \infty$. Such a potential is a step-like potential with zero reflection coefficient on the real axis. It is not included in the class of potentials studied in [EGKT13, CK85] because of the low decaying condition at $x\rightarrow -\infty$. Potentials with a low decay rate have appeared when studying rogue waves of infinite order of the focusing nonlinear Schrödinger Eq. [BM19], see also [BLM20].

Finally in Sects. 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, 0). The asymptotic behaviour depends on the quantity $\xi = x / 4 t$. 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; and (3) a region where the solution is approximated by a periodic travelling wave with modulated coefficients (see Fig. 1). More precisely:

• (1) for fixed $\xi > \eta _{2}^{2}$, there is a positive constant $C = C(\xi )$ so that

  $$\begin{aligned} u(x,t)= \mathcal {O}\left(e^{- C t}\right) \, . \end{aligned}$$

• (2) For $\xi <\xi _\mathrm{crit}$ we have

  $$\begin{aligned} u(x,t)=\eta _2^2-\eta _1^2-2\eta _2^2{\text {dn}}^2\left( \eta _2(x-2(\eta _1^2+\eta _2^2)t+\phi ) + K(m)\left| \, m\right. \right)+ \mathcal {O}\left(t^{-1}\right) \, , \end{aligned}$$

  (1.12)

  with $m=\eta _1/\eta _2$ and $\phi$ as in (1.11). The critical value $\xi _\mathrm{crit}$ is obtained from the equation

  $$\begin{aligned} \xi _\mathrm{crit}=\dfrac{\eta _2^2}{2}W(m)\, ,\quad W(m)=1+m^2+2\dfrac{m^2(1-m^2)}{1-m^2-\frac{E(m)}{K(m)}}\, ,\quad m=\frac{\eta _1}{\eta _2}\, . \end{aligned}$$

  (1.13)

• (3) For $\xi _\mathrm{crit}<\xi <\eta _2^2$ we have that

  $$\begin{aligned} u(x,t)=\eta _2^2-\alpha ^2-2\eta _2^2{\text {dn}}^2\left( \eta _2(x-2(\alpha ^2+\eta _2^2)t+\widetilde{\phi }) + K(m_\alpha ) \left| \, m_\alpha \right. \right)+ \mathcal {O}\left(t^{-1}\right) \, , \end{aligned}$$

  (1.14)

  where ${\text {dn}}\left( z\left| \, m_{\alpha }\right. \right)$ is the Jacobi elliptic function of modulus $m_\alpha =\alpha /\eta _2$,

  $$\begin{aligned} \widetilde{\phi }=\int _{\alpha }^{\eta _2}\dfrac{\log 2r_1(i\zeta )}{\sqrt{(\zeta ^2-\alpha ^2)(\zeta ^2-\eta _2^2)}}\dfrac{\mathrm{d}\zeta }{\pi i}\in \mathbb {R}\, \end{aligned}$$

  and the coefficient $\alpha =\alpha (\xi )$ is determined from the Whitham modulation equation [Whi74]

  $$\begin{aligned} \xi =\dfrac{x}{4t}=\dfrac{\eta _2^2}{2}W(m_\alpha ) \, , \end{aligned}$$

  (1.15)

  where W(m) has been defined in (1.13).

The Eq. (1.15) 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)=-\eta _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].

Fig. 1

Soliton gas behaviour at $t=10$ with endpoints $\eta _1=0.5$ and $\eta _2=1.5$ and reflection coefficient $r_1(\lambda ) \equiv 1$.

2 Soliton Gas as Limit of N Solitons as $N\rightarrow +\infty$

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

• (i) $M(\lambda )$ is meromorphic in $\mathbb {C}$, with simple poles at $\{\lambda _{j} \}_{j=1}^{N}$ in $i\mathbb {R}_{+}$, and at the corresponding conjugate points $\{\overline{\lambda _{j}} \}_{j=1}^{N}$ in $i\mathbb {R}_{-}$;

• (ii) M satisfies the residue conditions

  $$\begin{aligned} \mathop {\mathrm{res}}_{\lambda =\lambda _{j}} M (\lambda )= \lim _{\lambda \rightarrow \lambda _{j}} M(\lambda ) \begin{bmatrix} 0 &{} 0 \\ \displaystyle \frac{c_j e^{2 i \lambda _{j} x }}{N} &{} 0 \end{bmatrix}\, , \quad \mathop {\mathrm{res}}_{\lambda =\overline{\lambda _{j}}} M(\lambda ) = \lim _{\lambda \rightarrow \overline{\lambda _{j}}} M(\lambda ) \begin{bmatrix} 0 &{} \displaystyle \frac{- c_j e^{-2 i \overline{\lambda _{j}} x }}{N} \\ 0 &{} 0\end{bmatrix} \, , \end{aligned}$$

  (2.1)

  where $c_j\in i\mathbb {R}_+$;

• (iii) $\displaystyle M(\lambda )= \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right)$ as $\lambda \rightarrow \infty$,

• (iv) M satisfies the symmetry

  $$\begin{aligned} M(-\lambda )=M(\lambda )\begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}\,. \end{aligned}$$

The solution of the above Riemann–Hilbert problem is determined from the relation

$$\begin{aligned} M(\lambda ) = \left( 1 + \sum _{j=1}^{N} \frac{e^{i\lambda _j x}\alpha _{j}}{\lambda - \lambda _{j}} , \ 1 - \sum _{j=1}^{N} \frac{e^{i\lambda _j x}\alpha _{j}}{\lambda + \lambda _{j}} \right) \ , \end{aligned}$$

(2.2)

where the constants $\alpha _j$ are uniquely determined by the residue conditions (2.1). The N-soliton potential u(x) is determined from M via

$$\begin{aligned} u(x) = 2 \frac{\mathrm{d}}{\mathrm{d}x} \left( \lim _{\lambda \rightarrow \infty } \frac{ \lambda }{i} \left( M_{1}(\lambda ) - 1 \right) \right) \ , \end{aligned}$$

(2.3)

where $M_1(\lambda )$ is the first entry of the vector $M(\lambda )$. In particular, for a one-soliton potential, namely $N=1$, one recovers the expression (1.4) where the shift $x_0$ is given by

$$\begin{aligned} x_0=\dfrac{1}{4\eta _1}\log \frac{c_1}{2i\eta _1}\in \mathbb {R}. \end{aligned}$$

We are interested in the limit as $N \rightarrow + \infty$, under the additional assumptions:

• (i) The poles $\left\{ \lambda _{j}^{(N)} \right\} _{j=1}^{N}$ are sampled from a density function $\varrho (\lambda )$ so that $\int _{\eta _{1}}^{-i \lambda _{j}} \varrho ( \eta ) d \eta = j / N$, for $j = 1, \ldots , N$.

• (ii) The coefficients $\{c_j\}_{j=1}^N$ are purely imaginary (in fact $c_{j} \in i \mathbb {R}_{+}$) and are assumed to be discretizations of a given function:

  $$\begin{aligned} c_j = \frac{i (\eta _2-\eta _1) r_1(\lambda _j)}{\pi } \qquad j=1,\ldots , N \ . \end{aligned}$$

  (2.4)

  where $r_1(\lambda )$ is an analytic function for $\lambda$ near the intervals $(i \eta _{1}, i \eta _{2})$ and $(-i \eta _{2}, -i \eta _{1})$, with the symmetry $r_1(\overline{\lambda }) = r_1(\lambda )$, and is further assumed to be a real valued positive and non-vanishing function of $\lambda$ for $\lambda \in [i \eta _{1}, i \eta _{2}]$.

In the regime $x \rightarrow + \infty$, it is easy to notice that all residue conditions (2.1) contain only exponentially small terms and therefore, by a small norm argument, the potential is exponentially small.

On the other hand, for $x \rightarrow - \infty$ 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

$$\begin{aligned} A(\lambda )=M(\lambda ) \prod _{j=1}^{N} \left( \frac{\lambda - \lambda _j }{\lambda -\overline{\lambda _j}} \right) ^{\sigma _{3}} \ . \end{aligned}$$

(2.5)

Now the residue conditions are

$$\begin{aligned}&\mathop {\mathrm{res}}_{\lambda =\lambda _{j}} A(\lambda ) = \lim _{\lambda \rightarrow \lambda _{j}} A(\lambda ) \begin{bmatrix} 0 &{}\displaystyle \frac{N}{c_j} e^{-2 i \lambda _{j} x }(\lambda _{j}-\overline{\lambda _{j}})^{2}\prod _{k\ne j} \left( \frac{\lambda _{j}-\overline{\lambda _{k}}}{\lambda _{j}-\lambda _{k}}\right) ^{2} \\ 0 &{} 0 \end{bmatrix} \end{aligned}$$

(2.6)

$$\begin{aligned}&\mathop {\mathrm{res}}_{\lambda =\overline{\lambda _{j}}} A(\lambda ) = \lim _{\lambda \rightarrow \overline{\lambda _{j}}} A(\lambda ) \begin{bmatrix} 0 &{} 0 \\ \displaystyle \frac{-N}{c_j} e^{2 i \overline{\lambda _{j}} x } (\overline{\lambda _{j}}-\lambda _{j})^{2} \prod _{k\ne j} \left( \frac{\overline{\lambda _{j}}-\lambda _{k}}{\overline{\lambda _{j}}-\overline{\lambda _{k}}}\right) ^{2} &{} 0 \end{bmatrix} \end{aligned}$$

(2.7)

while the potential u(x) is still extracted from A via the same calculation:

$$\begin{aligned} u(x) = 2 \frac{\mathrm{d}}{\mathrm{d}x} \left( \lim _{\lambda \rightarrow \infty } \frac{\lambda }{i} \left( A_{1}(\lambda ) - 1 \right) \right) \, . \end{aligned}$$

The quantity $e^{ -2 i \lambda _{j} x}$ now decays exponentially as $x \rightarrow - \infty$, and this implies (again by a standard small-norm argument) exponential decay of the potential u(x) for $x \rightarrow - \infty$. On the other hand, the product term is exponentially large in N. One may show that there is $C>0$ so that

$$\begin{aligned} \left| \frac{N}{c_{j}} ( \lambda _{j} - \overline{\lambda _{j}}) \prod _{k\ne j } \left( \frac{\lambda _{j}-\overline{\lambda _{k}}}{\lambda _{j}-\lambda _{k}} \right) ^{2} \right| < D e^{CN} \qquad \text {for all } j=1,\ldots , N \ . \end{aligned}$$

(2.8)

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 \ll - C N$. 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 \rightarrow -\infty$. To be more precise, the above computations can be used to show the following lemma.

Lemma 2.1

For any $0< \tilde{k} < \eta _{1}/2$, there exists a constant $\tilde{C}$ so that if $\displaystyle x < - \tilde{C} N$, then

$$\begin{aligned} |u(x)| < e^{-\tilde{k} |x|} \ . \end{aligned}$$

(2.9)

In other words, for $x < - \tilde{C} N$, the potential u(x) is exponentially decreasing.

The proof of this lemma is straightforward: under the hypotheses of the lemma, all residue conditions are exponentially small. One may replace the residue conditions with jumps across small circles encircling the poles, and the jumps are all of the form $I + \mathcal {O}\left( e^{-\tilde{k} |x|} \right)$, so small norm existence theory applies.

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 \rightarrow +\infty$. First, we remove the poles by defining

$$\begin{aligned} Z(\lambda ) = M(\lambda ) \begin{bmatrix} 1 &{} 0 \\ \displaystyle - \frac{1}{N}\sum _{j=1}^{N}\frac{c_j e^{2 i \lambda x }}{\lambda - \lambda _{j}} &{} 1\end{bmatrix} \end{aligned}$$

(2.10)

within a closed curve $\gamma _+$ encircling the poles counterclockwise in the upper half plane $\mathbb {C}_{+}$, and

$$\begin{aligned} Z(\lambda ) = M(\lambda ) \begin{bmatrix}{1} &{} \displaystyle { \frac{1}{N}\sum _{k=1}^{N}\frac{c_j e^{-2 i \lambda x }}{\lambda - \overline{\lambda _{j}}}}\\ {0}&{}{1} \end{bmatrix} \end{aligned}$$

(2.11)

within a closed curve $\gamma _-$ surrounding the poles clockwise in the lower half plane $\mathbb {C}_{-}$. Outside these two sets, we take $Z(\lambda )=M(\lambda )$.

Then the jumps are

$$\begin{aligned} Z_+(\lambda ) = Z_-(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} 1 &{} 0 \\ \displaystyle - \frac{1}{N}\sum _{j=1}^{N}\frac{c_j e^{2 i \lambda x }}{\lambda - \lambda _{j}} &{} 1\end{bmatrix} &{} \lambda \in \gamma _+ \\ \begin{bmatrix}{1} &{} \displaystyle { -\frac{1}{N}\sum _{k=1}^{N}\frac{c_j e^{-2 i \lambda x }}{\lambda - \overline{\lambda _{j}}}}\\ {0}&{}{1} \end{bmatrix}&\lambda \in \gamma _- \ \end{array}\right. } \end{aligned}$$

(2.12)

where, for $\lambda \in \gamma _{+}$ or $\gamma _{-}$, the boundary values $Z_{+}(\lambda )$ are taken from the left side of the contour as one traverses it according to its orientation, and the boundary values $Z_{-}(\lambda )$ are taken from the right. The quantity $Z(\lambda )$ is normalized so that $Z(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\lambda ^{-1}\right)$ as $\lambda \rightarrow \infty$.

We assume now that in the limit as the number of poles goes to infinity, the poles are distributed according to some distribution $\varrho (\lambda )$ with density compactly supported in $(i\eta _1,i\eta _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\eta _1,i\eta _2)$ with distance between two poles equal to $| \Delta \lambda | = \frac{\eta _2-\eta _1}{N}$ and with (atomic) density:

$$\begin{aligned} \varrho _{N}(\lambda ) = \frac{1}{Z_N} \sum _{j=1}^N c_j \delta _{\lambda _j}(\lambda ) \qquad \lambda \in (i\eta _1,i\eta _2) \ , \end{aligned}$$

(2.13)

for some normalization constant $Z_N$.

Remark 2.2

In the case where the poles are distributed according to a more general measure $\varrho (\lambda )$, 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(\lambda )$.

As the number of poles increases within the support of the measure, the following result holds.

Proposition 2.3

For any open set $K_{+}$ containing the interval $[i \eta _{1}, i \eta _{2}]$, and any open set $K_{-}$ containing the interval $[-i \eta _{2}, - i \eta _{1}]$, the following limit holds uniformly for all $\lambda \in \mathbb {C }\backslash K_{+}$:

$$\begin{aligned} \lim _{N\rightarrow +\infty } \frac{1}{N}\sum _{j=1}^{N}\frac{c_j }{\lambda - \lambda _{j}} =\int _{i\eta _1}^{i\eta _2} \frac{2i r_1(\zeta )}{\lambda - \zeta }\frac{\mathrm{d}\zeta }{2\pi i} \,, \end{aligned}$$

(2.14)

and the following limit holds uniformly for all $\lambda \in \mathbb {C} \backslash K_{-}$:

$$\begin{aligned} \lim _{N\rightarrow +\infty } \frac{1}{N}\sum _{j=1}^{N}\frac{c_j }{\lambda - \overline{\lambda _{j}}} = \int _{-i\eta _2}^{-i\eta _1} \frac{2i r_1(\zeta )}{\lambda - \zeta }\frac{\mathrm{d}\zeta }{2\pi i} \ , \end{aligned}$$

(2.15)

where $r_1(\lambda )$ is an analytic function for $\lambda$ near the intervals $(i \eta _{1}, i \eta _{2})$ and $(-i \eta _{2}, -i \eta _{1})$, and it is assumed to be a positive real-valued and non-vanishing function of $\lambda$ for $\lambda \in [i \eta _{1}, i \eta _{2}]$.

Proof

Using (2.4), the expressions in the jumps can be rewritten as

$$\begin{aligned} \frac{1}{N}\sum _{j=1}^{N}\frac{c_j }{\lambda - \lambda _{j}} = \frac{1}{N}\sum _{j=1}^{N}\frac{1}{\lambda - \lambda _{j}} \frac{(\eta _2-\eta _1)i r_1(\lambda _j)}{\pi } = \frac{1}{2\pi i}\sum _{j=1}^{N}\frac{ 2i r_1(\lambda _j) }{\lambda - \lambda _{j}} \Delta \lambda \ .\nonumber \\ \end{aligned}$$

(2.16)

The convergence follows from the convergence of the Riemann sum to the Riemann–Stieltjes integral for $x \in K$ any compact subset of $\mathbb {R}$. Positivity of $r_1(\lambda )$ follows from the fact that $c_j\in i\mathbb {R}_+$.$\square$

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)

$$\begin{aligned}&Z_{+} (\lambda )= Z_{-}(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix}{1}&{} {0}\\ \displaystyle {e^{2i\lambda x}\int _{i\eta _1}^{ i\eta _2} \frac{ 2i r_1(\zeta )}{\zeta - \lambda }\frac{\mathrm{d}\zeta }{2\pi i}} &{} {1} \end{bmatrix} \qquad \lambda \in \gamma _+ \\ \begin{bmatrix} {1}&{}\displaystyle {e^{-2i\lambda x} \int _{-i\eta _2}^{-i\eta _1} \frac{ 2i r_1(\zeta )}{ \zeta - \lambda }\frac{\mathrm{d}\zeta }{2\pi i} } \\ {0}&{} {1} \end{bmatrix} \qquad \lambda \in \gamma _- \end{array}\right. } \end{aligned}$$

(2.17)

$$\begin{aligned}&Z(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned}$$

(2.18)

At this point it is important to point out to the reader that the contour $(i \eta _{1}, i \eta _{2})$ and $(-i \eta _{2}, -i \eta _{1})$ are both oriented upwards.

Next, we define

$$\begin{aligned} X(\lambda ) = Z(\lambda ) \begin{bmatrix} {1}&{} {0}\\ -\displaystyle { e^{ 2 i \lambda x} \int _{i\eta _1}^{ i\eta _2} \frac{2i r_1(\zeta )}{\zeta -\lambda } \frac{\mathrm{d}\zeta }{2\pi i} }&{} {1}\end{bmatrix} \end{aligned}$$

(2.19)

within the loop $\gamma _+$, and

$$\begin{aligned} X(\lambda )= Z(\lambda ) \begin{bmatrix} {1}&{}\displaystyle { e^{- 2 i \lambda x} \int _{-i\eta _2}^{-i\eta _1} \frac{2i r_1(\zeta )}{ \zeta -\lambda } \frac{\mathrm{d}\zeta }{2\pi i} }\\ {0}&{} {1} \end{bmatrix} \end{aligned}$$

(2.20)

within the loop $\gamma _-$. Outside these two curves, we define $X(\lambda ) = Z(\lambda )$.

The jumps across the curves are no longer present, but there are jumps across $(i \eta _{1}, i \eta _{2})$ and $(-i \eta _{2}, -i \eta _{1})$ because the integrals have jumps across those intervals. Using the Sokhotski-Plemelj formula, we arrive at a Riemann–Hilbert problem for X:

$$\begin{aligned}&X_{+}(\lambda ) = X_{-}(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} {1}&{} {0} \\ \displaystyle {-2i r_1(\lambda ) e^{ 2 i \lambda x} } &{} {1} \end{bmatrix} &{} \lambda \in (i\eta _1,i\eta _2) \\ \begin{bmatrix} {1}&{} \displaystyle {2i r_1(\lambda ) e^{ -2 i \lambda x} } \\ {0} &{} {1} \end{bmatrix}&\lambda \in (-i\eta _2,-i\eta _1) \end{array}\right. } \end{aligned}$$

(2.21)

$$\begin{aligned}&X(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ , \end{aligned}$$

(2.22)

$$\begin{aligned}&X(-\lambda )=X(\lambda )\begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}\,. \end{aligned}$$

(2.23)

This Riemann–Hilbert problem is equivalent to the one described in [DZZ16] with $r_2(\lambda ) = 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(\overline{\lambda }) = r_1(\lambda )$ for $\lambda \in (-i \eta _{1}, - i \eta _{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.

Since this Riemann–Hilbert problem has been derived through a limiting process, it is not at all clear that it actually possesses a solution. Although this will eventually follow for large x from the asymptotic analysis presented herein, we present a self-contained proof of existence and uniqueness in Appendix A of this paper. In fact in the appendix we establish the existence of a matrix-valued solution $\mathbf{Y}$. Equipped with that, it is straightforward to prove the following lemma.

Lemma 2.4

Let B be an arbitrary positive number. For all x such that $|x|<B$, the quantity Z defined in (2.10)–(2.11) satisfying the jump relation (2.12) converges as $N \rightarrow \infty$ to the solution of the Riemann–Hilbert problem (2.17)–(2.18), and the N-soliton potential u(x) converges to the potential determined by the solution to the soliton gas Riemann-Hilbert problem (2.21)–(2.23).

Note: A similar limiting procedure introduced in this section has already appeared in the literature when studying the focusing nonlinear Schrödinger equation [BLM20, BB19]. In those papers the limiting procedure $N\rightarrow \infty$ refers to the order of a soliton or a breather of the nonlinear Schrödinger equation. From a different point of view, our limiting procedure can be understood as replacing the reflection coefficient with its semiclassical limit, see for example [TVZ04, LL83a].

In what follows we have already taken the $N \rightarrow \infty$ limit, and we are considering the behavior for x (and later x and t) large for this soliton gas Riemann–Hilbert problem.

3 Behaviour of the Potential u(x, 0) as $x\rightarrow -\infty$

We consider a soliton gas Riemann–Hilbert problem as in (2.21)–(2.22) with $0< \eta _{1} < \eta _{2}$, and reflection coefficient $r_1(\lambda )$ defined on $( i \eta _{1}, i \eta _{2})$ such that it has an analytic extension to a neighbourhood of this interval. Furthermore, we assume that $r_1(-\lambda )=r_1(\lambda )$ on the imaginary axis. We set $\Sigma _1 = (\eta _1,\eta _2)$ and $\Sigma _2 = (-\eta _2,-\eta _1)$. The vector valued function X that will determine the KdV potential u(x) is the solution to the following Riemann–Hilbert problem:

$$\begin{aligned}&X(\lambda ) \text { is analytic for } \lambda \in \mathbb {C}\backslash \left\rbrace i\Sigma _1 \cup i\Sigma _2 \right\lbrace \nonumber \\&X_+(i\lambda ) = X_-(i\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} 1 &{} 0 \\ -2i r_1(i\lambda )e^{-2\lambda x} &{} 1 \end{bmatrix} &{}\quad \lambda \in \Sigma _1\\ \displaystyle \begin{bmatrix} 1 &{}2i r_1(i\lambda ) e^{2\lambda x} \\ 0 &{} 1 \end{bmatrix}&\quad \lambda \in \Sigma _2 \end{array}\right. } \nonumber \\&X (\lambda ) = \begin{bmatrix}1&1\end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned}$$

(3.1)

As explained in [DZZ16], we can recover the potential u(x) of the Schrödinger operator via the formula

$$\begin{aligned} u(x) = 2\frac{\mathrm{d}}{\mathrm{d}x} \left[ \lim _{\lambda \rightarrow \infty } \frac{\lambda }{i} (X_1(\lambda ;x) - 1) \right] \ , \end{aligned}$$

(3.2)

where $X_1(\lambda ;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

$$\begin{aligned} Y (\lambda )= X(i\lambda ) \, ,\quad r(\lambda )=2r_1(i\lambda )\, , \end{aligned}$$

(3.3)

the Riemann–Hilbert problem for Y reads as follows:

$$\begin{aligned}&Y_+ (\lambda ) = Y_-(\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} 1 &{}0 \\ -i r(\lambda ) e^{-2\lambda x} &{} 1 \end{bmatrix} &{}\quad \lambda \in \Sigma _1\\ \displaystyle \begin{bmatrix} 1 &{} i r(\lambda )e^{2\lambda x} \\ 0 &{} 1 \end{bmatrix}&\quad \lambda \in \Sigma _2 \end{array}\right. } \end{aligned}$$

(3.4)

$$\begin{aligned}&Y(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \end{aligned}$$

(3.5)

$$\begin{aligned}&Y(-\lambda ) = Y(\lambda )\begin{bmatrix}0&{}1\\ 1&{}0 \end{bmatrix} \,. \end{aligned}$$

(3.6)

The contours $\Sigma _1$ and $\Sigma _2$ are shown in Fig. 2. We can recover u(x) from

$$\begin{aligned} u(x) = 2\frac{\mathrm{d}}{\mathrm{d}x} \left[ \lim _{\lambda \rightarrow \infty } \lambda (Y_1(\lambda ;x) - 1) \right] \, . \end{aligned}$$

(3.7)

Fig. 2

Riemann–Hilbert problem for Y

3.1 Large x asymptotic

Introduce the following new vector function

$$\begin{aligned} T(\lambda )=Y(\lambda )e^{x g(\lambda )\sigma _3}f(\lambda )^{\sigma _3} \end{aligned}$$

where $g(\lambda )$ and $f(\lambda )$ are scalar functions to be determined below. We require that

• $g(\lambda )$ is analytic in $\mathbb {C}\backslash [-\eta _2, \eta _2]$ and

  $$\begin{aligned}&g_+(\lambda ) + g_-(\lambda ) = 2\lambda&\lambda \in \Sigma _1\cup \Sigma _2 \end{aligned}$$

  (3.8)
  $$\begin{aligned}&g_+ (\lambda ) - g_-(\lambda ) = \Omega&\lambda \in [-\eta _1,\eta _1] \end{aligned}$$

  (3.9)
  $$\begin{aligned}&g(\lambda ) =\mathcal {O}\left(\frac{1}{\lambda }\right)&\lambda \rightarrow \infty \, , \end{aligned}$$

  (3.10)

  where $\Omega$ is a constant independent of x, still to be determined, and

• $f(\lambda )$ is analytic in $\mathbb {C}\backslash [-\eta _2, \eta _2]$ and

  $$\begin{aligned} f(\lambda )=1+\mathcal {O}\left(\frac{1}{\lambda }\right) as \quad \lambda \rightarrow \infty . \end{aligned}$$

  (3.11)

In order to solve the scalar Riemann–Hilbert problem (3.8)–(3.10) for g we observe that

$$\begin{aligned}&g'_+(\lambda ) +g'_-(\lambda ) =2&\lambda \in \Sigma _1\cup \Sigma _2 \end{aligned}$$

(3.12)

$$\begin{aligned}&g'_+(\lambda ) -g'_-(\lambda ) =0&\lambda \in [-\eta _1,\eta _1] \end{aligned}$$

(3.13)

$$\begin{aligned}&g'(\lambda ) =\mathcal {O}\left(\frac{1}{\lambda ^2}\right)&\lambda \rightarrow \infty \ . \end{aligned}$$

(3.14)

From the above, we can write $g'(\lambda )$ as

$$\begin{aligned} g'(\lambda ) = 1 - \frac{\lambda ^2+\kappa }{R(\lambda )}\, , \end{aligned}$$

(3.15)

where

$$\begin{aligned} R(\lambda ) = \sqrt{(\lambda ^2-\eta _1^2)(\lambda ^2-\eta _2^2)}\, , \end{aligned}$$

(3.16)

is real and positive on $(\eta _2,+\infty )$ with branch cuts on the contours $\Sigma _1$ and $\Sigma _2$ and $\kappa$ is a constant to be determined. By integration we obtain

$$\begin{aligned} g(\lambda )=\lambda -\int _{\eta _2}^{\lambda }\frac{\zeta ^2+\kappa }{R(\zeta )}\mathrm{d}\zeta \, . \end{aligned}$$

(3.17)

Condition (3.8) implies that

$$\int _{-\eta _1}^{\eta _1}\frac{\zeta ^2+\kappa }{R(\zeta )}\mathrm{d}\zeta =0$$

and condition (3.9) implies that

$$\begin{aligned} \Omega =2\int _{\eta _1}^{\eta _2}\frac{\zeta ^2+\kappa }{R_+(\zeta )}\mathrm{d}\zeta \, . \end{aligned}$$

This gives

$$\begin{aligned} \Omega = \frac{2\pi i}{\int _{-\eta _1}^{\eta _1}\frac{\mathrm{d}\zeta }{R(\zeta )} } =-\dfrac{i\pi \eta _2}{K(m)}\in i \mathbb {R}_- \ ,\quad m=\frac{\eta _1}{\eta _2}\, , \end{aligned}$$

(3.18)

where $K(m)=\int _0^{\frac{\pi }{2}}\frac{\mathrm{d}\vartheta }{\sqrt{1-m^2\sin \vartheta }}$ is the complete elliptic integral of the first kind with modulus $m=\eta _1/\eta _2$ and

$$\begin{aligned} \kappa = - \int _{-\eta _1}^{\eta _1} \frac{ \zeta ^2 \Omega }{R(\zeta )} \frac{\mathrm{d}\zeta }{2\pi i}=\eta _2^2\left( \dfrac{E(m)}{K(m)}-1\right) \quad \in \mathbb {R}_- \ , \end{aligned}$$

(3.19)

where $E(m)=\int _0^{\frac{\pi }{2}}\sqrt{1-m^2\sin \vartheta }\, \mathrm{d}\vartheta$ is the complete elliptic integral of the second kind.

The Riemann–Hilbert problem for $T(\lambda )$ is

$$\begin{aligned}&T_+(\lambda )=T_-(\lambda )V_T(\lambda ) \end{aligned}$$

(3.20)

$$\begin{aligned}&V_T(\lambda )={\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} e^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )} &{} 0 \ i r(\lambda )f_+(\lambda )f_-(\lambda ) &{} e^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix} &{}\quad \lambda \in \Sigma _1\\ \displaystyle \begin{bmatrix} e^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )}&{}\dfrac{ i r(\lambda )}{f_+(\lambda )f_-(\lambda )} \\ 0 &{} e^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix} &{} \quad \lambda \in \Sigma _2 \\ \displaystyle \begin{bmatrix} e^{x\Omega }\dfrac{f_+(\lambda )}{f_-(\lambda )} &{} 0 \\ 0&{} e^{-x \Omega }\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix} &{}\quad \lambda \in [-\eta _1,\eta _1]\\ \end{array}\right. } \end{aligned}$$

(3.21)

$$\begin{aligned}&T(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned}$$

(3.22)

In order to solve the Riemann–Hilbert problem for $T(\lambda )$ we wish to obtain a constant jump matrix $J_T$. For this purpose we make the following ansatz on the function f

$$\begin{aligned}&f_+(\lambda ) f_-(\lambda ) = \frac{1}{r(\lambda )}&\lambda \in \Sigma _1 \end{aligned}$$

(3.23)

$$\begin{aligned}&f_+(\lambda ) f_-(\lambda ) = r(\lambda )&\lambda \in \Sigma _2\end{aligned}$$

(3.24)

$$\begin{aligned}&\frac{f_+(\lambda )}{ f_-(\lambda ) }=e^{\Delta }&\lambda \in [-\eta _1,\eta _1] \end{aligned}$$

(3.25)

$$\begin{aligned}&f(\lambda ) = 1+\mathcal {O}\left(\frac{1}{\lambda }\right)&\lambda \rightarrow \infty \, . \end{aligned}$$

(3.26)

It is easy to check that the function $f(\lambda )$ is given by

$$\begin{aligned} f(\lambda )={\text {exp}}\left\rbrace \frac{R(\lambda )}{2\pi i}\left[ \int _{ \Sigma _1} \frac{\log \frac{1}{r(\zeta )} }{R_+(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta + \int _{\Sigma _2} \frac{\log r(\zeta ) }{R_+(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta +\int ^{\eta _1}_{ -\eta _1} \frac{\Delta }{R(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta \right] \right\} . \end{aligned}$$

(3.27)

The inclusion of the constant jump (3.25) allows us to satisfy (3.26) by taking

$$\begin{aligned} \Delta&=\left[ \int _{\Sigma _1} \frac{\log r(\zeta )}{R_+(\zeta )} \mathrm{d}\zeta - \int _{\Sigma _2} \frac{\log r(\zeta )}{R_+(\zeta )} \mathrm{d}\zeta \right] \left[ \int ^{\eta _1}_{ -\eta _1} \frac{ \mathrm{d}\zeta }{R(\zeta )}\right] ^{-1} = -\dfrac{\eta _2}{K(m)}\int ^{ \eta _2}_{ \eta _1} \frac{\log r(\zeta )}{R_+(\zeta )} \mathrm{d}\zeta \, , \end{aligned}$$

(3.28)

where in the last equality in (3.28) we use the fact that $r(-\lambda )=r(\lambda )$. We remind the reader that we are assuming the function r to be real, positive, and non-vanishing on $\Sigma _{1}$ and $\Sigma _{2}$. The positivity of $r(\lambda )$ guarantees that $\Delta$ is pure imaginary.

3.2 Opening lenses

We start by defining the analytic continuation $\hat{r}(\lambda )$ of the function $r(\lambda )$ off the interval $(-\eta _2,-\eta _1)\cup (\eta _1,\eta _2)$ with the requirement that

$$\begin{aligned} \hat{r}_{\pm }(\lambda )=\pm r(\lambda )\, , \quad \lambda \in (-\eta _2,-\eta _1)\cup (\eta _1,\eta _2)\, . \end{aligned}$$

(3.29)

We can factor the jump matrix $J_T$ on $\Sigma _{1}$ as follows

$$\begin{aligned}&\begin{bmatrix} e^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )} &{} 0 \\ -i &{} e^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix} \\&\quad \quad =\begin{bmatrix}1&{}-\dfrac{ie^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}}{\hat{r}_-(\lambda )f_-^2 (\lambda )}\\ 0 &{}1 \end{bmatrix} \begin{bmatrix}0&{}-i \\ -i &{}0\end{bmatrix} \begin{bmatrix}1&{}\dfrac{ie^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)}}{ \hat{r}_+(\lambda )f_+^2(\lambda )}&{} \\ 0&{}1 \end{bmatrix} \end{aligned}$$

and on $\Sigma _2$ as

$$\begin{aligned}&\begin{bmatrix} e^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )}&{}i \\ 0 &{} e^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix}\\&\quad \quad =\begin{bmatrix}1&{}0\\ i\dfrac{f_-^2 (\lambda )}{\hat{r}_-(\lambda )}e^{-x\left(g_+(\lambda ) - g_-(\lambda )\right)} &{}1 \end{bmatrix} \begin{bmatrix}0&{}i\\ i&{}0\end{bmatrix} \begin{bmatrix}1&{}0\\ -i\dfrac{f_+^2(\lambda )}{ \hat{r}_+(\lambda )}e^{x\left(g_+(\lambda ) - g_-(\lambda )\right)}&{} &{}1 \end{bmatrix}\, . \end{aligned}$$

We can now proceed with “opening lenses". We define a new vector function S as follows

$$\begin{aligned} S(\lambda ) = {\left\{ \begin{array}{ll} \displaystyle T(\lambda )\begin{bmatrix}1&{} \dfrac{-i}{ \hat{r}(\lambda )f^2(\lambda )}e^{-2x(g(\lambda )-\lambda )}&{} \\ 0&{}1 \end{bmatrix} &{}\quad \text {in the upper lens, above } \Sigma _{1} \\ \displaystyle T(\lambda )\begin{bmatrix}1&{}\dfrac{-i}{\hat{r}(\lambda )f^2 (\lambda )}e^{-2x\left(g(\lambda ) - \lambda \right)}\\ 0 &{}1 \end{bmatrix} &{}\quad \text {in the lower lens, below } \Sigma _{1} \\ \displaystyle T(\lambda )\begin{bmatrix}1&{}0\\ i\dfrac{f^2(\lambda )}{ \hat{r}(\lambda )}e^{2x\left(g(\lambda )-\lambda \right)}&{} &{}1 \end{bmatrix} &{}\quad \text {in the upper lens, above } \Sigma _{2} \\ \displaystyle T(\lambda ) \begin{bmatrix}1&{}0\\ i\dfrac{f^2 (\lambda )}{\hat{r}(\lambda )}e^{2x\left( g(\lambda )-\lambda \right)} &{}1 \end{bmatrix} &{}\quad \text {in the lower lens, below } \Sigma _{2}\\ \displaystyle T(\lambda ) &{}\quad \text {outside the lenses}\,. \end{array}\right. } \nonumber \\ \end{aligned}$$

(3.30)

The vector $S(\lambda )$ satisfies

$$\begin{aligned} \begin{aligned} S_+(\lambda )&=S_-(\lambda )V_S(\lambda ),\\ S(\lambda )&=\begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned} \end{aligned}$$

(3.31)

where the matrix $V_S$ for the jumps of $S(\lambda )$ is depicted in Fig.  3. In order to proceed we need the following lemma

Lemma 3.1

The following inequalities are satisfied

$$\begin{aligned}&{\text {Re}} \left(g(\lambda )-\lambda \right) <0\, ,\quad \lambda \in {\mathcal C}_1\backslash \{\eta _1,\eta _2\} \end{aligned}$$

(3.32)

$$\begin{aligned}&{\text {Re}} \left(g(\lambda )-\lambda \right) >0\, ,\quad \lambda \in {\mathcal C}_2\backslash \{-\eta _1,-\eta _2\}\, , \end{aligned}$$

(3.33)

where $\mathcal {C}_1$ and $\mathcal {C}_2$ are the contours defining the lenses as shown in Fig.  3.

Proof

Given $\lambda = x+iy$, we write $g_+(\lambda )-\lambda = u(x,y)+ iv(x,y)$. From the formula (3.17) for g, it follows that $g_+(\lambda )-\lambda$ is purely imaginary on $\Sigma _1 \cup \Sigma _2$; furthermore, for $\lambda \in \Sigma _1$

$$\begin{aligned} v_x = {\text {Im}} \left(g'_+(\lambda )-1\right) = \frac{\lambda ^2+\kappa }{\left|R_+(\lambda )\right|} = \frac{ \Omega }{\left|R_+(\lambda )\right|} \int _{-\eta _1}^{\eta _1} \frac{ \lambda ^2- \zeta ^2 }{R(\zeta )} \frac{\mathrm{d}\zeta }{2\pi i} > 0 \ . \end{aligned}$$

(3.34)

Using the Cauchy–Riemann equation it follows that $u_y = - v_x <0$ for $\lambda \in \Sigma _1$ and thus ${\text {Re}}\left( g(\lambda )-\lambda \right) <0$ for $\lambda$ above $\Sigma _1$ and $\lambda \in {\mathcal C}_1$. Repeating the same reasoning for the function $g_-(\lambda )-\lambda$ we obtain that ${\text {Re}}\left( g(\lambda )-\lambda \right) <0$ for $\lambda$ below $\Sigma _1$ and $\lambda \in {\mathcal C}_1$. In a similar way the inequality (3.33) can be obtained.

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 \rightarrow -\infty$, therefore those jump matrices are asymptotically close to the identity outside small neighbourhoods of $\pm \eta _1$ and $\pm \eta _2$. We are left with the model problem

$$\begin{aligned} S^{\infty }_+(\lambda )= & {} S^{\infty }_-(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} e^{x\Omega +\Delta } &{}0 \\ 0 &{} e^{-x\Omega -\Delta }\end{bmatrix} &{} \lambda \in [-\eta _1,\eta _1] \\ \begin{bmatrix}0 &{} -i\\ -i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{1}\\ \begin{bmatrix}0 &{} i\\ i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{2}\\ \end{array}\right. } \end{aligned}$$

(3.35)

$$\begin{aligned} S^{\infty }(\lambda )= & {} \begin{bmatrix} 1&1\end{bmatrix}+\mathcal {O}\left(\frac{1}{\lambda }\right),\quad \lambda \rightarrow \infty \,. \end{aligned}$$

(3.36)

The Riemann–Hilbert problem for $S^{\infty }$ has previously appeared in the study of long time asymptotics for KdV with step-like initial data [EGKT13]. Below we follow the lines of [EGKT13] to obtain the solution.

Fig. 3

Riemann–Hilbert problem for $S(\lambda )$ defined in (3.30). Opening lenses: the entries in gray in the jump matrices on the contours ${\mathcal C}_1$ and ${\mathcal C}_2$ are exponentially small in the regime as $x \rightarrow -\infty$

3.3 The outer parametrix $S^{\infty }$

To solve the Riemann–Hilbert problem (3.35) and (3.36) we introduce a two-sheeted Riemann surface $\mathfrak {X}$ of genus 1 associated to the multivalued function $R(\lambda )$, namely

$$\begin{aligned} \mathfrak {X}=\left\rbrace (\lambda ,\eta )\in \mathbb {C}^2\;|\; \eta ^2=R^{2}(\lambda )=(\lambda ^2-\eta _1^2)(\lambda ^2-\eta _2^2)\right\lbrace \, . \end{aligned}$$

The first sheet of the surface is identified with the sheet where $R(\lambda )$ is real and positive for $\lambda \in (\eta _2,+\infty )$. We introduce a canonical homology basis with the B cycle encircling $\Sigma _1$ clockwise on the first sheet and the A cycle going from $\Sigma _2$ to $\Sigma _1$ on the first sheet and coming back to $\Sigma _2$ on the second sheet. The points at infinity on the surface are denoted by $\infty ^{\pm }$ where $\infty ^+$ is on the first sheet and $\infty ^-$ on the second sheet of $\mathfrak {X}$. See Fig.  4. We introduce the holomorphic differential

$$\begin{aligned} \omega = \frac{\Omega }{R(\lambda )} \frac{\mathrm{d}\lambda }{4\pi i} \end{aligned}$$

(3.37)

so that

$$\begin{aligned} \oint _{A}\omega =1 \, . \end{aligned}$$

We also have

$$\begin{aligned} \tau =\oint _{B}\omega =\frac{i}{2}\dfrac{K(\sqrt{1-m^2})}{K(m)}\, ,\quad m=\dfrac{\eta _1}{\eta _2} \, . \end{aligned}$$

Next, we introduce the Jacobi elliptic function

$$\begin{aligned} \vartheta _3 (z;\tau )= \sum _{n\in \mathbb {Z}} e^{2\pi i \, nz + \pi n^2 i \tau } \ , \qquad z \in \mathbb {C} \ , \end{aligned}$$

(3.38)

which is an even function of z and satisfies the periodicity conditions

$$\begin{aligned} \vartheta _3 (z+h+k\tau ;\tau )=e^{-\pi i k^2 \tau -2\pi i kz}\vartheta _3 (z;\tau )\, ,\quad h,k\in \mathbb {Z}\, . \end{aligned}$$

(3.39)

We also recall that the Jacobi elliptic function with half-period ratio $\tau$ vanishes on the half period $\frac{\tau }{2}+\frac{1}{2}$. Finally, we define the integral

$$\begin{aligned} w(\lambda )=\int _{\eta _2}^{\lambda } \omega \end{aligned}$$

(3.40)

and we observe that

$$\begin{aligned} w(+\infty )=-\frac{1}{4},\quad w_+(\eta _1)=-\frac{\tau }{2},\quad w_+(-\eta _1)=-\frac{\tau }{2}-\frac{1}{2}\, \ , \end{aligned}$$

(3.41)

and

$$\begin{aligned} w(-\lambda ) = -w(\lambda ) - 1/2 \ \text{ for } \lambda \in \mathbb {C} \setminus \mathbb {R} \ . \end{aligned}$$

(3.42)

We introduce the following functions

$$\begin{aligned} \psi _{1}(\lambda )= & {} \frac{\vartheta _3 \left(2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3 \left(2w(\lambda ) -\frac{1}{2};2\tau \right)}\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3( \frac{x\Omega +\Delta }{2\pi i};2\tau )}\, ,\\ \psi _{2}(\lambda )= & {} \frac{\vartheta _3 \left(-2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2} ;2\tau \right)}{\vartheta _3 \left(-2w(\lambda ) -\frac{1}{2};2\tau \right)}\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3( \frac{x\Omega +\Delta }{2\pi i};2\tau )}\,, \end{aligned}$$

and we observe that

$$\begin{aligned} {\left\{ \begin{array}{ll} \vartheta _3 \left(\pm 2w_+(\eta _1) -\frac{1}{2};2\tau \right)=\vartheta _3 \left({\mp }\tau -\frac{1}{2};2\tau \right)=0,\\ \vartheta _3 \left(\pm 2w_+(-\eta _1) -\frac{1}{2};2\tau \right)=\vartheta _3 \left({\mp }\tau {\mp } 1 -\frac{1}{2};2\tau \right)=0. \end{array}\right. } \end{aligned}$$

It follows that the functions $\psi _1$ and $\psi _2$ are analytic except at $\lambda =\pm \eta _1$ where they admit, at most, square root singularities. Furthermore, the following jump relations are satisfied:

$$\begin{aligned}&w_+(\lambda )-w_-(\lambda )=0\quad \lambda \in [\eta _2,+\infty ) \end{aligned}$$

(3.43)

$$\begin{aligned}&w_+(\lambda )+w_-(\lambda )=0,\quad \lambda \in \Sigma _1\end{aligned}$$

(3.44)

$$\begin{aligned}&w_+(\lambda )-w_-(\lambda )=-\tau ,\quad \lambda \in (-\eta _1,\eta _1)\end{aligned}$$

(3.45)

$$\begin{aligned}&w_+(\lambda )+w_-(\lambda )=-1,\quad \lambda \in \Sigma _2 \, . \end{aligned}$$

(3.46)

Therefore for $\lambda \in \Sigma _1\cup \Sigma _2$ we have

$$\begin{aligned} \psi _{1+}(\lambda )=\psi _{2-}(\lambda )\, ,\quad \psi _{2+}(\lambda )=\psi _{1-}(\lambda )\, , \end{aligned}$$

(3.47)

while for $\lambda \in (-\eta _1,\eta _1)$

$$\begin{aligned} \begin{aligned}&\psi _{1+}(\lambda )=\psi _{1-}(\lambda )e^{x\Omega +\Delta }\, ,\quad \psi _{2+}(\lambda )=\psi _{2-}(\lambda )e^{-x\Omega -\Delta }\,. \end{aligned} \end{aligned}$$

(3.48)

Next we introduce the quantity

$$\begin{aligned} \gamma (\lambda )=\left( \dfrac{\lambda ^2-\eta _1^2}{\lambda ^2-\eta _2^2}\right) ^{\frac{1}{4}}\,, \end{aligned}$$

analytic in $\mathbb {C}\setminus (\Sigma _1 \cup \Sigma _2)$ and normalized such that $\gamma (\lambda ) \rightarrow 1$ as $\lambda \rightarrow \infty$. Then,

$$\begin{aligned} \gamma _+(\lambda )=-i\gamma _-(\lambda )\, , \ \ \text {for } \lambda \in \Sigma _1 \quad \text {and} \quad \gamma _+(\lambda )=i\gamma _-(\lambda )\, ,\ \ \text {for }\in \Sigma _2 \,. \end{aligned}$$

(3.49)

We are now ready to construct the solution of the Riemann–Hilbert problem (3.35)–(3.36).

Theorem 3.2

The vector $S^{\infty }(\lambda )$ given by

$$\begin{aligned} \begin{aligned} S^{\infty }(\lambda )=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3\left( \frac{x\Omega +\Delta }{2\pi i};2\tau \right)} \begin{bmatrix} \displaystyle \frac{\vartheta _3 \left(2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3 \left(2w(\lambda ) -\frac{1}{2};2\tau \right)}&\displaystyle \frac{\vartheta _3 \left(-2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3 \left(-2w(\lambda ) -\frac{1}{2};2\tau \right)} \end{bmatrix} \end{aligned} \end{aligned}$$

(3.50)

solves the Riemann–Hilbert problem (3.35).

Proof

We observe that $S^{\infty }(\lambda )$ has at most fourth root singularities at the branch points and it is bounded everywhere else on the complex plane. Because of (3.39) and (3.41) we have $S^{\infty }(\infty )=\begin{bmatrix} 1&1\end{bmatrix}$, namely the condition (3.36) is satisfied. Combining (3.47), (3.48) and (3.49), we conclude that the jump conditions (3.35) are satisfied.

Fig. 4

Construction of the genus-1 Riemann surface $\mathfrak {X}$ and its basis of cycles

This vector solution provides the asymptotic behaviour of the solution S to Riemann–Hilbert problem depicted in Fig. 3, for all $\lambda$ 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^{\infty }(\lambda )$. The matrix solution we construct has a pole at $\lambda =0$, however this pole does not affect the vector behaviour of our local and outer parametrices.

This will be accomplished in the next two subsections, by creating a second, independent vector solution.

3.4 The outer matrix parametrix $P^{\infty }$

Consider the 1-form $\mathrm{d}p(\lambda )=(1-g'(\lambda ))\ \mathrm{d}\lambda =\dfrac{\lambda ^2+\kappa }{R(\lambda )}\mathrm{d}\lambda$ and the Abelian integral

$$\begin{aligned} p(\lambda )=\int _{\eta _2}^{\lambda }\mathrm{d}p,\ \end{aligned}$$

(3.51)

which satisfies the relations

$$\begin{aligned}&p_+(\lambda )+p_-(\lambda )=0\quad \lambda \in \Sigma _1\cup \Sigma _2 \ , \end{aligned}$$

(3.52)

$$\begin{aligned}&p_+(\lambda )-p_-(\lambda )=-\Omega \quad \lambda \in (-\eta _1,\eta _1)\ , \end{aligned}$$

(3.53)

and for $\lambda \in \mathbb {C} \setminus \mathbb {R}$,

$$\begin{aligned} p(-\lambda ) =- p(\lambda ) \ . \end{aligned}$$

(3.54)

Then the vector function

$$\begin{aligned} \Psi := \begin{bmatrix} \varphi _1&\varphi _2 \end{bmatrix} = \begin{bmatrix} S^{\infty }_1&S^{\infty }_2 \end{bmatrix} e^{xp(\lambda )\sigma _3} \end{aligned}$$

(3.55)

solves a Riemann–Hilbert problem with constant jumps (independent from x).

Indeed on $\Sigma _1\cup \Sigma _2$ we have

$$\begin{aligned} \Psi _+(\lambda )&=\begin{bmatrix} S^{\infty }_{1+}&S^{\infty }_{2+} \end{bmatrix} e^{xp_+(\lambda )\sigma _3}= \begin{bmatrix} S^{\infty }_{1-}&S^{\infty }_{2-} \end{bmatrix} \begin{bmatrix}0 &{} {\mp } i\\ {\mp } i &{} 0 \end{bmatrix} e^{xp_+(\lambda )\sigma _3}\\&=\Psi _-(\lambda )e^{-xp_-(\lambda )\sigma _3}\begin{bmatrix}0 &{} {\mp } i\\ {\mp } i &{} 0 \end{bmatrix} e^{xp_+(\lambda )\sigma _3}\\&=\Psi _-(\lambda )\begin{bmatrix}0 &{} {\mp } i\\ {\mp } i &{} 0 \end{bmatrix}, \end{aligned}$$

where the ${\mp }$ signs correspond to $\Sigma _1$ and $\Sigma _2$ respectively, and the last identity has been obtained using (3.52). On $(-\eta _1,\eta _1)$ we have

$$\begin{aligned} \Psi _+(\lambda )&=\begin{bmatrix} S^{\infty }_{1+}&S^{\infty }_{2+} \end{bmatrix} e^{xp_+(\lambda )\sigma _3}= \begin{bmatrix} S^{\infty }_{1-}&S^{\infty }_{2-} \end{bmatrix} e^{(x\Omega +\Delta )\sigma _3} e^{xp_+(\lambda )\sigma _3}\\&=\Psi _-(\lambda )e^{-xp_-(\lambda )\sigma _3} e^{(x\Omega +\Delta )\sigma _3} e^{xp_+(\lambda )\sigma _3}\\&=\Psi _-(\lambda )e^{\Delta \sigma _3}, \end{aligned}$$

where the last identity has been obtained from (3.53). Therefore, the x derivative of $\Psi$ in (3.55), namely

$$\begin{aligned} \Psi _x(\lambda )=\begin{bmatrix} \varphi _{1x}&\varphi _{2x} \end{bmatrix} = \begin{bmatrix} S^{\infty }_{1x}&S^{\infty }_{2x} \end{bmatrix} e^{xp(\lambda )\sigma _3}+ \begin{bmatrix}p(\lambda ) S^{\infty }_{1}&-p(\lambda )S^{\infty }_{2} \end{bmatrix} e^{xp(\lambda )\sigma _3}, \end{aligned}$$

has the same jumps on $(-\eta _2,\eta _2)$ as $\Psi (\lambda )$. For this reason we consider the matrix function [Min] (see also [CG09])

$$\begin{aligned} \Phi (\lambda ):&=\begin{bmatrix} \varphi _1(\lambda )&{}\varphi _2(\lambda )\\ \varphi _{1x}(\lambda )&{}\varphi _{2x}(\lambda ) \end{bmatrix}\end{aligned}$$

(3.56)

$$\begin{aligned}&=\begin{bmatrix} S^{\infty }_1(\lambda )&{}S^{\infty }_2(\lambda )\\ p(\lambda )S^{\infty }_1(\lambda )+S^{\infty }_{1x}(\lambda )&{}-p(\lambda )S^{\infty }_2(\lambda )+S^{\infty }_{2x}(\lambda ) \end{bmatrix}e^{xp(\lambda )\sigma _3}\end{aligned}$$

(3.57)

$$\begin{aligned}&=\dfrac{1}{2}\begin{bmatrix} 1&{}1\\ \lambda &{}-\lambda \end{bmatrix} \begin{bmatrix} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_1+\dfrac{1}{\lambda }S^{\infty }_{1x} &{}(1-\frac{p(\lambda )}{\lambda })S^{\infty }_2+\dfrac{1}{\lambda }S^{\infty }_{2x}\\ (1-\frac{p(\lambda )}{\lambda })S^{\infty }_1-\dfrac{1}{\lambda }S^{\infty }_{1x}&{} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_2-\dfrac{1}{\lambda }S^{\infty }_{2x} \end{bmatrix}e^{xp(\lambda )\sigma _3}\,, \end{aligned}$$

(3.58)

where the last expression is an algebraic manipulation that can be verified by performing the matrix multiplication. It follows that for $\lambda \ne 0$, the matrix function $\begin{bmatrix} 1&{}1\\ \lambda &{}-\lambda \end{bmatrix}^{-1}\Phi (\lambda )$ has the same jumps as the vector $\Psi (\lambda )$ on the interval $(-\eta _2,\eta _2)$ and the matrix function $\begin{bmatrix} 1&{}1\\ \lambda &{}-\lambda \end{bmatrix}^{-1}\Phi (\lambda )e^{-xp(\lambda )\sigma _3}$ has the same jumps as the vector $S^{\infty }$ defined in (3.50). For this reason we take as a matrix solution for the exterior parametrix

$$\begin{aligned} P^{\infty }(\lambda )=\frac{1}{2}\begin{bmatrix} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_1+\dfrac{1}{\lambda }S^{\infty }_{1x} &{}(1-\frac{p(\lambda )}{\lambda })S^{\infty }_2+\dfrac{1}{\lambda }S^{\infty }_{2x}\\ (1-\frac{p(\lambda )}{\lambda })S^{\infty }_1-\dfrac{1}{\lambda }S^{\infty }_{1x}&{} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_2-\dfrac{1}{\lambda }S^{\infty }_{2x} \end{bmatrix}\,. \end{aligned}$$

(3.59)

It satisfies the following Riemann–Hilbert problem:

$$\begin{aligned}&P^{\infty }(\lambda )\hbox { is analytic for } \lambda \in \mathbb {C}\backslash [-\eta _2,\eta _2]\hbox { with a singularity at }\lambda =0, \end{aligned}$$

(3.60)

$$\begin{aligned}&P^{\infty }_+(\lambda ) = P^{\infty }_-(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} e^{x\Omega +\Delta } &{}0 \\ 0 &{} e^{-x\Omega -\Delta }\end{bmatrix} &{} \lambda \in [-\eta _1,\eta _1] \\ \begin{bmatrix}0 &{} -i\\ -i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{1} \\ \begin{bmatrix}0 &{} i\\ i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{2}, \\ \end{array}\right. } \end{aligned}$$

(3.61)

$$\begin{aligned}&P^{\infty }(\lambda )=\begin{bmatrix} 1&{} 0\\ 0&{}1\end{bmatrix}+\mathcal {O}\left(\frac{1}{\lambda }\right),\quad \lambda \rightarrow \infty \,. \end{aligned}$$

(3.62)

Despite the singularity of the matrix $P^{\infty }(\lambda )$ at $\lambda =0$, its determinant is equal to one. Before proving this fact we first make a slight change of notation that will be relevant in the next sections. We observe that the x-derivative of $S_1^\infty$ and $S_2^\infty$ can be written in the form

$$\begin{aligned} S^{\infty }_{1x}(\lambda )=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3 \left(2w(\lambda ) -\frac{1}{2};2\tau \right)} \dfrac{\Omega }{2\pi i}\dfrac{\mathrm{d}}{\mathrm{d}z}\left[ \frac{\vartheta _3 \left(z+2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3\left(z+ \frac{x\Omega +\Delta }{2\pi i};2\tau \right)}\right] \bigg |_{z=0} \end{aligned}$$

and similarly for $S^{\infty }_2(\lambda )$. Since in the next sections we will use similar formulas where the quantity $\Omega$ is replaced by $\tilde{\Omega }$ that is dependent on x and t, it is important to distinguish the operation of derivative with respect to x from the operation on the right hand side of the above expression. For this reason we introduce the notation

$$\begin{aligned} \nabla _{\Omega }S^{\infty }_1(\lambda ):=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3 \left(2w(\lambda ) -\frac{1}{2};2\tau \right)} \dfrac{\Omega }{2\pi i}\dfrac{\mathrm{d}}{\mathrm{d}z}\left[ \frac{\vartheta _3 \left(z+2w(\lambda ) +\frac{x\Omega +\Delta }{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3\left(z+ \frac{x\Omega +\Delta }{2\pi i};2\tau \right)}\right] \bigg |_{z=0} \end{aligned}$$

and similarly for $S^{\infty }_2(\lambda )$. Clearly, when $\Omega$ is x-independent then $\nabla _{\Omega }S^{\infty }(\lambda )\equiv S_{x}(\lambda )$. Therefore the exterior parametrix $P^{\infty }$ in (3.59) will be written in the form

$$\begin{aligned} P^{\infty }(\lambda )=\frac{1}{2}\begin{bmatrix} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_1(\lambda )+\dfrac{1}{\lambda }\nabla _{\Omega }S^{\infty }_1(\lambda ) &{}(1-\frac{p(\lambda )}{\lambda })S^{\infty }_2(\lambda )+\dfrac{1}{\lambda }\nabla _{\Omega }S^{\infty }_2(\lambda )\\ (1-\frac{p(\lambda )}{\lambda })S^{\infty }_1(\lambda )-\dfrac{1}{\lambda }\nabla _{\Omega }S^{\infty }_1(\lambda )&{} (1+\frac{p(\lambda )}{\lambda })S^{\infty }_2(\lambda )-\dfrac{1}{\lambda }\nabla _{\Omega }S^{\infty }_2(\lambda ) \end{bmatrix}\,. \end{aligned}$$

(3.63)

Lemma 3.3

We have

$$\begin{aligned} \det P^{\infty }(\lambda ) \equiv 1\,. \end{aligned}$$

(3.64)

Proof

We observe that

$$\begin{aligned} \det P^{\infty }(\lambda )=-\dfrac{1}{2\lambda }\Big ( -2p(\lambda )S^{\infty }_2(\lambda )S^{\infty }_1(\lambda )+S^{\infty }_1(\lambda ) \nabla _{\Omega } S^{\infty }_{2}(\lambda )-\nabla _{\Omega }S^{\infty }_{1}(\lambda )S^{\infty }_2(\lambda ) \Big ) \end{aligned}$$

(3.65)

does not have any jumps on the complex plane and therefore it is a meromorphic function on the complex plane. Considering the behaviour near $\lambda =\eta _2$, we have

$$\begin{aligned} S^{\infty }_1(\lambda )=\dfrac{1}{(\lambda -\eta _2)^{\frac{1}{4}}}(\sum _{k=0}^{\infty }\gamma _k(\lambda -\eta _2)^k)\left( \sum _{j=0}^{\infty }S_{j}^{(1)}(\lambda -\eta _2)^j+\sqrt{\lambda -\eta _2}\sum _{j=0}^{\infty }S_j^{(2)}(\lambda -\eta _2)^j\right) \end{aligned}$$

where $\gamma _k$ are the coefficients of the Puiseux expansion of $\gamma (\lambda )$ and $S_{j}^{(1)}$ and $S_{j}^{(2)}$ are the coefficients of the Puiseux expansion of the $\vartheta _3$ function terms of $S^{\infty }_1(\lambda )$ near $\lambda =\eta _2$; in particular, $\gamma _0\ne 0$ and $S_{0}^{(1,2)}\ne 0$. In a similar way we obtain

$$\begin{aligned} S^{\infty }_2(\lambda )=\dfrac{1}{(\lambda -\eta _2)^{\frac{1}{4}}}(\sum _{k=0}^{\infty }\gamma _k(\lambda -\eta _2)^k)\left( \sum _{j=0}^{\infty }S_{j}^{(1)}(\lambda -\eta _2)^j-\sqrt{\lambda -\eta _2}\sum _{j=0}^{\infty }S_j^{(2)}(\lambda -\eta _2)^j\right) \end{aligned}$$

and

$$\begin{aligned} p(\lambda )=2\sqrt{\lambda -\eta _2}\sum _{k=0}^{\infty }c_k(\lambda -\eta _2)^k. \end{aligned}$$

Plugging the above three expansions into (3.65) it is straightforward to check that $\det P^{\infty }(\lambda )$ has a Taylor expansion at the point $\lambda =\eta _2$. It can be checked similarly that $\det P^{\infty }(\lambda )$ has a Taylor expansion at $\lambda =\pm \eta _1$ and $\lambda =-\eta _2$. Regarding the point $\lambda =0$, we consider the Abelian integral $p(\lambda )$ defined in (3.51) and denote by $p_{\pm }(\lambda )$ the boundary values of $p(\lambda )$ on the real axis. We have

$$\begin{aligned} p_\pm (0)=\int _{\eta _2}^0 \mathrm{d}p_{\pm }(\xi )=\int _{\eta _2}^{\eta _1}\mathrm{d}p_{\pm }(\xi )+\int _{\eta _1}^0 \mathrm{d}p(\xi )={\mp }\frac{\Omega }{2}, \end{aligned}$$

(3.66)

where, in the last relation we use the identity

$$\begin{aligned} 0=\int _{-\eta _1}^{\eta _1}\mathrm{d}p(\xi )=\int _{-\eta _1}^{0}\mathrm{d}p(\xi )+\int _{0}^{\eta _1}\mathrm{d}p(\xi )=\int _{\eta _1}^{0}\mathrm{d}p(-\xi )+\int _{0}^{\eta _1}\mathrm{d}p(\xi )= 2\int _{0}^{\eta _1}\mathrm{d}p(\xi ). \end{aligned}$$

In this last line we do not use $\mathrm{d}p_\pm (\lambda )$ because the function $R(\lambda )$ is analytic off the contours $\Sigma _1$ and $\Sigma _2$ and the same property holds for $\mathrm{d}p$. Using the periodicity properties of the Jacobi elliptic function

$$\begin{aligned} \vartheta _3 (z+h+k\tau ;\tau )=e^{-\pi i k^2 \tau -2\pi i kz}\vartheta _3 (z;\tau )\, ,\quad h,k\in \mathbb {Z}\,, \end{aligned}$$

we have

$$\begin{aligned} S^{\infty }_{\pm }(0)= \gamma (0)\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3 \left( \pm \tau ;2\tau \right)} \frac{\vartheta _3 \left(\pm \tau +\frac{x\Omega +\Delta }{2\pi i};2\tau \right)}{ \vartheta _3( \frac{x\Omega +\Delta }{2\pi i};2\tau )} \begin{bmatrix} \displaystyle e^{ \pm (x\Omega +\Delta )}&1\end{bmatrix} \end{aligned}$$

(3.67)

and

$$\begin{aligned} \begin{aligned} \nabla _{\Omega }S^{\infty }_{\pm }(0)&= \dfrac{\Omega }{2\pi i}S^{\infty }_{\pm }(0)\left. \partial _z\left[ \log \frac{\vartheta _3 \left(z\pm \tau +\frac{x\Omega +\Delta }{2\pi i};2\tau \right)}{\vartheta _3 \left(z+\frac{x\Omega +\Delta }{2\pi i};2\tau \right) } \right] \right| _{{z=0}} \, +\begin{bmatrix}\pm \Omega S_{1\pm }^{\infty }(0)\;0\end{bmatrix}.\end{aligned} \end{aligned}$$

(3.68)

We conclude that

$$\begin{aligned}&\Big ( -2p(\lambda )S^{\infty }_2(\lambda )S^{\infty }_1(\lambda )+S^{\infty }_1(\lambda )\nabla _{\Omega }S^{\infty }_{2}(\lambda )-\nabla _{\Omega }S^{\infty }_{1}(\lambda )S^{\infty }_2(\lambda ) \Big )_{\pm }\\&\quad =\pm \Omega S^{\infty }_{2\pm }(0)S^{\infty }_{1\pm }(0) \\&\quad +S^{\infty }_{1\pm }(0) \dfrac{\Omega }{2\pi i}S^{\infty }_{2\pm }(0)\left. \partial _z\left[ \log \frac{\vartheta _3 \left(z\pm \tau +\frac{x\Omega +\Delta }{2\pi i};2\tau \right)}{\vartheta _3 \left(z+\frac{x\Omega +\Delta }{2\pi i};2\tau \right) } \right] \right| _{{z=0}}\\&\quad -\dfrac{\Omega }{2\pi i}S^{\infty }_{1\pm }(0)S^{\infty }_{2\pm }(0)\left. \partial _z\left[ \log \frac{\vartheta _3 \left(z\pm \tau +\frac{x\Omega +\Delta }{2\pi i};2\tau \right)}{\vartheta _3 \left(z+\frac{x\Omega +\Delta }{2\pi i};2\tau \right) } \right] \right| _{{z=0}}\\&\quad {\mp }\Omega S ^{\infty }_{1\pm }(0)S^{\infty }_{2\pm }(0)+O(\lambda )=O(\lambda ) \end{aligned}$$

as $\lambda \rightarrow 0$. Therefore

$$\begin{aligned} \det P^{\infty }(\lambda )=-\dfrac{1}{2\lambda } \Big ( -2p(\lambda )S^{\infty }_2(\lambda )S^{\infty }_1(\lambda )+S^{\infty }_1(\lambda )\nabla _{\Omega }S^{\infty }_{2}(\lambda )-\nabla _{\Omega }S^{\infty }_{1}(\lambda )S^{\infty }_2(\lambda ) \Big ) \end{aligned}$$

is a holomorphic function of $\lambda$ near $\lambda =0$. Since

$$\begin{aligned} \det P^{\infty }(\lambda )=1+O(\lambda ^{-1}),\quad \text{ as } \lambda \rightarrow \infty , \end{aligned}$$

it follows by Liouville’s theorem that

$$\begin{aligned} \det P^{\infty }(\lambda ) \equiv 1\,. \end{aligned}$$

Remark: The reader may verify using the definition (3.63), along with the symmetry relations (3.42) and (3.54), that $P^{\infty }$ satisfies the symmetry

$$\begin{aligned} P^{\infty }(-\lambda ) = \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} P^{\infty }(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \ . \end{aligned}$$

(3.69)

3.5 The local parametrix $P^{\pm \eta _j}$ at the endpoints

Thanks to Lemma 3.1, the off-diagonal entries of the jump matrices for S exponentially vanish as $x\rightarrow -\infty$ along the upper and lower lenses, while near the endpoints the g function has a square-root-vanishing behaviour

$$\begin{aligned} g_+(\lambda ) -g_-(\lambda ) = \mathcal {O}\left( \sqrt{\lambda {\mp } \eta _2}\right) \qquad \text {as } \lambda \rightarrow \pm \eta _2 \ , \ \end{aligned}$$

(3.70)

and

$$\begin{aligned} g_+(\lambda ) -g_-(\lambda ) - \Omega = \mathcal {O}\left( \sqrt{\lambda {\mp } \eta _1}\right) \qquad \text {as } \lambda \rightarrow \pm \eta _1 \ . \ \end{aligned}$$

(3.71)

Additionally, the original solution Y of the Riemann–Hilbert problem (3.4)–(3.6) has a logarithmic singularity in those points. Therefore, the jump matrices for S are bounded in a neighbourhood of those points (but they are not close to the identity).

On the other hand, the outer parametrix $P^{\infty }$ is a good approximation of the solution S to the Riemann–Hilbert problem away from the endpoints $\lambda = \pm \eta _2, \pm \eta _1$, where $P^{\infty }$ exhibits a fourth-root singularity. So, we need to introduce four local parametrices $P^{\pm \eta _j}$ ($j=1,2$) in a suitable neighbourhood of each endpoint.

3.5.1 Local parametrix near $\lambda = \eta _2$.

We show here the construction of a (matrix) local parametrix $P^{\eta _2}$ around $\lambda = \eta _2$.

Performing the same calculations as in [KMAV04, Section 6], we will construct a local parametrix $P^{\eta _2}$ with the help of modified Bessel functions. We fix a small disc $B^{(\eta _2)}_{\rho } = \left\rbrace \lambda \in \mathbb {C} \left| \, \left|\lambda - \eta _2\right|< \rho \right. \right\lbrace$ centered at $\eta _2$ of radius $\rho$, and we define the (local) conformal map

$$\begin{aligned} \zeta = \frac{1}{4} \left[ x \left(g(\lambda ) - \lambda \right)\right]^2, \qquad \lambda \in B^{(\eta _2)}_\rho \ . \end{aligned}$$

(3.72)

To define the local parametrix $P^{\eta _{2}}$ in $B^{(\eta _2)}_\rho$, we consider

$$\begin{aligned} P(\lambda ) = S(\lambda ) \left( \frac{ e^{ i \pi /4}}{\sqrt{\pm \hat{r}} f} \right) ^{\sigma _{3}}&\lambda \in B^{(\eta _2)}_{\rho } \cap \mathbb {C}_{\pm } , \end{aligned}$$

and then, using the inverse of the transformation $\zeta (\lambda )$, we define

$$\begin{aligned} P^{(1)}(\zeta ) = P(\lambda (\zeta )) e^{-2 \zeta ^{\frac{1}{2}} \sigma _3} \begin{bmatrix}0&{}1\\ 1&{}0 \end{bmatrix}\ , \quad \zeta \in \mathbb {C}\, , \end{aligned}$$

with branch cut $(-\infty ,0]$. By construction, $P^{(1)}$ satisfies a Riemann–Hilbert problem with jumps

$$\begin{aligned} P^{(1)}_+(\zeta ) = P^{(1)}_-(\zeta ) {\left\{ \begin{array}{ll} \begin{bmatrix} 1 &{} 0 \\ 1 &{} 1 \end{bmatrix} &{} \text {on \{upper and lower lenses\}}\cap B^{(\eta _2)}_\rho \\ \begin{bmatrix} 0 &{} 1 \\ -1\ &{} 0 \end{bmatrix}&\text {on } (-\infty ,0] \cap B^{(\eta _2)}_\rho \ . \end{array}\right. } \end{aligned}$$

(3.73)

We introduce now the model parametrix $\Psi _\mathrm {Bes}(\zeta )$ as in [KMAV04, formulæ (6.16)–(6.20)]). The Riemann–Hilbert problem for $\Psi _\mathrm {Bes}$ is the following:

• (a) $\Psi _\mathrm {Bes}$ is analytic for $\zeta \in \mathbb {C} \backslash \Gamma _{\Psi }$, where $\Gamma _\Psi$ is the union of the three contours $\Gamma _\pm = \left\rbrace \arg \zeta = \pm \frac{2\pi }{3}\right\lbrace$ and $\Gamma _0 = \left\rbrace \arg \zeta = \pi \right\lbrace$;

• (b) $\Psi _{\mathrm {Bes}}$ satisfies the following jump relations

  $$\begin{aligned} \Psi _{\mathrm {Bes}\, +}(\zeta ) = \Psi _{\mathrm {Bes}\, -}(\zeta ) {\left\{ \begin{array}{ll} \begin{bmatrix} 1 &{} 0 \\ 1 &{} 1 \end{bmatrix} &{} \text {on } \Gamma _+ \cup \Gamma _-\\ \begin{bmatrix} 0 &{} 1 \\ -1\ &{} 0 \end{bmatrix}&\text {on }\Gamma _0, \end{array}\right. }\, \end{aligned}$$

  (3.74)

• (c) as $\zeta \rightarrow 0$

  $$\begin{aligned} \Psi _\mathrm {Bes}(\zeta ) = \begin{bmatrix} \mathcal {O}\left( \ln |\zeta | \right) &{} \mathcal {O}\left( \ln |\zeta | \right) \\ \mathcal {O} \left( \ln |\zeta | \right)&{} \mathcal {O} \left(\ln |\zeta |\right) \end{bmatrix} \ . \end{aligned}$$

  (3.75)

The solution is the following

$$\begin{aligned} \Psi _\mathrm {Bes}(\zeta ) = {\left\{ \begin{array}{ll} \begin{bmatrix} \displaystyle I_0 (2\zeta ^{\frac{1}{2}}) &{}\displaystyle \frac{i}{\pi } K_0 (2\zeta ^{\frac{1}{2}})\\ \displaystyle 2\pi i \zeta ^{\frac{1}{2}}I'_0 (2\zeta ^{\frac{1}{2}}) &{} \displaystyle - 2\zeta ^{\frac{1}{2}}K_0(2\zeta ^{\frac{1}{2}}) \end{bmatrix} &{} |\arg \zeta |< \frac{2\pi }{3} \\ &{} \\ \begin{bmatrix} \displaystyle \frac{1}{2}H^{(1)}_0 (2(-\zeta )^{\frac{1}{2}}) &{}\displaystyle \frac{1}{2} H^{(2)}_0 (2(-\zeta )^{\frac{1}{2}})\\ \displaystyle \pi \zeta ^{\frac{1}{2}} \left[H^{(1)}_0 (2(-\zeta )^{\frac{1}{2}})\right]' &{} \displaystyle \pi \zeta ^{\frac{1}{2}} \left[H^{(2)}_0(2(-\zeta )^{\frac{1}{2}})\right]' \end{bmatrix} &{} \frac{2\pi }{3}< |\arg \zeta |< \pi \\ &{} \\ \begin{bmatrix} \displaystyle \frac{1}{2}H^{(2)}_0 (2(-\zeta )^{\frac{1}{2}}) &{}\displaystyle -\frac{1}{2} H^{(1)}_0 (2(-\zeta )^{\frac{1}{2}})\\ \displaystyle -\pi \zeta ^{\frac{1}{2}} \left[H^{(2)}_0 (2(-\zeta )^{\frac{1}{2}})\right]' &{} \displaystyle \pi \zeta ^{\frac{1}{2}} \left[H^{(1)}_0(2(-\zeta )^{\frac{1}{2}})\right]' \end{bmatrix}&-\pi< |\arg \zeta |< - \frac{2\pi }{3} \end{array}\right. } \nonumber \\ \end{aligned}$$

(3.76)

with asymptotic behaviour at infinity

$$\begin{aligned} \Psi _\mathrm {Bes}(\zeta ) = \left( 2\pi \zeta ^{\frac{1}{2}}\right)^{-\frac{1}{2}\sigma _3} \frac{1}{\sqrt{2}} \begin{bmatrix} \displaystyle 1&{} i \displaystyle \\ \displaystyle i&{} \displaystyle 1 \end{bmatrix} \left(I + \mathcal {O} \left(\frac{1}{\zeta ^{\frac{1}{2}}}\right)\right)e^{2\zeta ^{\frac{1}{2}}\sigma _3} \end{aligned}$$

(3.77)

uniformly as $\zeta \rightarrow \infty$ everywhere in the complex plane aside from the jumps.

In the above formulæ $I_0(\zeta )$, $K_0(\zeta )$ are the modified Bessel functions of first and second kind, respectively, and $H^{(j)}(\zeta )$ the Hankel functions.

In conclusion, the local parametrix around the endpoint $\lambda = \eta _2$ is

$$\begin{aligned} P^{\eta _{2}}(\lambda ) = A(\lambda ) \Psi _\mathrm {Bes}(\zeta (\lambda )) \begin{bmatrix} 0&{}1\\ 1&{}0\end{bmatrix} e^{2\zeta (\lambda )^{\frac{1}{2}}\sigma _3} \left( \frac{ e^{ i \pi /4}}{\sqrt{\pm \hat{r}(\lambda )} f(\lambda )} \right) ^{-\sigma _{3}} \quad \lambda \in B^{(\eta _2)}_{\rho } \cap \mathbb {C}_{\pm } \ ,\nonumber \\ \end{aligned}$$

(3.78)

where A is a prefactor that is determined by imposing that

$$\begin{aligned} P^{\eta _{2}}(\lambda ) \left(P^{\infty }(\lambda ) \right)^{-1} = I + \mathcal {O}\left(|x|^{-1}\right) \qquad \text {as } x \rightarrow - \infty , \ \text {for } \lambda \in \partial B^{(\eta _2)}_\rho \backslash \Sigma _\Psi \ .\nonumber \\ \end{aligned}$$

(3.79)

Therefore, we set

$$\begin{aligned} A(\lambda ) = P^{\infty }(\lambda ) \left( \frac{ e^{ i \pi /4}}{\sqrt{\pm \hat{r}(\lambda )} f(\lambda )} \right) ^{\sigma _{3}} \frac{1}{\sqrt{2}} \begin{bmatrix}-i&{} 1 \\ 1 &{} -i\end{bmatrix} \left(2\pi \zeta ^{\frac{1}{2}}\right)^{\frac{1}{2}\sigma _3} \quad \lambda \in B^{(\eta _2)}_{\rho } \cap \mathbb {C}_{\pm }\ .\nonumber \\ \end{aligned}$$

(3.80)

By construction, A is well-defined and analytic in a neighbourhood of $\eta _2$, minus the cut $(-\infty , \eta _2]$; additionally, it is easy to see that A is invertible ($\det A(\lambda ) \equiv 1$).

Lemma 3.4

$A(\lambda )$ is analytic everywhere in the neighbourhood $B^{(\eta _2)}_\rho$ of $\eta _2$.

Proof

To prove the statement, one needs to check that A has no jumps across the interval ${\Sigma _1} \cap B^{(\eta _2)}_{\rho }$ and that it has at most a removable singularity at $\lambda = \eta _2$. Starting from (3.80) we observe from (3.23) and (3.29) that for $\lambda \in \Sigma _1$, $\sqrt{ \hat{r}(\lambda )} f_+(\lambda ) = \left( \sqrt{-\hat{r}(\lambda )} f_-(\lambda ) \right) ^{-1}$. Using this and the jump (3.613.62) of $P^\infty$ on $\Sigma _1$ we have

$$\begin{aligned} A_+(\lambda )&= P_-^{\infty }(\lambda ) \begin{bmatrix} 0&-i\mathrm{i}&0\end{bmatrix} \left( \frac{ e^{ i \pi /4}}{\sqrt{\hat{r}_-(\lambda )} f_-(\lambda )} \right) ^{-\sigma _{3}} i^{\sigma _3}\frac{1}{\sqrt{2}} \begin{bmatrix}-i&{} 1 \\ 1 &{} -i\end{bmatrix} i^{\sigma _3} \left(2\pi \zeta _-^{\frac{1}{2}}\right)^{\frac{1}{2}\sigma _3} \\&= P_-^{\infty }(\lambda ) \left( \frac{ e^{ i \pi /4}}{\sqrt{\hat{r}_-(\lambda )} f_-(\lambda )} \right) ^{\sigma _{3}} \begin{bmatrix} 0&-i\mathrm{i}&0\end{bmatrix} i^{\sigma _3}\frac{1}{\sqrt{2}} \begin{bmatrix}-i&{} 1 \\ 1 &{} -i\end{bmatrix} i^{\sigma _3} \left(2\pi \zeta _-^{\frac{1}{2}}\right)^{\frac{1}{2}\sigma _3} \\&= P_-^{\infty }(\lambda ) \left( \frac{ e^{ i \pi /4}}{\sqrt{\hat{r}_-(\lambda )} f_-(\lambda )} \right) ^{\sigma _{3}} \frac{1}{\sqrt{2}} \begin{bmatrix}-i&{} 1 \\ 1 &{} -i\end{bmatrix} \left(2\pi \zeta _-^{\frac{1}{2}}\right)^{\frac{1}{2}\sigma _3} = A_-(\lambda ) \end{aligned}$$

Next, we notice that $\zeta (\lambda )$ has a simple zero at $\eta _2$ by construction, thus $\zeta (\lambda )^{\frac{1}{4}\sigma _3}$ has at most a fourth-root singularity at the point $\lambda = \eta _2$. Also the outer parametrix $P^{\infty }(\lambda )$ has at most a fourth-root singularity near $\eta _2$ and consequently all the entries of $A(\lambda )$ have at most a square root singularity at $\lambda = \eta _2$.

On the other hand $A(\lambda )$ is analytic in $B^{(\eta _2)}_\rho \backslash \{\eta _2\}$, therefore the point $\lambda = \eta _2$ is a removable singularity and $A(\lambda )$ is indeed analytic everywhere in $B^{(\eta _2)}_\rho$.

3.5.2 Local parametrix near other branch points.

The construction of the parametrix in a vicinity $B^{(\eta _1)}_\rho$ of $\eta _{1}$ is quite similar, and has also been carried out in [KMAV04, Section 6], so we will not present the formula here.

For the parametrices near $- \eta _{2}$ and $- \eta _{1}$, it will prove convenient to construct them explicitly via the underlying $\lambda \mapsto - \lambda$ symmetry, as follows:

$$\begin{aligned}&P^{-\eta _{2}} := \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} P^{\eta _{2}}(-\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}, \end{aligned}$$

(3.81)

$$\begin{aligned}&P^{-\eta _{1}} := \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} P^{\eta _{1}}(-\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}. \end{aligned}$$

(3.82)

First, the reader may verify that, if $P^{\eta _{j}}$ satisfies the appropriate jump relationships along the contours within the disc centered at $\eta _{j}$, then $P^{-\eta _{j}}$ satisfies the appropriate jump relationships along the contours within the disc (of the same radius) centered at $- \eta _{j}$. Along the way, the following symmetry relations are needed (and are easy to establish) for $\lambda \in \mathbb {C}\backslash (-\eta _2,\eta _2)$,

$$\begin{aligned}&\hat{r}(-\lambda ) = \hat{r}(\lambda ), \end{aligned}$$

(3.83)

$$\begin{aligned}&f^{2}(-\lambda ) = f^{-2}(\lambda ) ,\end{aligned}$$

(3.84)

$$\begin{aligned}&g(-\lambda ) = - g(\lambda ) \ . \end{aligned}$$

(3.85)

Moreover, since $P^{\eta _{j}}$ has been constructed to satisfy

$$\begin{aligned} P^{\eta _{j}}(\lambda )P^{\infty }(\lambda )^{-1} = I + \mathcal {O} \left( \frac{1}{x} \right) \quad {\text {as}}\, x\rightarrow -\infty , \end{aligned}$$

for $\lambda$ on the boundary of $B_{\rho }^{(\eta _{j})}$ (the small disc of radius $\rho$ centered at $\eta _{j}$), it follows that $P^{-\eta _{j}}$ satisfies

$$\begin{aligned} P^{-\eta _{j}}(\lambda )P^{\infty }(\lambda )^{-1} = I + \mathcal {O} \left( \frac{1}{x} \right) \quad {\text {as}}\, x\rightarrow -\infty , \end{aligned}$$

for $\lambda$ on the boundary of an analogous small disc $B^{(-\eta _j)}_{\rho }$ centered at $-\eta _{j}$.

3.6 Small norm argument and determination of u(x, 0) for large negative x

Define the error vector

$$\begin{aligned} \mathcal {E}(\lambda ) = S(\lambda ) \left( P(\lambda ) \right) ^{-1}. \end{aligned}$$

(3.86)

where the global parametrix $P(\lambda )$ is defined by

$$\begin{aligned} P(\lambda ) = {\left\{ \begin{array}{ll} P^{\infty }(\lambda ) &{} \lambda \in \mathbb {C}\backslash \cup _{j=1,2}B^{(\pm \eta _j)}_{\rho } \\ P^{\eta _{2}}(\lambda ) &{} \lambda \in B^{(\eta _2)}_{\rho } \\ P^{\eta _{1}} (\lambda ) &{} \lambda \in B^{(\eta _1)}_{\rho } \\ P^{-\eta _{1}}(\lambda ) &{} \lambda \in B^{(-\eta _1)}_{\rho } \\ P^{-\eta _{2}}(\lambda ) &{} \lambda \in B^{(-\eta _2)}_{\rho } \ . \end{array}\right. } \end{aligned}$$

(3.87)

Fig. 5

The Riemann–Hilbert problem for the remainder $\mathcal E$

Then across any contour where either S is non-analytic or any boundary in the definition of P the matrix $\mathcal {E}$ has a jump given by

$$\begin{aligned} \mathcal E_{+} (\lambda ) = \mathcal E_{-} (\lambda ) V_{\mathcal E}(\lambda ) \end{aligned}$$

(3.88)

with

$$\begin{aligned} V_{\mathcal E}(\lambda )= & {} \left( \mathcal E_{-} (\lambda ) \right) ^{-1} \mathcal E_{+} (\lambda ) = P_-(\lambda ) \left( S_-(\lambda ) \right) ^{-1} S_+(\lambda ) \left( P_+(\lambda ) \right) ^{-1} \nonumber \\= & {} P_-(\lambda ) V_S(\lambda ) \left( V_P(\lambda ) \right) ^{-1} \left( P_-(\lambda ) \right) ^{-1} \, \end{aligned}$$

(3.89)

where the jump matrix $V_S$ is as defined in Fig. 3, $V_P$ is the jump of P, and both jumps are understood to be the identity matrix anywhere S or P is analytic respectively. Observe that both within the discs $B^{(\pm \eta _j)}_{\rho }$, and across any component of $(-\eta _{2}, \eta _{2})$ outside the discs $B^{(\pm \eta _j)}_\rho$, $j=1,2$, the quantities $V_{S}$ and $V_P$ coincide, and hence $\mathcal E$ has no jump across those contours. Across the lens boundaries (outside the discs) we have $V_{P} = I$, and hence

$$\begin{aligned}&V_{\mathcal E}(\lambda )=\left( P^{\infty }(\lambda ) \right) V_{S}(\lambda ) \left( P^{\infty } (\lambda )\right) ^{-1} = \left( I + \mathcal {O} \left( e^{- c x } \right) \right) ,&\lambda \in {\mathcal C}_j,\, j=1,2, \end{aligned}$$

(3.90)

while across the circles centered at $\pm \eta _{j}$ (which we have chosen to orient counter-clockwise), we have

$$\begin{aligned}&V_{\mathcal E}(\lambda )=\left( P^{\infty }(\lambda ) \right) ^{-1} P^{\pm \eta _{j}}(\lambda ) = \left( I + \mathcal {O} \left( x^{-1} \right) \right) ,&\lambda \in \partial B^{(\pm \eta _j)}_\rho ,\;j=1,2. \end{aligned}$$

(3.91)

Finally, since $P = P^\infty (\lambda )$ near $\lambda =0$ and $P^\infty (\lambda )$ is singular there, we need to check the behaviour of $\mathcal {E}$ at $\lambda = 0$.

Lemma 3.5

The error vector $\mathcal {E}$ defined by (3.86)–(3.87) is regular at $\lambda =0$.

Proof

Near $\lambda =0$, $\mathcal {E}(\lambda ) = S(\lambda ) \left( P^{\infty }(\lambda ) \right) ^{-1}$. We have $S(\lambda )=Y(\lambda )e^{xg(\lambda )\sigma _3}f(\lambda )^{\sigma _3}$ with the functions $g(\lambda )$ and $f(\lambda )$ defined in (3.17) and (3.27) respectively and where $Y(\lambda )$ is the solution of the Riemann–Hilbert problem defined by (3.4)-(3.6) whose existence is established in the Appendix. We need to prove that

$$\begin{aligned} Y(\lambda )e^{xg(\lambda )\sigma _3}f(\lambda )^{\sigma _3}(P^{\infty }(\lambda ))^{-1} \end{aligned}$$

is regular at $\lambda =0$ where $Y(\lambda )$ satisfies the symmetry (3.6) so that $Y_1(0)=Y_2(0)$. We observe that $g(\lambda )=\lambda -p(\lambda )$ so that by (3.66)

$$\begin{aligned} g_\pm (0)=\pm \frac{\Omega }{2}. \end{aligned}$$

We conclude that

$$\begin{aligned} e^{xg_{\pm }(\lambda )\sigma _3}=e^{\pm \frac{x\Omega }{2}\sigma _3}(1+O(\lambda ))\quad \text{ as }\quad \lambda \rightarrow 0. \end{aligned}$$

In a similar way it can be proved that that $f_{\pm }(\lambda )=e^{\pm \frac{\Delta }{2}}(1+O(\lambda ))$ as $\lambda \rightarrow 0$. Using the above expansion we have that as $\lambda \rightarrow 0_+$

$$\begin{aligned} \begin{aligned}&Y_+(\lambda )e^{xg_+(\lambda )\sigma _3}f_+(\lambda )^{\sigma _3}(P_+^{\infty }(\lambda ))^{-1}=-\dfrac{Y_1(0)}{2\lambda } \begin{bmatrix} e^{\frac{x\Omega +\Delta }{2}}&e^{-\frac{x\Omega +\Delta }{2}} \end{bmatrix} \\&\quad \times \left( \begin{bmatrix} p_+(0)S^{\infty }_{2+}(0)-\nabla _{\Omega }S^{\infty }_{2+}(0) &{}p_+(0)S^{\infty }_{2+}(0)-\nabla _{\Omega }S^{\infty }_{2+}(0)\\ p_+(0)S^{\infty }_{1+}(0)+\nabla _{\Omega }S^{\infty }_{1+}(0)&{} p_+(0)S^{\infty }_{1+}(0)+\nabla _{\Omega }S^{\infty }_{1+}(0) \end{bmatrix} +O(\lambda )\right) \,. \end{aligned} \end{aligned}$$

(3.92)

Using the relations (3.67) and (3.68) we obtain

$$\begin{aligned} Y_+(\lambda )e^{xg_+(\lambda )\sigma _3}f_+(\lambda )^{\sigma _3}(P_+^{\infty }(\lambda ))^{-1}=O(1),\quad \text{ as }\quad \lambda \rightarrow 0_+ . \end{aligned}$$

In a similar way it can be verified the regular behaviour at $0_-$ which concludes the proof of the lemma.

Let $\Sigma _\mathcal {E}$ be the system of contours shown in Fig. 5. The above arguments show that the error vector $\mathcal {E}$ satisfies the following Riemann–Hilbert problem

$$\begin{aligned} {\mathcal E}_+(\lambda )={\mathcal E}_-(\lambda ) V_\mathcal {E}(\lambda ) \qquad \lambda \in \Sigma _\mathcal {E} \end{aligned}$$

and as $\lambda \rightarrow \infty$

$$\begin{aligned} \mathcal E(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O} \left( \lambda ^{-1} \right) \ . \end{aligned}$$

(3.93)

where the jump matrix $V_{\mathcal {E}}$ satisfies

$$\begin{aligned} V_\mathcal {E} (\lambda ) = {\left\{ \begin{array}{ll} I + \mathcal {O} \left( e^{-c{|x|}} \right) &{} \lambda \in \mathcal {C}_j, \ j=1,2, \\ I + \mathcal {O}\left( |x|^{-1} \right) &{} \lambda \in \partial B_\rho ^{\pm \eta _j}, \ j=1,2. \end{array}\right. } \end{aligned}$$

(3.94)

Therefore, by a standard small norm argument (see, for example [Its11, Section 5.1.3]) there exists a unique solution $\mathcal E$, which possesses an asymptotic expansion for large negative x and large $\lambda$ of the form:

$$\begin{aligned} \mathcal E (\lambda ) = \begin{bmatrix} 1&1 \end{bmatrix} \ + \ \frac{\mathcal {E}_{1}(x)}{x \lambda } \ + \mathcal {O}\left( \frac{1}{\lambda ^{2}} \right) \ , \end{aligned}$$

(3.95)

where $\mathcal {E}_{1}(x)$ possesses bounded derivatives in 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.

We also notice that the solution $\mathcal E$ which we have constructed obeys the symmetry

$$\begin{aligned} \mathcal E(-\lambda ) = \mathcal E(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \ . \end{aligned}$$

(3.96)

Indeed, the jump matrices $V_{\mathcal E}$ for $\mathcal E$ all satisfy the symmetry

$$\begin{aligned} V_{\mathcal E}(-\lambda ) = \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} V_{\mathcal E}(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}\ \end{aligned}$$

(3.97)

where $V_{\mathcal {E}}$ is given in (3.90) and (3.91). Properly, to see this, one must ensure that the contours for the Riemann–Hilbert problem for $\mathcal E$ are symmetric with respect to the mapping $\lambda \mapsto - \lambda$, and then verify that $V_{\mathcal E}$ satisfies (3.97). We have already specified in Sect. 3.5.2 that the circular contours should possess this symmetry, and it is clear that the lens boundaries may be chosen to satisfy this symmetry.

The verification of (3.97) for $\lambda$ in any of the four circles follows from the definitions (3.81) and (3.82). The verification of (3.97) for $\lambda$ in any of the lens boundaries follows by inspection of the jump matrices for S (only those defined on the lens boundaries) as described in Fig. 3, and using (3.69). The fact that (3.97) implies (3.96) is a straighforward exercise from the theory of Riemann–Hilbert problems.

Because $\mathcal E$ is analytic in a vicinity of $\lambda =0$, the symmetry relation (3.96) implies that $\mathcal {E}(\lambda )$ has the expansion

$$\begin{aligned} \mathcal {E}(\lambda ) = \hat{c}_{0} \begin{bmatrix} 1&1\end{bmatrix} + \lambda \hat{c}_{1} \begin{bmatrix} 1&-1\end{bmatrix} \ , \text{ for } \lambda \ \text{ near } \ 0 , \end{aligned}$$

(3.98)

for some constants $\hat{c}_{0}$ and $\hat{c}_{1}$.

Taking 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:

$$\begin{aligned} \begin{aligned} Y(\lambda ) =\,&T(\lambda ) e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3} = S(\lambda ) e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3} \\ =\,&\mathcal E(\lambda ) P(\lambda ) e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3} =\left( \begin{bmatrix} 1&1\end{bmatrix} + \frac{\mathcal {E}_{1}(x)}{x \lambda } \ + \mathcal {O}\left( \frac{1}{\lambda ^{2}} \right) \right) P(\lambda ) e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3} \ , \end{aligned} \end{aligned}$$

(3.99)

where $P(\lambda )$ is the global parametrix defined by (3.87).

In particular, for $\lambda$ near 0, $\mathcal {E}(\lambda ) P(\lambda )$ appearing in (3.99) is actually $\mathcal {E} (\lambda ) P^{\infty }(\lambda )$. The reader will recall that $P^{\infty }$ has a pole at $\lambda = 0$ (see (3.63)). However, because of the behavior of $\mathcal {E}(\lambda )$ for $\lambda$ near 0 shown in (3.98), the product $\mathcal {E}(\lambda ) P^{\infty }(\lambda )$ has no pole at $\lambda =0$.

We recall that the potential u(x) can be calculated from the solution $Y(\lambda )$ as

$$\begin{aligned} u(x) = 2\frac{\mathrm{d}}{\mathrm{d}x} \left[ \lim _{\lambda \rightarrow \infty } \lambda (Y_1(\lambda ;x) -1) \right] , \end{aligned}$$

(3.100)

where $Y_1(\lambda ; x)$ is the first entry of the vector Y.

Theorem 3.6

In the regime $x \rightarrow - \infty$, the potential u(x) has the following asymptotic behaviour

$$\begin{aligned} u(x) = \eta _2^2-\eta _1^2 -2\eta _2^2\dfrac{E(m)}{K(m)} -2\dfrac{\partial ^2 }{\partial x^2}\log \vartheta _3 \left(\frac{\eta _2}{2K(m)}(x+\phi );2\tau \right)+ \mathcal {O}\left(|x|^{-1}\right)\nonumber \\ \end{aligned}$$

(3.101)

where E(m) and K(m) are the complete elliptic integrals of the first and second kind respectively with modulus $m=\eta _1/\eta _2$, $\phi$ is given by

$$\begin{aligned} \phi =\int _{\eta _1}^{\eta _2}\dfrac{\log r(\zeta )}{R_+(\zeta )}\dfrac{\mathrm{d}\zeta }{\pi i}\in \mathbb {R}\, \end{aligned}$$

(3.102)

and $2\tau =i\frac{K(m')}{K(m)}$, $m'=\sqrt{1-m^2}$. The formula (3.101) can be written in the equivalent form

$$\begin{aligned} u(x)=\eta _2^2-\eta _1^2-2\eta _2^2{\text {dn}}^2\left( \eta _2(x+\phi ) +K(m)\left| \, m\right. \right)+ \mathcal {O}\left(|x|^{-1}\right) \end{aligned}$$

(3.103)

where ${\text {dn}}\left(z \left| \, m\right. \right)$ is the Jacobi elliptic function of modulus m.

Proof

We are interested in the first entry of the vector $Y(\lambda )$ (for $\lambda$ large), and we have, from (3.99),

$$\begin{aligned} Y(\lambda )&=\left( \begin{bmatrix} 1&1\end{bmatrix} +\displaystyle \frac{\mathcal {E}_{1}(x)}{x \lambda } \ + \mathcal {O}\left( \frac{1}{\lambda ^{2}} \right) \right) P^{\infty }(\lambda )e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3}\\&=\left( S^{\infty }(\lambda )+ \displaystyle \frac{\mathcal {E}_{1}(x)}{x \lambda } \ + \mathcal {O}\left( \frac{1}{\lambda ^{2}} \right) \right) e^{-xg(\lambda ) \sigma _3}f(\lambda )^{-\sigma _3}\,. \end{aligned}$$

Hence

$$\begin{aligned} Y_1(\lambda ) =\left[ S^{\infty }_{1}(\lambda ) + \frac{\left( \mathcal {E}_{1}\right) _{1}(x)}{x \lambda } \ + \mathcal {O}\left( \frac{1}{\lambda ^{2}} \right) \right] \frac{e^{-xg(\lambda ) }}{f(\lambda )} \, . \end{aligned}$$

(3.104)

From the expression (3.17) for the g function, we have

$$\begin{aligned} e^{-xg(\lambda ) } = 1 - \frac{x}{\lambda } \left[ \frac{\eta _1^2+\eta _2^2}{2} +\eta _2^2\left( \dfrac{E(m)}{K(m)}-1\right) \right] + \mathcal {O}\left(\frac{1}{\lambda ^2}\right) \, . \end{aligned}$$

(3.105)

From the formula of $f(\lambda )$ in (3.27) we have

$$\begin{aligned} f(\lambda )=1+\frac{f_1}{\lambda }+\mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, , \end{aligned}$$

where $f_1$ is independent of x. Starting from the vector $S^{\infty }(\lambda )$ in (3.50) we observe that $\gamma (\lambda ) = 1 + \mathcal {O} \left(\lambda ^{-2} \right)$, using (3.37) and (3.40) we have

$$\begin{aligned} 2w(\lambda ) = -\frac{1}{2} - \dfrac{1}{\lambda }\dfrac{\Omega }{2\pi i}+\mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, ,\;\quad \frac{\Omega }{2\pi i}=-\dfrac{\eta _2}{2K(m)} \end{aligned}$$

so expanding (3.50) gives

$$\begin{aligned} S^{\infty }_{1}(\lambda )&= 1 - \frac{1}{\lambda } \frac{\Omega }{2\pi i} \left[ \frac{ \vartheta _3' \left( \frac{x\Omega +\Delta }{2\pi i}; 2 \tau \right) }{ \vartheta _3 \left( \frac{x\Omega +\Delta }{2\pi i}; 2 \tau \right) } - \frac{ \vartheta _3'(0; 2\tau )}{\vartheta _3(0;2\tau )} \right] + \left(\frac{1}{\lambda ^2}\right) \\&= 1 - \dfrac{1}{\lambda }\dfrac{\partial }{\partial x}\log \vartheta _3 \left(\frac{x\Omega +\Delta }{2\pi i };2\tau \right)+ \mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, , \end{aligned}$$

where we have used the property that $\vartheta _3'(0; 2\tau )=0$ because $\vartheta _3(z; 2\tau )$ is an even function of z. Therefore

$$\begin{aligned} Y_1(\lambda )= & {} 1+\dfrac{1}{\lambda }\left( f_1-x \left[ \frac{\eta _1^2+\eta _2^2}{2} +\eta _2^2\left( \dfrac{E(m)}{K(m)}-1\right) \right] \right. \\&\left. -\dfrac{\partial }{\partial x}\log \vartheta _3 \left(\frac{x\Omega +\Delta }{2\pi i };2\tau \right)+ \frac{\left( \mathcal {E}_{1}(x) \right) _{1}}{x} \right) +\mathcal {O}\left(\frac{1}{\lambda ^2}\right) \,. \end{aligned}$$

From the above expansions, using (3.100), and the explicit expression of $\Delta$ in (3.28), we obtain the expression of u(x) in (3.101). In order to obtain the expression (3.103) we need the following identity (see e.g. [Law89] pg. 45 exercise 16 and 3.5.5)

$$\begin{aligned} \dfrac{1}{4K^2(m)}\dfrac{\mathrm{d}^2}{\mathrm{d}z^2} \log \vartheta _3(z;2\tau )=-\dfrac{E(m)}{K(m)}+{\text {dn}}^2\left( 2K(m) z+K(m)\left| \, m\right. \right)\, , \end{aligned}$$

(3.106)

where ${\text {dn}}\left( z\left| \, m\right. \right)$ is the Jacobi elliptic function of modulus m and period 2K(m) and we recall that $2\tau =iK(m')/K(m)$. Then we can write

$$\begin{aligned} \dfrac{\partial ^2 }{\partial x^2}\log \vartheta _3 \left(\frac{x\Omega +\Delta }{2\pi i };2\tau \right)=-\eta _2^2\dfrac{E(m)}{K(m)}+\eta _2^2{\text {dn}}^2\left( \eta _2(x+\phi ) + K(m)\left| \, m\right. \right)\, , \end{aligned}$$

so that the expression for u(x) in (3.101) can be written in the form (3.103).

4 Behaviour of the Potential u(x, t) as $t \rightarrow +\infty$

Letting the potential u(x, t) evolve in time according to the KdV equation, the reflection coefficient evolves as ${r}(\lambda ; t) = {r}(\lambda ) e^{-8\lambda ^3t}$. This will lead to the study of a Riemann–Hilbert problem Y for the soliton gas described as follows

$$\begin{aligned}&Y_+ (\lambda ) = Y_-(\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} 1 &{}0 \\ -i r(\lambda ) e^{8\lambda t \left( \lambda ^2 - \frac{x}{4t} \right)} &{} 1 \end{bmatrix} &{}\quad \lambda \in \Sigma _{1}\\ \\ \displaystyle \begin{bmatrix} 1 &{} ir(\lambda ) e^{-8\lambda t \left( \lambda ^2 - \frac{x}{4t} \right)} \\ 0 &{} 1 \end{bmatrix}&\quad \lambda \in \Sigma _{2} \end{array}\right. } \end{aligned}$$

(4.1)

$$\begin{aligned}&Y(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned}$$

(4.2)

We are interested in the asymptotic behaviour of $Y(\lambda )$ in the long-time regime ($t\rightarrow + \infty$).

The phase appearing in the exponents in the jump matrix shows different sign depending on the value of the quantity

$$\begin{aligned} \xi = \frac{x}{4t} \in \mathbb {R}\ . \end{aligned}$$

(4.3)

It is clear that in the case $\xi > \eta _2^2$, the phases in the jumps are exponentially decaying in the regime $t \rightarrow + \infty$, therefore by a straightforward small norm argument we conclude

$$\begin{aligned} Y(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left( {e^{-8\eta _1(\xi ^2-\eta _2^2)t}} \right) \qquad \text {as } t\rightarrow + \infty {\text { with } \xi ^2>\eta _2^2,} \end{aligned}$$

(4.4)

and the potential u(x, t) becomes trivial.

The more interesting case $\xi \le \eta _2^2$ will be studied below. It will become clear that we will observe the presence of a critical value $\xi _\mathrm{crit}$ at which a phase transition occurs when passing from $\xi > \xi _\mathrm{crit}$ (the “super-critical" case) to $\xi \le \xi _\mathrm{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.

5 Super-Critical Case: The $\alpha$-Dependency

We first consider the case

$$\begin{aligned} \xi _\mathrm{crit}< \xi < \eta _2^2 \end{aligned}$$

(5.1)

where the value of $\xi _\mathrm{crit} \in \mathbb {R}$ will be defined in (5.18).

In order to study the Riemann–Hilbert problem for Y in this setting we need to split the contours in the following way: let $\alpha \in (\eta _1, \eta _2)$ and define the sub intervals

$$\begin{aligned} \Sigma _{1,\alpha } = (\alpha , \eta _2) \subseteq \Sigma _1 \qquad \text {and} \qquad \Sigma _{2,\alpha }= (-\eta _2, -\alpha ) \subseteq \Sigma _2 \ . \end{aligned}$$

(5.2)

The value of $\alpha$ will be determined in Eq. (5.16) as a function of $\xi$.

We introduce again scalar functions $g(\lambda )$ and $f(\lambda )$ (in a slight abuse of notation, we are using the same letter g and f to denote these functions, though properly we should probably use $g_{\alpha }$ and $f_{\alpha }$). We make the first transformation $Y(\lambda )\mapsto T(\lambda )$ given by

$$\begin{aligned} T(\lambda ) = Y(\lambda ) e^{t g(\lambda )\sigma _3}f(\lambda )^{\sigma _3} \end{aligned}$$

(5.3)

such that

$$\begin{aligned}&g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda = 0&\lambda \in \Sigma _{1,\alpha }\cup \Sigma _{2,\alpha } \end{aligned}$$

(5.4)

$$\begin{aligned}&g_+(\lambda )-g_-(\lambda ) = \widetilde{\Omega }&\lambda \in [-\alpha , \alpha ] \end{aligned}$$

(5.5)

$$\begin{aligned}&g(\lambda ) = \mathcal {O}\left(\frac{1}{\lambda }\right)&\lambda \rightarrow \infty \ . \end{aligned}$$

(5.6)

We further require that $g(\lambda )-4\lambda ^3+4\xi \lambda - \widetilde{\Omega }$ behaves as $(\lambda {\mp } \alpha )^{\frac{3}{2}}$ near $\lambda =\pm \alpha$. In addition, there are two types of inequalities that must be satisfied by this function in order to have a successful Riemann–Hilbert analysis. First we will need inequalities satisfied on the complement (relative to $\Sigma _1\cup \Sigma _2$) of the sets $\Sigma _{1,\alpha }$ and $\Sigma _{2,\alpha }$:

$$\begin{aligned}&{\text {Re}}\left[ g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda \right] < 0 \ \ \ \ \lambda \in (\eta _{1},\alpha ) \end{aligned}$$

(5.7)

$$\begin{aligned}&{\text {Re}}\left[ g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda \right] > 0 \ \ \ \ \lambda \in (-\alpha , -\eta _{1}) \ . \end{aligned}$$

(5.8)

Second, we will require monotonicity properties on $\Sigma _{1}$ and $\Sigma _{2}$:

$$\begin{aligned}&-i( g_{+}(\lambda ) - g_{-}(\lambda )) \text{ is } \text{ purely } \text{ real } \text{ and } \text{ monotonically } \text{ decreasing } \text{ on } (\alpha , \eta _{2}) \ \end{aligned}$$

(5.9)

$$\begin{aligned}&-i( g_{+}(\lambda ) - g_{-}(\lambda )) \text{ is } \text{ purely } \text{ real } \text{ and } \text{ monotonically } \text{ increasing } \text{ on } (-\eta _{2}, - \alpha ) \ .\nonumber \\ \end{aligned}$$

(5.10)

It is well-known that there is a unique function g satisfying all these properties, which we will define explicitly here (we will actually define $g'$, which of course determines g). We define

$$\begin{aligned} g'(\lambda ) = - 12\lambda ^2 + 4\xi + 12 \dfrac{Q_2(\lambda )}{R_\alpha (\lambda ) }-4\xi \dfrac{Q_1(\lambda )}{R_\alpha (\lambda ) }\, , \end{aligned}$$

(5.11)

where

$$\begin{aligned} R_\alpha (\lambda ) =\sqrt{(\lambda ^2-\alpha ^2)(\lambda ^2-\eta _2^2)} \end{aligned}$$

(5.12)

is taken to be analytic in $\mathbb {C}\backslash \left\rbrace \Sigma _{1,\alpha }\cup \Sigma _{2,\alpha } \right\lbrace$ and real and positive on $(\eta _2,+\infty )$; moreover, let

$$\begin{aligned} Q_1(\lambda )=\lambda ^2+c_1\, ,\quad \text{ and } \ Q_2(\lambda )=\lambda ^4-\frac{1}{2}\lambda ^2(\alpha ^2+\eta _2^2)+c_2\, . \end{aligned}$$

(5.13)

The constants $c_1$ and $c_2$ are chosen so that

$$\begin{aligned} \int _{-\alpha }^{\alpha }\dfrac{Q_2(\zeta )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta = 0 \ , \quad \int _{-\alpha }^{\alpha }\dfrac{Q_1(\zeta )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta = 0 \ . \end{aligned}$$

(5.14)

Explicitly, we find

$$\begin{aligned} c_1=- \eta _2^2+\eta _2^2\dfrac{E(m_{\alpha })}{K(m_{\alpha })}\, ,\quad c_2=\dfrac{1}{3}\alpha ^2\eta _2^2+\dfrac{1}{6}(\eta _2^2+\alpha ^2)c_1 \ , \quad m_{\alpha }=\dfrac{\alpha }{\eta _2}\, , \end{aligned}$$

(5.15)

where $K(m_{\alpha })$ and $E(m_{\alpha })$ are, respectively, the complete elliptic integrals of the first and second kind.

The parameter $\alpha$ is determined by requiring that the function $g(\lambda )-4\lambda ^3 + 4\xi \lambda - \widetilde{\Omega }$ has a zero at $\lambda =\pm \alpha$, which yields the equation

$$\begin{aligned} \xi =3\dfrac{Q_2(\pm \alpha )}{Q_1(\pm \alpha )}=\dfrac{1}{2}(\alpha ^2+\eta _2^2)+\dfrac{\alpha ^2(\alpha ^2-\eta _2^2)}{\alpha ^2-\eta _2^2+\eta _2^2\frac{E(m_{\alpha })}{K(m_{\alpha })}}\, , \end{aligned}$$

(5.16)

and this determines the constant $\alpha$ implicitly as a function of $\xi$.

Before continuing our analysis we want to comment on Eq. (5.16). We can rewrite it in the form

$$\begin{aligned} \xi =\dfrac{x}{4t}=\dfrac{\eta _2^2}{2}W(m_\alpha )\, ,\quad W(m_{\alpha })=\left[ 1+m_{\alpha }^2+2\dfrac{m_{\alpha }^2(1-m_{\alpha }^2)}{1-m_{\alpha }^2-\frac{E(m_{\alpha })}{K(m_{\alpha })}}\right] \, . \end{aligned}$$

(5.17)

This relation describes the modulation of the parameter $\alpha$ as a function of $\xi$. The quantity $\eta _2^2W(m_{\alpha })$ 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 $\eta _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)=-\eta _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

$$\begin{aligned} \dfrac{E(m_\alpha )}{K(m_\alpha )}&=1-\dfrac{1}{2}m_\alpha ^2+\mathcal {O}(m_\alpha ^4)\, ,\ \text {as } m_\alpha \rightarrow 0 \quad \text {and}\quad \\&\dfrac{E(m_\alpha )}{K(m_\alpha )}\simeq \dfrac{2}{\log (8/(1-m_{\alpha }))}\, ,\ \text {as } m_\alpha \rightarrow 1\, , \end{aligned}$$

so that

$$\begin{aligned} \lim _{\alpha \rightarrow 0 }\dfrac{3Q_2(\alpha )}{Q_1(\alpha )}=-\frac{3\eta _2^2}{2}\, , \quad \text {and} \quad \lim _{\alpha \rightarrow \eta _2}\dfrac{3Q_2(\alpha )}{Q_1(\alpha )}=\eta _2^2 \, . \end{aligned}$$

The Whitham equations are strictly hyperbolic ([Lev88]), so that $\dfrac{\partial }{\partial \alpha }W(m_\alpha )>0$ for $0<\alpha <\eta _2$. Hence by the implicit function theorem, the equation (5.17) defines $\alpha$ as a monotone increasing function of $\xi$ for $\xi \in [\xi _\mathrm{crit},\eta _2^2]$ where $\xi _\mathrm{crit}$ is given by

$$\begin{aligned} \xi _\mathrm{crit}=\dfrac{3Q_2(\eta _1)}{Q_1(\eta _1)}=\dfrac{1}{2}(\eta _1^2+\eta _2^2)+\dfrac{\eta _1^2(\eta _1^2-\eta _2^2)}{\eta _1^2-\eta _2^2+\eta _2^2\frac{E(m)}{K(m)}}\, , \quad m = \frac{\eta _1}{\eta _2}\, . \end{aligned}$$

(5.18)

Then, clearly $\xi _\mathrm{crit}>-\frac{3\eta _2^2}{2}$.

From $g'(\lambda )$, we also have a representation of $g(\lambda )$:

$$\begin{aligned} g(\lambda )=-4\lambda ^3+4\xi \lambda +12\int _{\eta _2}^\lambda \dfrac{Q_2(\zeta )}{R_\alpha (\zeta ) }\mathrm{d}\zeta -4\xi \int _{\eta _2}^\lambda \dfrac{Q_1(\zeta )}{R_\alpha (\zeta ) }\mathrm{d}\zeta \, . \end{aligned}$$

(5.19)

This, together with (5.5), yields the formula

$$\begin{aligned} \widetilde{\Omega }=24\int _{\eta _2}^{\alpha }\dfrac{Q_2(\zeta )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta -8\xi \int _{\eta _2}^{\alpha }\dfrac{Q_1(\zeta )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta \, . \end{aligned}$$

(5.20)

For future use we will need the x derivatives of $t g(\lambda )$ and $t\widetilde{\Omega }$. Before calculating them, let us observe that

$$\begin{aligned} \widetilde{\Omega }=24\int _{\eta _2}^{\alpha }\dfrac{Q_2(\zeta )-Q_2(\alpha )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta -8\xi \int _{\eta _2}^{\alpha }\dfrac{Q_1(\zeta )-Q_1(\alpha )}{R_{\alpha +}(\zeta )}\mathrm{d}\zeta \, , \end{aligned}$$

which gives, using the Riemann bilinear relations (see e.g. [Spr57]),

$$\begin{aligned} \widetilde{\Omega }= 2\pi i \frac{4\xi - 2(\alpha ^2+\eta _2^2 )}{\int _{-\alpha }^{\alpha } \frac{\mathrm{d}\zeta }{R_\alpha (\zeta ) }}=2\pi i \eta _2 \frac{\alpha ^2+\eta _2^2 -2\xi }{K(m_{\alpha })} \in i\mathbb {R}\ ,\quad m_{\alpha }=\dfrac{\alpha }{\eta _2}\, . \end{aligned}$$

(5.21)

Fig. 6

Construction of the genus-1 Riemann surface $\mathfrak {X}_\alpha$ and its basis of cycles

Lemma 5.1

The following identities are satisfied

$$\begin{aligned}&\dfrac{\partial }{\partial x}tg'(\lambda )=1-\dfrac{Q_1(\lambda )}{R_\alpha (\lambda )}\, , \end{aligned}$$

(5.22)

$$\begin{aligned}&\dfrac{\partial }{\partial x}t\widetilde{\Omega }=- \frac{ \pi i \eta _2 }{K(m_{\alpha })}\, . \end{aligned}$$

(5.23)

Proof

We observe that $g'(\lambda )\mathrm{d}\lambda$ defined in (5.11) is a meromorphic one-form on the Riemann surface $\mathfrak {X}_{\alpha }$ defined as

$$\begin{aligned} \mathfrak {X}_{\alpha }=\left\rbrace (\eta ,\lambda )\in \mathbb {C}^2\;|\; \eta ^2=R_\alpha ^2(\lambda )=(\lambda ^2-\alpha ^2)(\lambda ^2-\eta _2^2)\right\lbrace \, . \end{aligned}$$

We define a homology basis on $\mathfrak {X}_{\alpha }$ in the following way: the B cycle encircles the cut $[\alpha ,\eta _2]$ clockwise and the A cycle starts on the cut $[-\eta _2,-\alpha ]$ on the upper semi-plane, goes to the cut $[\alpha ,\eta _2]$ and then goes back to $[-\eta _2,-\alpha ]$ on the second sheet of $\mathfrak {X}_{\alpha }$. See Fig. 6. Then we have

$$\begin{aligned} \oint _{A}g'(\zeta )\, \mathrm{d}\zeta =0\, ,\quad \oint _{B}g'(\zeta )\, \mathrm{d}\zeta =-\widetilde{\Omega }\, . \end{aligned}$$

(5.24)

Regarding the first relation in (5.22) we have

$$\begin{aligned} \dfrac{\partial }{\partial x}tg'(\lambda )\mathrm{d}\lambda&=\dfrac{\partial }{\partial x}\left[ - 12t\lambda ^2 \mathrm{d}\lambda + x\mathrm{d}\lambda + 12t \dfrac{Q_2(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda -x \dfrac{Q_1(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda \right] \\&=\mathrm{d}\lambda -\dfrac{Q_1(\lambda )}{R_\alpha (\lambda )}\mathrm{d}\lambda +\dfrac{\partial }{\partial \alpha }\left[ 12t \dfrac{Q_2(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda -x \dfrac{Q_1(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda \right] \dfrac{\partial \alpha }{\partial x}\\&=\mathrm{d}\lambda -\dfrac{Q_1(\lambda )}{R_\alpha (\lambda )}\mathrm{d}\lambda \, , \end{aligned}$$

(5.25)

because the term $\dfrac{\partial }{\partial \alpha }\left[ 12t \dfrac{Q_2(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda -x \dfrac{Q_1(\lambda )}{R_\alpha (\lambda ) }\mathrm{d}\lambda \right]$ vanishes since it is a holomorphic one-form (no singularity at $\pm \alpha$ or infinity) which is normalized to zero on the A cycle because of (5.24); therefore it is identically zero [Kri88] (see also [Gra02, GT02]). An alternative proof is to calculate the derivative and use the explicit formulæ of the constants $c_1$ and $c_2$ in (5.15). We conclude that

$$\begin{aligned} \dfrac{\partial }{\partial x}e^{-tg(\lambda ) }=- \frac{1}{\lambda } \left[ \frac{\alpha ^2+\eta _2^2}{2} +\eta _2^2\left( \dfrac{E(m_\alpha )}{K(m_\alpha )}-1\right) \right] + \mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, . \end{aligned}$$

Regarding the relation (5.23), by (5.22) and (5.24) we have

$$\begin{aligned} \dfrac{\partial }{\partial x}(t \widetilde{\Omega })=-\dfrac{\partial }{\partial x} \oint _{B}tg'(\lambda )\, \mathrm{d}\lambda =- \oint _{B}\dfrac{\partial }{\partial x}(tg'(\lambda )\, \mathrm{d}\lambda )=- \frac{ \pi i \eta _2 }{K(m_{\alpha })} \, . \end{aligned}$$

As we did in Sect. 3, we choose the function f to simplify the jumps on $\Sigma _{1,\alpha }$ and $\Sigma _{2,\alpha }$ via

$$\begin{aligned}&f_+(\lambda ) f_-(\lambda ) =\frac{1}{r(\lambda )}&\lambda \in \Sigma _{1,\alpha } \end{aligned}$$

(5.26)

$$\begin{aligned}&f_+(\lambda ) f_-(\lambda ) = r(\lambda )&\lambda \in \Sigma _{2,\alpha } \end{aligned}$$

(5.27)

$$\begin{aligned}&\frac{f_+(\lambda )}{ f_-(\lambda )} =e^{\widetilde{\Delta }}&\lambda \in [-\alpha ,\alpha ] \end{aligned}$$

(5.28)

$$\begin{aligned}&f(\lambda ) = 1+\mathcal {O}\left(\frac{1}{\lambda }\right)&\lambda \rightarrow \infty \ . \end{aligned}$$

(5.29)

It is easy to check that the function $f(\lambda )$ is given by

$$\begin{aligned} f(\lambda )= & {} {\text {exp}}\left\{ \frac{R_{\alpha }(\lambda )}{2\pi i}\left[ \int _{ \Sigma _{1,\alpha }} \frac{\log \frac{1}{r(\zeta )} }{R_{\alpha +}(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta + \int _{\Sigma _{2,\alpha }} \frac{\log r(\zeta )}{R_{\alpha +}(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta \right. \right. \nonumber \\&\left. \left. +\int ^{\alpha }_{ -\alpha } \frac{\widetilde{\Delta }}{R_{\alpha }(\zeta )(\zeta -\lambda )} \mathrm{d}\zeta \right] \right\} \, , \end{aligned}$$

(5.30)

where the constraint (5.31) determines $\widetilde{\Delta }$ as

$$\begin{aligned} \widetilde{\Delta }&=\left[ \int _{\Sigma _{1, \alpha }} \frac{\log r(\zeta )}{R_{\alpha +}(\zeta )} \mathrm{d}\zeta - \int _{\Sigma _{2, \alpha }} \frac{\log r(\zeta )}{R_{\alpha +}(\zeta )} \mathrm{d}\zeta \right] \left[ \int ^{\alpha }_{ -\alpha } \frac{ \mathrm{d}\zeta }{R_{\alpha }(\zeta )}\right] ^{-1}\nonumber \\&=2\left[ \int _{\Sigma _{1, \alpha }} \frac{\log r(\zeta )}{R_{\alpha +}(\zeta )} \mathrm{d}\zeta \right] \left[ \int ^{\alpha }_{ -\alpha } \frac{ \mathrm{d}\zeta }{R_{\alpha }(\zeta )}\right] ^{-1}\, , \end{aligned}$$

(5.31)

where in the last relation in (5.33) we have used the fact that $r(-\lambda )=r(\lambda )$.

As a consequence, T satisfies the following Riemann–Hilbert problem:

$$\begin{aligned}&T_+(\lambda ) = T_-(\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} e^{t\left(g_+(\lambda ) - g_-(\lambda )\right)} \frac{f_{+}(\lambda )}{f_{-}(\lambda )}&{} 0 \\ -i&{} e^{-t\left(g_+(\lambda ) - g_-(\lambda )\right)} \frac{f_{-}(\lambda )}{f_{+}(\lambda )}\end{bmatrix} &{}\quad \lambda \in \Sigma _{1,\alpha }\\ \displaystyle \begin{bmatrix} e^{t\left(g_+(\lambda ) - g_-(\lambda )\right)} \frac{f_{+}(\lambda )}{f_{-}(\lambda )}&{}i \\ 0 &{} e^{-t\left(g_+(\lambda ) - g_-(\lambda )\right)} \frac{f_{-}(\lambda )}{f_{+}(\lambda )} \end{bmatrix} &{} \quad \lambda \in \Sigma _{2,\alpha } \\ \displaystyle \begin{bmatrix} e^{\widetilde{\Omega }t + \widetilde{ \Delta } } &{} 0 \\ - i r(\lambda ) f_{+}(\lambda )f_{-}(\lambda ) e^{t \left(g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda \right)}&{} e^{-\widetilde{\Omega }t - \widetilde{ \Delta } } \end{bmatrix} &{}\quad \lambda \in [\eta _1,\alpha ]\\ \displaystyle \begin{bmatrix} e^{\widetilde{\Omega }t +\widetilde{ \Delta } } &{} e^{-t \left(g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda \right)} \frac{ i r(\lambda ) }{f_{+}(\lambda )f_{-}(\lambda )} \\ 0 &{} e^{-\widetilde{\Omega }t- \widetilde{ \Delta } } \end{bmatrix} &{} \quad \lambda \in [-\alpha , -\eta _1] \\ \displaystyle \begin{bmatrix} e^{\widetilde{\Omega }t+\widetilde{ \Delta } } &{}0 \\ 0 &{} e^{-\widetilde{\Omega }t-\widetilde{ \Delta } } \end{bmatrix}&\quad \lambda \in [-\eta _1, \eta _1] \end{array}\right. } \end{aligned}$$

(5.32)

$$\begin{aligned}&T(\lambda ) = \begin{bmatrix}1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \, . \end{aligned}$$

(5.33)

5.1 Opening lenses

It is useful to provide representations of the entries appearing in the jump matrix for $T(\lambda )$ in either $\Sigma _{1,\alpha }$ or $\Sigma _{2,\alpha }$, that clearly demonstrate their analytic continuation off these intervals, as was done in Sect. 3. The following formulæ are valid on both intervals:

$$\begin{aligned}&g_{+}(\lambda ) - g_{-}(\lambda ) = 2 g_{+}(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \ , \end{aligned}$$

(5.34)

$$\begin{aligned}&g_{+}(\lambda ) - g_{-}(\lambda ) = -\left( 2 g_{-}(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right) \ . \end{aligned}$$

(5.35)

The following formulæ are valid on $\Sigma _{1,\alpha }$:

$$\begin{aligned}&\frac{f_{+}(\lambda )}{f_{-}(\lambda )} =- \frac{1}{ f_{-}^{2}(\lambda ) \hat{r}_{-}(\lambda ) } \quad \text {and} \quad \frac{f_{-}(\lambda )}{f_{+}(\lambda )} = \frac{1}{ f_{+}^{2}(\lambda ) \hat{r}_{+}(\lambda )} \ . \end{aligned}$$

(5.36)

And the following ones are valid on $\Sigma _{2,\alpha }$:

$$\begin{aligned}&\frac{f_{-}(\lambda )}{f_{+}(\lambda )} = - \frac{f_{-}^{2}(\lambda )}{ \hat{r}_{-}(\lambda )} \quad \text {and} \quad \frac{f_{+}(\lambda )}{f_{-}(\lambda )} = \frac{f_{+}^{2}(\lambda )}{ \hat{r}_{+}(\lambda )} \ . \end{aligned}$$

(5.37)

As it was done in Sect. 3, we can factor the jump matrix on $\Sigma _{1,\alpha }$ as follows

$$\begin{aligned}&\begin{bmatrix} e^{t\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )} &{} 0 \\ -i &{} e^{-t\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix}\\&\quad \quad =\begin{bmatrix}1&{}-\dfrac{ie^{ -t \left( 2 g_{-}(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right) )}}{f_-^2 (\lambda )\hat{r}_-(\lambda )}\\ 0 &{}1 \end{bmatrix} \begin{bmatrix}0&{}-i \\ -i &{}0\end{bmatrix} \begin{bmatrix}1&{}\dfrac{ie^{ -t \left( 2 g_{+}(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right) }}{ \hat{r}_+(\lambda )f_+^2(\lambda )}&{} \\ 0&{}1 \end{bmatrix} \end{aligned}$$

and on $\Sigma _{2,\alpha }$ as

$$\begin{aligned}&\begin{bmatrix} e^{t\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_+(\lambda )}{f_-(\lambda )}&{}i \\ 0 &{} e^{-t\left(g_+(\lambda ) - g_-(\lambda )\right)}\dfrac{f_-(\lambda )}{f_+(\lambda )} \end{bmatrix}\\&\quad \quad =\begin{bmatrix}1&{}0\\ i\dfrac{f_-^2 (\lambda )}{\hat{r}_-(\lambda )}e^{ t \left( 2 g_{-} (\lambda )+ 8 \lambda ^{3} - 8 \xi \lambda \right) )} &{}1 \end{bmatrix} \begin{bmatrix}0&{}i\\ i&{}0\end{bmatrix} \begin{bmatrix}1&{}0\\ -i\dfrac{f_+^2(\lambda )}{ \hat{r}_+(\lambda )}e^{ t\left( 2 g_{+}(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right) }&{} &{}1 \end{bmatrix}\, . \end{aligned}$$

These factorizations motivate us to open lenses $\mathcal {C}_j$ around each $\Sigma _{j,\alpha }$, $j=1,2$. We use these lenses to define the transformation

(5.38)

The lens contours and the resulting jump relations for S(z) are shown in Fig. 7.

Fig. 7

Opening lenses: the gray entries in the jumps represent exponentially small quantities in the limit $t \rightarrow + \infty$. The contours $\mathcal {C}_{1}$ and $\mathcal {C}_{2}$ are the lens boundaries

Lemma 5.2

The following inequalities are satisfied:

$$\begin{aligned}&{\text {Re}}\left[2 g(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right] >0 \ \text{ for } \lambda \in \mathcal {C}_{1} \backslash \{ \alpha , \eta _{2} \} \ , \end{aligned}$$

(5.39)

$$\begin{aligned}&{\text {Re}}\left[2 g(\lambda ) + 8 \lambda ^{3} - 8 \xi \lambda \right] <0 \ \text{ for } \lambda \in \mathcal {C}_{2} \backslash \{ - \eta _{2}, - \alpha \} \ , \end{aligned}$$

(5.40)

$$\begin{aligned}&{\text {Re}}\left[ g_+(\lambda ) + g_-(\lambda ) + 8\lambda ^3 - 8\xi \lambda \right] < 0 \text{ for } \lambda \in [\eta _{1}, \alpha ) \ , \end{aligned}$$

(5.41)

$$\begin{aligned}&{\text {Re}}\left[ g_+(\lambda ) + g_-(\lambda ) + 8\lambda ^3 - 8\xi \lambda \right] > 0 \text{ for } \lambda \in (-\alpha , - \eta _{1}] \ . \end{aligned}$$

(5.42)

Proof

Using (5.16) the function $g'(\lambda )$ in (5.11) can be written in the form

$$\begin{aligned} g'(\lambda )=- 12\lambda ^2 + 4\xi + 12 \dfrac{Q_2(\lambda )-Q_2(\alpha )}{R_\alpha (\lambda ) }-4\xi \dfrac{Q_1(\lambda )-Q_1(\alpha )}{R_\alpha (\lambda ) }\, , \end{aligned}$$

so that we have

$$\begin{aligned} g'_+(\lambda ) - g'_-(\lambda ) = - i24 \frac{\sqrt{\lambda ^2-\alpha ^2}}{ \sqrt{\eta _2^2-\lambda ^2}} \left[ \lambda ^2 - \left(\frac{\eta _2^2 -\alpha ^2}{2} + \frac{\xi }{3} \right) \right] \end{aligned}$$

(5.43)

and from (5.14) we deduce that the quadratic polynomial has one root $\rho _+$ in the interval $[0, \alpha ]$ which is positive for $\lambda > \alpha$. Therefore, for $\lambda \in \Sigma _{1,\alpha }$

$$\begin{aligned} {\text {Im}}\left[g'_+(\lambda ) - g'_-(\lambda )\right] =- 24 \frac{\sqrt{\lambda ^2-\alpha ^2}}{ \sqrt{\eta _2^2-\lambda ^2}} \left[ \lambda ^2 - \left(\frac{\eta _2^2 -\alpha ^2}{2} + \frac{\xi }{3} \right) \right]<0 \ . \end{aligned}$$

(5.44)

From the formula (5.19) for g we also have that for $\lambda \in [\eta _1,\alpha ]$

$$\begin{aligned} g_+(\lambda ) + g_-(\lambda ) + 8\lambda ^3 - 8\xi \lambda = - 24\int _\lambda ^\alpha \frac{\sqrt{\alpha ^2-\zeta ^2}}{\sqrt{\eta _2^2-\zeta ^2}} \left[ \zeta ^2- \left(\frac{\eta _2^2 -\alpha ^2}{2} + \frac{\xi }{3} \right) \right] \mathrm{d}\zeta \ .\nonumber \\ \end{aligned}$$

(5.45)

Setting

$$\begin{aligned} h_{\alpha , \xi }(\zeta ) = \frac{\sqrt{\alpha ^2-\zeta ^2}}{\sqrt{\eta _2^2-\zeta ^2}} \left[ \zeta ^2 - \left(\frac{\eta _2^2 -\alpha ^2}{2} + \frac{\xi }{3} \right) \right] \ ,\end{aligned}$$

(5.46)

we need to show that the function

$$\begin{aligned} H_{\alpha , \xi }(\lambda ) =\int _{\lambda }^\alpha - h_{\alpha ,\xi }(\zeta ) <0 \qquad \text {for } \lambda \in [\eta _1,\alpha ] \ . \end{aligned}$$

(5.47)

It is easy to check that $H_{\alpha ,\xi } (\alpha ) =0$ and $H_{\alpha ,\xi } (0) = 0$ (see (5.14)). Next, $H_{\alpha ,\xi }'(\lambda ) = h_{\alpha ,\xi }(\lambda )$ is negative on $[0,\rho _+]$ and positive on $[\rho _+,\alpha ]$. This implies that indeed the inequality (5.49) is satisfied on $[\eta _1,\alpha ]$.

Because of Lemma 5.2, letting $t\rightarrow +\infty$, the jump matrices (as depicted in Fig. 7) will converge to constant jumps exponentially fast outside neighbourhoods of $\pm \alpha$ and $\pm \eta _2$. We then obtain the following model Riemann–Hilbert problem for $\widetilde{S}^{\infty }$:

$$\begin{aligned} \widetilde{S}^{\infty }_+(\lambda ) = \widetilde{S}^{\infty }_-(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} e^{t\widetilde{\Omega }+\widetilde{ \Delta } } &{}0 \\ 0 &{} e^{-t\widetilde{\Omega }-\widetilde{ \Delta } }\end{bmatrix} &{} \lambda \in [-\alpha ,\alpha ] \\ \begin{bmatrix}0 &{} -i \\ -i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{1,\alpha } \\ \begin{bmatrix}0 &{} i \\ i &{} 0 \end{bmatrix}&\lambda \in \Sigma _{2,\alpha } \end{array}\right. } \end{aligned}$$

(5.48)

$$\begin{aligned} \widetilde{ S}^{\infty }(\lambda )=\begin{bmatrix} 1&1\end{bmatrix}+\mathcal {O}\left(\frac{1}{\lambda }\right),\quad \lambda \rightarrow \infty \, . \end{aligned}$$

(5.49)

5.2 The outer parametrix $\widetilde{P}^{\infty }$

Along the same lines as we did in Sect. 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.

We therefore seek a matrix valued function $\widetilde{P}^{\infty }$ that is analytic in $\mathbb {C}\backslash (-\eta _2,\eta _2)$ and satisfies the following Riemann–Hilbert problem

$$\begin{aligned} \widetilde{P}^{\infty }_+(\lambda ) =\widetilde{ P}^{\infty }_-(\lambda ) {\left\{ \begin{array}{ll} \begin{bmatrix} e^{t\widetilde{\Omega }+\widetilde{ \Delta } } &{}0 \\ 0 &{} e^{-t\widetilde{\Omega }-\widetilde{ \Delta } }\end{bmatrix} &{} \lambda \in [-\alpha ,\alpha ] \\ \begin{bmatrix}0 &{} -i \\ -i &{} 0 \end{bmatrix} &{}\lambda \in \Sigma _{1,\alpha } \\ \begin{bmatrix}0 &{} i \\ i &{} 0 \end{bmatrix}&\lambda \in \Sigma _{2,\alpha } \end{array}\right. } \end{aligned}$$

(5.50)

$$\begin{aligned} \widetilde{P}^{\infty }(\lambda )=\begin{bmatrix}1&{}0\\ 0&{}1\end{bmatrix}+\mathcal {O}\left(\frac{1}{\lambda }\right),\quad \lambda \rightarrow \infty \, . \end{aligned}$$

(5.51)

In order to get the solution of the above Riemann-Hilbert problem, let us introduce in analogy to Sect. 3.4 the vector

$$\begin{aligned} \begin{aligned} \widetilde{S}^{\infty }(\lambda )=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3\left( \frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i};2\tau \right)} \begin{bmatrix} \displaystyle \frac{\vartheta _3 \left(2\widetilde{w}(\lambda ) +\frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3 \left(2\widetilde{w}(\lambda ) -\frac{1}{2};2\tau \right)}&\displaystyle \frac{\vartheta _3 \left(-2\widetilde{w}(\lambda ) +\frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3 \left(-2\widetilde{w}(\lambda ) -\frac{1}{2};2\tau \right)} \end{bmatrix} \end{aligned} \ , \end{aligned}$$

(5.52)

with $\widetilde{w}(\lambda )$ defined as

$$\begin{aligned} \widetilde{w}(\lambda ) = \int _{\eta _{2}}^{\lambda } \frac{ \Omega _\alpha }{R_{\alpha }(\lambda )} \frac{d \lambda }{4 \pi i} \, \end{aligned}$$

(5.53)

where $\Omega _\alpha =- \frac{ \pi i \eta _2 }{K(m_{\alpha })}$. Further from (5.21) we have

$$\begin{aligned} \widetilde{\Omega }=2\pi i \eta _2 \frac{\alpha ^2+\eta _2^2}{K(m_{\alpha })}+4\xi \Omega _\alpha \end{aligned}$$

and

$$\begin{aligned} p_{\alpha }(\lambda )=\int _{\eta _2}^\lambda \dfrac{Q_1(\zeta )}{R_{\alpha }(\zeta )}\mathrm{d}\zeta ,\quad \Omega _{\alpha }=-2p_{\alpha +}(\alpha ) \end{aligned}$$

where $Q_1$ has been defined in (5.13). We note that for $\lambda \in (-\alpha , \alpha )$, we have

$$\begin{aligned} p_{\alpha +}(\lambda ) - p_{\alpha _-}(\lambda ) = - \Omega _{\alpha } \ . \end{aligned}$$

(5.54)

Then the solution $\widetilde{P}^{\infty }(\lambda )$ to the Riemann-Hilbert problem (5.52) and (5.53) is given explicitly by

$$\begin{aligned} \widetilde{P}^{\infty }(\lambda )=\frac{1}{2}\begin{bmatrix} (1+\frac{p_{\alpha }(\lambda )}{\lambda })\widetilde{S}^{\infty }_1(\lambda )+\dfrac{1}{\lambda }\nabla _{\Omega _{\alpha }}\widetilde{S}^{\infty }_1(\lambda ) &{}(1-\frac{p_{\alpha }(\lambda )}{\lambda })\widetilde{S}^{\infty }_2(\lambda )+\dfrac{1}{\lambda }\nabla _{\Omega _{\alpha }}\widetilde{S}^{\infty }_2(\lambda )\\ (1-\frac{p_{\alpha }(\lambda )}{\lambda })\widetilde{S}^{\infty }_1(\lambda )-\dfrac{1}{\lambda }\nabla _{\Omega _{\alpha }}\widetilde{S}^{\infty }_1(\lambda )&{} (1+\frac{p_{\alpha }(\lambda )}{\lambda })\widetilde{S}^{\infty }_2(\lambda )-\dfrac{1}{\lambda }\nabla _{\Omega _\alpha }\widetilde{S}^{\infty }_2(\lambda ) \end{bmatrix}\,, \end{aligned}$$

(5.55)

where $\widetilde{S}_{1}^{\infty }$ and $\widetilde{S}_{2}^{\infty }$ are the entries of the row vector $\widetilde{S}^\infty$ defined in (5.54), and

$$\begin{aligned}&\nabla _{\Omega _{\alpha }}\widetilde{S}^{\infty }_1(\lambda ):=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3 \left(2\widetilde{w}(\lambda ) -\frac{1}{2};2\tau \right)} \dfrac{\Omega _{\alpha }}{2\pi i}\dfrac{\mathrm{d}}{\mathrm{d}z}\left[ \frac{\vartheta _3 \left(z+2\widetilde{w}(\lambda ) +\frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3\left(z+ \frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i};2\tau \right)}\right] \bigg |_{z=0} \ ,\\&\nabla _{\Omega _{\alpha }}\widetilde{S}^{\infty }_2(\lambda ):=\gamma (\lambda )\dfrac{\vartheta _3(0;2\tau )}{\vartheta _3 \left(-2\widetilde{w}(\lambda ) -\frac{1}{2};2\tau \right)} \dfrac{\Omega _{\alpha }}{2\pi i}\dfrac{\mathrm{d}}{\mathrm{d}z}\left[ \frac{\vartheta _3 \left(z-2\widetilde{w}(\lambda ) +\frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i}-\frac{1}{2};2\tau \right)}{\vartheta _3\left(z+ \frac{t \widetilde{\Omega }+\widetilde{\Delta }}{2\pi i};2\tau \right)}\right] \bigg |_{z=0} \ . \end{aligned}$$

The above construction has been obtained by modifying the construction of $P^{\infty }$ in (3.63), in such a way that $\widetilde{P}^{\infty }(\lambda )$ solves the Riemann–Hilbert problem (5.52)–(5.53), with $\text{ det }\widetilde{P}^{\infty }(\lambda ) = 1$, and $\widetilde{P}^{\infty }(-\lambda ) = \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \widetilde{P}^{\infty }(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix}$.

5.3 The local parametrix $P^{\pm \alpha }$

We will construct now a (matrix) local parametrix around the points $\lambda = \pm \alpha$. The construction of the local parametrices near $\lambda =\pm \eta _2$ is the same one as in Sect. 3.5.

We focus again on a small but fixed neighbourhood $B^{(-\alpha )}_{\rho } = \left\rbrace \lambda \in \mathbb {C} \left| \, \left|\lambda + \alpha \right|< \rho \right. \right\lbrace$ of the endpoint $\lambda = -\alpha$. We define the conformal map

$$\begin{aligned} \zeta= & {} \left(\frac{3}{4}\right)^{\frac{2}{3}} \left[ t\int _{-\alpha }^\lambda g'_+(s) - g'_-(s) \mathrm{d}s \right]^{\frac{2}{3}}\nonumber \\= & {} \left[ 18 t\int _{-\alpha }^\lambda \left(\frac{\sqrt{\alpha ^2-s^2}}{\sqrt{\eta _2^2-s^2}}\right)_+ \left(s^2 - \frac{\eta _2^2-\alpha ^2}{2} - \frac{\xi }{3}\right) \mathrm{d}s \right]^{\frac{2}{3}} \end{aligned}$$

(5.56)

locally in $B^{(-\alpha )}_\rho$.

To define the local parametrix $P^{-\alpha }$ in $B^{(-\alpha )}_\rho$, we consider

$$\begin{aligned} P(\lambda ) = S (\lambda )e^{ \frac{\pi i }{4} \sigma _{3}}\left( \frac{\sqrt{ \pm \hat{r}(\lambda )}}{f(\lambda )} \right) ^{\sigma _{3}} e^{ {\mp } \frac{1}{2} \left( \widetilde{\Omega } t + \widetilde{\Delta } \right) \sigma _{3}}\ , \ \lambda \in B^{(-\alpha )}_{\rho } \cap \mathbb {C}_{\pm }, \end{aligned}$$

and then, using the inverse of the transformation $\zeta (\lambda )$, we define

$$\begin{aligned} P^{(1)}(\zeta ) = P(\lambda (\zeta )) e^{-\frac{2}{3} \zeta ^{\frac{3}{2}} \sigma _3} \ , \quad \zeta \in \mathbb {C}\end{aligned}$$

with branch cut $(-\infty ,0]$. By construction, $P^{(1)}$ satisfies a Riemann–Hilbert problem with jumps in a neighbourhood of $\zeta = 0$ as shown in Fig. 8.

Fig. 8

The contour setting under the conformal map $\zeta$ in a neighbourhood of 0

We introduce the (local) Airy parametrix (see [Dei99, DKM+99]): let $\Psi _\mathrm {Ai}(\zeta )$ be the solution to the following Riemann–Hilbert problem

• (a) $\Psi _\mathrm {Ai}$ is analytic for $\zeta \in \mathbb {C} \backslash \Gamma _{\Psi }$, where the contours $\Gamma _\Psi$ are defined as $\Gamma _\pm = \left\rbrace \arg \zeta = \pm \frac{2\pi }{3}\right\lbrace$, $\Gamma _{0, -} =\left\rbrace \arg \zeta = \pi \right\lbrace$ and $\Gamma _{0,+} =\left\rbrace \arg \zeta = 0 \right\lbrace$;

• (b) $\Psi$ satisfies the following jump relations

  $$\begin{aligned} \Psi _{\mathrm {Ai}\, +} (\zeta ) = \Psi _{\mathrm {Ai}\, -}(\zeta ) {\left\{ \begin{array}{ll} \begin{bmatrix} 1 &{} 0 \\ 1 &{} 1 \end{bmatrix} &{} \text {on } \Gamma _+ \text { and } \Gamma _-\\ \begin{bmatrix} 0 &{} 1 \\ -1\ &{} 0 \end{bmatrix} &{} \text {on }\Gamma _{0,-} \\ \begin{bmatrix} 1 &{} 1 \\ 0\ &{} 1 \end{bmatrix}&\text {on }\Gamma _{0,+}\, ; \end{array}\right. }\, \end{aligned}$$

  (5.57)

• (c) as $\zeta \rightarrow \infty$

  $$\begin{aligned} \Psi _\mathrm {Ai}(\zeta ) = \zeta ^{-\frac{1}{4}\sigma _3} \frac{1}{\sqrt{2}}\begin{bmatrix} 1&{} i \\ i&{}1 \end{bmatrix} \left(I + \mathcal {O}\left(\frac{1}{\zeta ^{\frac{3}{2}}}\right)\right) e^{-\frac{2}{3}\zeta ^{\frac{3}{2}}\sigma _3} \ , \end{aligned}$$

  (5.58)

• (d) $\Psi _\mathrm {Ai}$ remains bounded as $\zeta \rightarrow 0$, $\zeta \in \mathbb {C} \backslash \Gamma _{\Psi }$.

The solution to this Riemann–Hilbert problem is constructed with the help of Airy functions. Setting $\omega = e^{\frac{2\pi i}{3}}$, we have

$$\begin{aligned}&\Psi _\mathrm {Ai}(\zeta ) = \sqrt{2\pi }\begin{bmatrix} \displaystyle \mathrm {Ai}(\zeta )&{}\displaystyle -\omega ^2 \mathrm {Ai}(\omega ^2\zeta ) \\ \displaystyle -i\mathrm {Ai}'(\zeta ) &{} \displaystyle i\omega \mathrm {Ai}'(\omega ^2\zeta ) \end{bmatrix}&\text {for } 0< \arg \zeta < \frac{2\pi }{3} \end{aligned}$$

(5.59)

$$\begin{aligned}&\Psi _\mathrm {Ai}(\zeta ) = \sqrt{2\pi } \begin{bmatrix} \displaystyle -\omega \mathrm {Ai}(\omega \zeta )&{}\displaystyle -\omega ^2 \mathrm {Ai}(\omega ^2\zeta ) \\ \displaystyle i\omega ^2 \mathrm {Ai}'(\zeta ) &{} \displaystyle i\omega \mathrm {Ai}'(\omega ^2\zeta ) \end{bmatrix}&\text {for } \frac{2\pi }{3}< \arg \zeta < \pi \end{aligned}$$

(5.60)

$$\begin{aligned}&\Psi _\mathrm {Ai}(\zeta ) = \sqrt{2\pi } \begin{bmatrix} \displaystyle -\omega ^2 \mathrm {Ai}(\omega ^2\zeta ) &{} \displaystyle \omega \mathrm {Ai}(\omega \zeta )\\ \displaystyle i\omega \mathrm {Ai}'(\omega ^2\zeta ) &{} \displaystyle -i\omega ^2 \mathrm {Ai}'(\zeta ) \end{bmatrix}&\text {for } -\pi< \arg \zeta < - \frac{2\pi }{3} \end{aligned}$$

(5.61)

$$\begin{aligned}&\Psi _\mathrm {Ai}(\zeta ) = \sqrt{2\pi } \begin{bmatrix} \displaystyle \mathrm {Ai}(\zeta ) &{} \displaystyle \omega \mathrm {Ai}(\omega \zeta )\\ \displaystyle - i \mathrm {Ai}'(\zeta ) &{} \displaystyle -i\omega ^2 \mathrm {Ai}'(\zeta ) \end{bmatrix}&\text {for } - \frac{2\pi }{3}< \arg \zeta < 0 \ , \end{aligned}$$

(5.62)

where $\mathrm {Ai}(\zeta )$ is the Airy function.

In conclusion, our local parametrix is then defined as

$$\begin{aligned} P^{-\alpha }(\zeta (\lambda )) = A(\lambda ) \Psi _{\mathrm {Ai}}(\zeta (\lambda )) e^{\frac{2}{3}\zeta ^{\frac{3}{2}} \sigma _3} e^{\pm \frac{1}{2} \left( \widetilde{\Omega } t + \widetilde{\Delta } \right) } \left( \frac{f(\lambda )}{\sqrt{ \pm \hat{r}(\lambda ) } } \right) ^{\sigma _{3}}e^{ - \frac{\pi i }{4} \sigma _{3}} \ , \ \lambda \in B^{(-\alpha )}_{\rho } \cap \mathbb {C}_{\pm } \, ,\nonumber \\ \end{aligned}$$

(5.63)

where A is an analytic prefactor whose expression is determined by imposing that

$$\begin{aligned} P^{-\alpha } (\lambda ) \left(\widetilde{P}^{\infty }(\lambda ) \right)^{-1} = I + \mathcal {O}\left(t^{-1}\right) \qquad \text {as } t \rightarrow + \infty \, , \ \text {for } \lambda \in \partial B^{(-\alpha )}_\rho \backslash \Gamma _\Psi \ . \end{aligned}$$

(5.64)

In light of this asymptotic behaviour we set

$$\begin{aligned} A(\lambda ) = \widetilde{P}^{\infty }(\lambda ) e^{ {\mp } \frac{1}{2} \left( \widetilde{\Omega } t + \widetilde{\Delta } \right) \sigma _{3}} e^{ \frac{\pi i}{4} \sigma _{3}} \left( \frac{ \sqrt{ \pm \hat{r}(\lambda ) } }{f(\lambda )} \right) ^{\sigma _{3} } \frac{1}{\sqrt{2}} \begin{bmatrix}1&{} -i \\ - i &{} 1\end{bmatrix} \zeta (\lambda )^{\frac{1}{4}\sigma _3} \ , \ \text{ for } \lambda \in B^{(-\alpha )}_{\rho } \cap \mathbb {C}_{\pm } \ .\nonumber \\ \end{aligned}$$

(5.65)

By construction, A is defined and analytic in a neighbourhood of $-\alpha$, minus the cuts $(-\infty , -\alpha ] \cup [-\alpha ,+ \infty )$; moreover, A is invertible ($\det A(\lambda ) \equiv 1$).

Lemma 5.3

$A(\lambda )$ is analytic in an open neighbourhood of $-\alpha$.

Proof

The proof entails verifying that A has no jumps across the interval $(-\alpha -\rho , - \alpha +\rho )$ and that it has at most a removable singularity at $\lambda = -\alpha$. We leave the verification that $A_{+}(\lambda ) = A_{-}(\lambda )$ across the interval $(-\alpha -\rho ,- \alpha +\rho )$ to the reader, using the jump relations satisfied by $\widetilde{P}^{\infty }$ and the above definitions.

The conformal map $\zeta (\lambda )$ has a simple zero at $\lambda = -\alpha$ (by construction), therefore $\zeta (\lambda )^{-\frac{1}{4}\sigma _3}$ has at most a fourth-root singularity at $-\alpha$. Similarly, $\widetilde{P}^{\infty }(\lambda )$ has a fourth-root singularity at $- \alpha$, as well; therefore, all the entries of $A(\lambda )$ have at most a square-root singularity at $\lambda = -\alpha$, and $A(\lambda )$ is analytic in $B^{(-\alpha )}_\rho \backslash \{ -\alpha \}$. The point $\lambda = - \alpha$ is a removable singularity. This implies that $A(\lambda )$ is indeed analytic everywhere in $B^{(-\alpha )}_\rho \$.

The construction of the local parametrix for $\lambda$ near $\alpha$ is obtained from the parametrix near $-\alpha$ as follows. We fix a disk of the same radius as the radius of the disk used for the parametrix near $-\alpha$, and within that disk, define

$$\begin{aligned}&P^{\alpha } := \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} P^{-\alpha }(-\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \ . \end{aligned}$$

(5.66)

Remark. We note that, as with the analysis presented for $t=0$ and $x \rightarrow - \infty$, we may choose the contours so that they are preserved under the transformation $\lambda \mapsto - \lambda$. Moreover, the construction of $\widetilde{P}^{\infty }$ continues to enjoy the symmetry relation

$$\begin{aligned} \widetilde{P}^{\infty }(-\lambda ) = \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \widetilde{P}^{\infty }(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \ . \end{aligned}$$

(5.67)

5.4 Small norm argument and determination of u(x, t) as $t\rightarrow + \infty$

As before, we define a global (matrix) parametrix P replacing each model $P^\xi$ in (3.87) with $\widetilde{P}^\xi$, $\xi = \infty ,\pm \alpha , \pm \eta _2$. Then we define the following “remainder” Riemann–Hilbert problem:

$$\begin{aligned} \mathcal E(\lambda ) = S(\lambda ) P(\lambda )^{-1}\ . \end{aligned}$$

(5.68)

For some $c>0$, the vector $\mathcal E$ satisfies

$$\begin{aligned} \mathcal E_+(\lambda ) = {\left\{ \begin{array}{ll} \mathcal E_-(\lambda ) \left( I + \mathcal {O}\left( \displaystyle { e^{-ct} }\right) \right) \qquad \text {on the upper and lower lenses, outside the discs } \\ \mathcal E_-(\lambda ) \left( I + \mathcal {O}\left(\displaystyle t^{-1} \right) \right) \qquad \text {on the circles around the endpoints} \\ \end{array}\right. } \end{aligned}$$

(5.69)

and

$$\begin{aligned} \mathcal E(\lambda ) = \begin{bmatrix} 1&1 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \text {as }\lambda \rightarrow \infty \ . \end{aligned}$$

(5.70)

Furthermore, as in Sect. 3.6, the solution $\mathcal {E}(\lambda )$ is analytic in a neighborhood of $\lambda =0$. Moreover, the jumps $V_{\mathcal {E}}$ for $\mathcal {E}$ satisfy the symmetry

$$\begin{aligned} V_{\mathcal {E}}(-\lambda ) = \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} V_{\mathcal {E}}(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} \ . \end{aligned}$$

(5.71)

Therefore, by a small norm argument (see [Its11, Section 5.1.3]), we learn that there is a unique $\mathcal {E}$ solving the Riemann–Hilbert problem, and (as in Sect. 3.6), the solution satisfies the symmetry relation

$$\begin{aligned} \mathcal {E}(-\lambda ) = \mathcal {E}(\lambda ) \begin{bmatrix}0&{}1\\ 1&{}0\end{bmatrix} , \end{aligned}$$

(5.72)

and has a complete asymptotic expansion, satisfying

$$\begin{aligned} \mathcal E(\lambda ) = \begin{bmatrix} 1&1 \end{bmatrix} + \frac{\mathcal {E}_{1}( x, t)}{ \lambda t}\ + \mathcal {O} \left( \frac{1}{\lambda ^{2}} \right) , \end{aligned}$$

(5.73)

where $\mathcal {E}_{1}(x,t)$ and its derivatives are bounded.

Unraveling the transformations, we can again get back to the potential. Our original Riemann–Hilbert problem, for the unknown Y, satisfies

$$\begin{aligned} Y(\lambda )= & {} T(\lambda ) e^{-tg(\lambda ) \sigma _3} f(\lambda )^{-\sigma _3} = S(\lambda ) e^{-tg(\lambda ) \sigma _3} f(\lambda )^{-\sigma _3}\\= & {} \left( \begin{bmatrix} 1&1\end{bmatrix} + \frac{\mathcal {E}_{1}( x, t)}{ \lambda t}\ + \mathcal {O} \left( \frac{1}{\lambda ^{2}} \right) \right)P(\lambda ) e^{-tg(\lambda ) \sigma _3} f(\lambda )^{-\sigma _3}. \end{aligned}$$

In particular we are interested in the vector $Y(\lambda )$ for large $\lambda$

$$\begin{aligned} Y(\lambda )= & {} \left( \begin{bmatrix} 1&1\end{bmatrix} + \frac{\mathcal {E}_{1}( x, t)}{ \lambda t}\ + \mathcal {O} \left( \frac{1}{\lambda ^{2}} \right) \right)\widetilde{P}^{\infty }(\lambda ) e^{-tg(\lambda ) \sigma _3} f(\lambda )^{-\sigma _3}\\= & {} \left( \widetilde{S}^{\infty }(\lambda ) + \frac{\mathcal {E}_{1}( x, t)}{ \lambda t}\ + \mathcal {O} \left( \frac{1}{\lambda ^{2}} \right) \right)e^{-tg(\lambda ) \sigma _3} f(\lambda )^{-\sigma _3}, \end{aligned}$$

so that

$$\begin{aligned} Y_1(\lambda ) = \left[ \widetilde{S}^{\infty }_{1}(\lambda )+ \frac{\mathcal {E}_{1}( x, t)}{ \lambda t}\ + \mathcal {O} \left( \frac{1}{\lambda ^{2}} \right) \right]e^{-tg(\lambda ) }f(\lambda )^{-1} , \end{aligned}$$

(5.74)

where $\widetilde{S}^{\infty }_1$ refers to the the first row vector solution $\widetilde{S}^{\infty }$ (5.54). Since

$$\begin{aligned} u(x,t) = 2\frac{\mathrm{d}}{\mathrm{d}x} \left[ \lim _{\lambda \rightarrow \infty } \lambda (Y_1(\lambda ;x,t) -1) \right] \ , \end{aligned}$$

(5.75)

we have the following theorem.

Theorem 5.4

Given $\xi = \frac{x}{4t}$, in the region $\xi _\mathrm{crit}<\xi <\eta _2^2$ the solution of the KdV equation in the large time limit is

$$\begin{aligned} u(x,t)= \eta _2^2-\alpha ^2 -2\eta _2^2\dfrac{E(m_\alpha )}{K(m_\alpha )} -2\dfrac{\partial ^2}{\partial x^2} \log \vartheta _3 \left( \frac{\eta _2}{2K(m_\alpha )}(x-2(\alpha ^2+\eta _2^2)t+\widetilde{\phi });2\tau _\alpha \right) +\mathcal {O}(t^{-1}) \end{aligned}$$

(5.76)

where $E(m_\alpha )$ and $K(m_\alpha )$ are the complete elliptic integrals of first and second kind respectively, with modulus $m_\alpha = \frac{\alpha }{\eta _2}$; $2\tau _\alpha =i\dfrac{K(m'_\alpha )}{K(m_\alpha )}$, with $m'_\alpha = \sqrt{1-m_\alpha ^2}$,

$$\begin{aligned} \widetilde{\phi }=\int _{\alpha }^{\eta _2}\dfrac{\log r(\zeta )}{R_{\alpha +}(\zeta )}\dfrac{\mathrm{d}\zeta }{\pi i}\in \mathbb {R}\end{aligned}$$

and the parameter $\alpha =\alpha (\xi )$ is determined from the equation

$$\begin{aligned} \xi =\dfrac{\eta _2^2}{2}\left[ 1+m_{\alpha }^2+2\dfrac{m_{\alpha }^2(1-m_{\alpha }^2)}{1-m_{\alpha }^2-\frac{E(m_{\alpha })}{K(m_{\alpha })}}\right] . \end{aligned}$$

The error term $\mathcal {O}(t^{-1})$ is uniform for t sufficiently large.

Alternatively,

$$\begin{aligned} u(x,t)=\eta _2^2-\alpha ^2-2\eta _2^2{\text {dn}}^2\left( \eta _2(x-2(\alpha ^2+\eta _2^2)t+\widetilde{\phi }) + K(m_\alpha ) \left| \, m_\alpha \right. \right)+ \mathcal {O}\left(t^{-1}\right) \end{aligned}$$

(5.77)

where ${\text {dn}}\left( z\left| \, m \right.\right)$ is the Jacobi elliptic function.

Proof

Starting from (5.76) we expand each term of $Y_1(\lambda )$ in a neighbourhood of infinity. Regarding $f(\lambda )$ defined in (5.32) we have

$$\begin{aligned} f(\lambda )=1+\dfrac{f_1(\alpha ,\eta _2)}{\lambda }+\mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, , \end{aligned}$$

where

$$\begin{aligned} f_1(\alpha ,\eta _2)=\left[ \int _{\alpha }^{\eta _2}\dfrac{\zeta ^2 \log r(\zeta )}{R_\alpha (\zeta )}\dfrac{\mathrm{d}\zeta }{\pi i }-\widetilde{\Delta }\int _{-\alpha }^{\alpha }\dfrac{\zeta ^2}{R_\alpha (\zeta )}\dfrac{\mathrm{d}\zeta }{2\pi i }\right] \, . \end{aligned}$$

Regarding $e^{-tg(\lambda ) }$ we are interested in the x derivative of this expression. Using (5.22) we have

$$\begin{aligned} \dfrac{\partial }{\partial x}e^{-tg(\lambda ) }=- \frac{1}{\lambda } \left[ \frac{\alpha ^2+\eta _2^2}{2} +\eta _2^2\left( \dfrac{E(m_\alpha )}{K(m_\alpha )}-1\right) \right] + \mathcal {O}\left(\frac{1}{\lambda ^2}\right) \, . \end{aligned}$$

Regarding $\widetilde{S}^{\infty }_{1}(\lambda )$, we have

$$\begin{aligned} \widetilde{S}^{\infty }_{1}(\lambda )=1+ \dfrac{1}{\lambda }\left[ \left( \log \vartheta _3 \left(\frac{t\widetilde{\Omega }+\widetilde{\Delta }}{2\pi i };2\tau \right)\right) ' - \frac{\vartheta _{3}'(0)}{\vartheta _{3}(0)}\right] \dfrac{\eta _2}{2K(m_\alpha )}+ \mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, , \end{aligned}$$

where $'$ stands for the derivative with respect to the argument of the $\vartheta$-function. By (5.23) we have

$$\begin{aligned} \dfrac{\partial }{\partial x}\widetilde{S}^{\infty }_{1}(\lambda )= \dfrac{1}{\lambda }\left[ \log \vartheta _3 \left(\frac{t\widetilde{\Omega }+\widetilde{\Delta }}{2\pi i };2\tau \right) \right] ''\dfrac{\eta _2}{2K(m_\alpha )}\left( -\dfrac{\eta _2}{2K(m_\alpha )}+\dfrac{\partial }{\partial x}\dfrac{\widetilde{\Delta }}{2\pi i }\right) + \mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, , \end{aligned}$$

where $''$ stands for second derivative with respect to the argument of the $\vartheta$-function. Taking into account that $\Delta =\Delta ( \alpha (\xi ),\eta _2)$ and therefore the quantity $\frac{\partial }{\partial x} \Delta (\alpha (\xi ),\eta _2)=\mathcal {O}(t^{-1})$ by (5.16), we can write the above expression in the form

$$\begin{aligned} \dfrac{\partial }{\partial x}\widetilde{S}^{\infty }_{1}(\lambda )= -\dfrac{1}{\lambda }\left[ \dfrac{\partial }{\partial x^2}\log \vartheta _3 \left(\frac{t\widetilde{\Omega }+\widetilde{\Delta }}{2\pi i };2\tau \right)+\mathcal {O}(\frac{1}{t})\right] + \mathcal {O}\left(\frac{1}{\lambda ^2}\right)\, . \end{aligned}$$

Gathering the above expansions, using the fact that $\mathcal {E}_{1}(x,t)$ and its derivatives are bounded, and using the explicit expression of $\widetilde{\Omega }$ and $\widetilde{\Delta }$ in (5.21) and (5.33) respectively, we obtain (5.78). Also in this case, using the same calculations as in Theorem 3.6 we can reduce the expression of u(x, t) to the form (5.79).

The equivalence of formulas (5.78) and (5.79) is slightly more delicate than in the case of Theorem 3.6. It follows from (3.106) and the relation (5.23), that is a particular case of the more general relations obtained in [Kri88] in the context of Whitham modulation theory.

The equivalence of the formulas (5.78) and (5.79) is a well known fact in the theory of dispersive shock waves for the KdV equation where modulated travelling waves are developed (see e.g. the review [Gra12]). Such equivalence is valid for any solution of Whitham’s modulation equations. For the particular solution of the Whitham modulation equations appearing in the long time asymptotic analysis of KdV this equivalence was obtained in [G.16].

6 Sub-Critical Case

As the parameter $\xi < \eta _2^2$ decreases, we proved that there is a critical value $\xi _\mathrm{crit}$ (see Sect. 5, equation (5.18)) such that

$$\begin{aligned} \alpha (\xi _\mathrm{crit}) = \eta _1 \ .\end{aligned}$$

(5.78)

For $\xi < \xi _\mathrm{crit}$, we define

$$\begin{aligned} g'(\lambda ) = - 12\lambda ^2 + 4\xi + 12 \dfrac{Q_2(\lambda )}{R(\lambda ) }-4\xi \dfrac{Q_1(\lambda )}{R(\lambda ) } \end{aligned}$$

(5.79)

where R is defined in (3.16), specifically $R(\lambda ) =\sqrt{(\lambda ^2 - \eta _{1}^2) (\lambda ^2-\eta ^2_2)}$, and

$$\begin{aligned} Q_1(\lambda )=\lambda ^2+c_1\, ,\quad Q_2(\lambda )=\lambda ^4-\frac{1}{2}\lambda ^2(\eta _{1}^2+\eta _2^2)+c_2\, , \end{aligned}$$

(6.1)

with the constants $c_1$ and $c_2$ chosen so that

$$\begin{aligned} \int _{0}^{\eta _{1}}\dfrac{Q_2(\zeta )}{R_{+}(\zeta )}\mathrm{d}\zeta = 0 \ , \quad \int _{0}^{\eta _{1}}\dfrac{Q_1(\zeta )}{R_{+}(\zeta )}\mathrm{d}\zeta = 0 \ . \end{aligned}$$

(6.2)

Integration yields

$$\begin{aligned} g(\lambda ) = - 4 \lambda ^{3} + 4\xi \lambda + \int _{\eta _1}^{\lambda } \frac{12 Q_{2}(\zeta ) - 4 \xi Q_{1}(\zeta )}{R(\zeta )} \mathrm{d}\zeta \ . \end{aligned}$$

(6.3)

By construction, g satisfies the following constraints:

$$\begin{aligned}&g_+(\lambda ) + g_-(\lambda ) +8\lambda ^3 - 8\xi \lambda = 0&\lambda \in \Sigma _{1}\cup \Sigma _{2} \end{aligned}$$

(6.4)

$$\begin{aligned}&g_+(\lambda )-g_-(\lambda ) = \overline{ \Omega }&\lambda \in [-\eta _1,\eta _1] \end{aligned}$$

(6.5)

$$\begin{aligned}&g(\lambda ) = \mathcal {O}\left(\frac{1}{\lambda }\right)&\lambda \rightarrow \infty \ . \end{aligned}$$

(6.6)

with

$$\begin{aligned} \overline{\Omega }= 2\pi i \eta _2 \frac{2\xi - (\eta _1^2 +\eta _2^2) }{K(m)} \in i\mathbb {R}\ . \end{aligned}$$

(6.7)

Remark 6.1

The reader may verify that for $\xi = \xi _\mathrm{crit}$ the above function $g(\lambda ;\eta _1,\eta _2)$ in (6.2) agrees with the function $g(\lambda ;\alpha =\eta _1,\eta _2)$ in (5.19).

In order to show that the usual contour deformations can be carried out, as they were in Sects. 3 and 5, we need to verify that the quantity ${\text {Re}}\left[2 g(\lambda ) + 8 \lambda ^{3} - 8 \xi ^{2} \lambda \right]$ is positive on the contour $\mathcal {C}_{1}$, and negative on the contour $\mathcal {C}_{2}$, where these contours are as shown in Fig. 3.

To accomplish this, we consider the quadratic polynomial

$$\begin{aligned} q(r; \xi ) = 12\left( r^2-\frac{1}{2}r(\eta _{1}^2+\eta _2^2)+c_2\right) - 4 \xi (r + c_{1}) \ , \end{aligned}$$

(6.8)

with $r\in [0,\eta _1^2]$. A quick inspection shows that $q(\eta _1^2; \xi _\mathrm{crit}) = 0$ and $q(0; \xi _\mathrm{crit}) >0$, and moreover, for all $\xi \in \mathbb {R}$

$$\begin{aligned} \frac{\partial q}{\partial \xi }(0; \xi ) > 0 \quad \text {and} \quad \frac{\partial q}{\partial \xi }(\eta _1^2; \xi ) <0\ ; \end{aligned}$$

(6.9)

therefore, $0= q(\eta _1^2; \xi _\mathrm{crit}) < q(\eta _1^2; \xi )$ for all $\xi <\xi _\mathrm{crit}$. So, for all $\xi < \xi _\mathrm{crit}$, there are two roots of $q(r; \xi )$ within $(0, \eta _{1}^{2})$, and the polynomial is strictly positive on $[\eta _{1}^{2}, \eta _{2}^{2}]$.

This in turn implies, using arguments nearly identical to those used to prove Lemma 5.2, that

$$\begin{aligned}&{\text {Re}}\left[2 g(\lambda ) + 8 \lambda ^{3} - 8 \xi ^{2} \lambda \right] > 0 \ \text{ for } \lambda \in \mathcal {C}_{1}\backslash \{\eta _1,\eta _2\} \, , \end{aligned}$$

(6.10)

$$\begin{aligned}&{\text {Re}}\left[2 g(\lambda ) + 8 \lambda ^{3} - 8 \xi ^{2} \lambda \right] < 0 \ \text{ for } \lambda \in \mathcal {C}_{2} \backslash \{-\eta _1,-\eta _2\}\ . \end{aligned}$$

(6.11)

The use of this function, and the sequence of steps in the Riemann–Hilbert analysis which have been carried out for $t=0$ in Sect. 3, may be applied directly to the present situation, and we use the same outer model $P^\infty (\lambda )$ (cf. (3.63)) as was used in Sect. 3, with $x\Omega$ replaced by $t \overline{\Omega }$, with $\overline{\Omega }$ as defined by (6.9), along with the same local parametrices near each of the endpoints $\pm \eta _{1}, \pm \eta _{2}$. Therefore we arrive at the following result.

Theorem 6.2

In the regime $t \rightarrow + \infty$, $\xi < \xi _\mathrm{crit}$, the potential u(x, t) has the following asymptotic expansion

$$\begin{aligned} u(x,t)=\eta _2^2-\eta _1^2-2\eta _2^2{\text {dn}}^2 \left( \eta _2(x-2(\eta _1^2+\eta _2^2)t+\phi ) + K(m)\left| \, m\right. \right)+ \mathcal {O}\left(t^{-1}\right) \ ,\nonumber \\ \end{aligned}$$

(6.12)

where $m = \eta _{1}/\eta _{2}$, and

$$\begin{aligned} \phi = \int _{\eta _{1}}^{\eta _{2}} \frac{ \log {r( \zeta )}}{R_{+}(\zeta )} \frac{\mathrm{d}\zeta }{\pi i} \ . \end{aligned}$$

(6.13)

7 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 $+\infty$. 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 \rightarrow - \infty$, 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 behaviours depending on the value of the parameter $\xi = x/(4t)$: for $\xi > \eta ^{2}_2$ the solution decays exponentially, while for $\xi < \xi _{crit}$ the solution is described by the periodic travelling wave solution of KdV with fixed parameters; between, for $\xi \in (\xi _{\text {crit}}, \eta _2^2)$ these two extreme asymptotic states are connected by a periodic travelling wave solution of KdV with slowly varying parameters.

Several challenges remain, like the asymptotic analysis when there are two nontrivial reflection coefficients or 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.

Existence of Solution to the Soliton gas Riemann–Hilbert Problem

We will provide a proof of existence and uniqueness for the Riemann–Hilbert problem

$$\begin{aligned}&Y(\lambda ) \text { is analytic for } \lambda \in \mathbb {C}\backslash \left\rbrace \Sigma _1 \cup \Sigma _2 \right\lbrace \nonumber \\&Y_+ (\lambda ) = Y_-(\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} 1 &{}0 \\ -i r(\lambda ;x,t) &{} 1 \end{bmatrix} &{}\quad \lambda \in \Sigma _1\\ \displaystyle \begin{bmatrix} 1 &{} i r(\lambda ;x,t)\\ 0 &{} 1 \end{bmatrix}&\quad \lambda \in \Sigma _2 \end{array}\right. } \end{aligned}$$

(6.14)

$$\begin{aligned}&Y(\lambda ) = \begin{bmatrix}a&b \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \end{aligned}$$

(6.15)

$$\begin{aligned}&Y(-\lambda ) = Y(\lambda )\begin{bmatrix}0&{}\frac{b}{a}\\ \frac{a}{b}&{}0 \end{bmatrix} \,, \end{aligned}$$

(A.1)

with parameters $a>0$ and $b>0$, and with the contours shown in Fig. 2, and $r(\lambda ;x,t)=r(\lambda ) e^{8\lambda t \left( \lambda ^2 - \frac{x}{4t} \right)}$, and where we as usual seek a solution with at worst logarithmic singularities at the endpoints $\pm \eta _{j}$.

To establish that there is a solution to the Riemann–Hilbert problem (A.1)–(A.3), we seek $y_{1}$ with the following representation:

$$\begin{aligned} y_{1} = a + \frac{1}{2 \pi i } \int _{\eta _{1}}^{\eta _{2}} \frac{ \sqrt{r(s;x,t)} f(s)}{s - \lambda } \mathrm{d}s . \end{aligned}$$

(A.2)

This is consistent with the jump relations in (A.1), in which the first entry is analytic across $(-\eta _{2}, -\eta _{1})$, with a jump across $(\eta _{1}, \eta _{2})$. Plugging (A.4) into the jump relation across $(\eta _{1}, \eta _{2})$, we find

$$\begin{aligned} f(\lambda ) + \frac{b}{a} \frac{\sqrt{r( \lambda ;x,t)}}{2\pi } \int _{\eta _{1}}^{\eta _{2}} \frac{\sqrt{r(s;x,t)}f(s)}{s + \lambda } \mathrm{d}s = - i b \sqrt{r(\lambda ;x,t)} \ . \end{aligned}$$

(A.3)

The reader may verify that this integral equation appears (after some manipulation) in both entries of the jump relationships.

Now the integral operator appearing on the left hand side of (A.5) is compact (since it can obviously be approximated by a sequence of finite dimensional operators) and hence the index is zero. It is also positive definite, which shows that this integral equation is uniquely invertible.

To see that the integral operator is positive definite, we follow the classic [KM56, Formula 2.9], starting with the simple identity

$$\begin{aligned} \frac{1}{s+\lambda } =\int _{-\infty }^{0} e^{ ( s + \lambda ) z } \mathrm{d}z \ . \end{aligned}$$

(A.4)

We have

$$\begin{aligned}&\int _{\eta _{1}}^{\eta _{2}} \sqrt{r(\lambda ;x,t)} \ \overline{f(\lambda )} \int _{\eta _{1}}^{\eta _{2}} \frac{\sqrt{r(s;x,t)}f(s)}{s + \lambda } \mathrm{d}s \mathrm{d}\lambda \nonumber \\&\quad = \int _{-\infty }^{0} \int _{\eta _{1}}^{\eta _{2}} \int _{\eta _{1}}^{\eta _{2}} \sqrt{r(\lambda ;x,t)} \ \overline{f(\lambda )}\sqrt{r(s;x,t)}f(s) e^{(s + \lambda ) z} \mathrm{d}s \mathrm{d}\lambda \mathrm{d}z \nonumber \\&\quad = \int _{-\infty }^{0} \left| \int _{\eta _{1}}^{\eta _{2}} \sqrt{r(s;x,t)} f(s) e^{ s z } \mathrm{d}s \right| ^{2} \mathrm{d}z \ > \ 0 \ , \end{aligned}$$

(A.5)

provided f is not identically equal to 0.

Regarding uniqueness, if $\tilde{Y}$ is a solution to the Riemann–Hilbert problem (A.1)–(A.3), then setting $\tilde{f}(s) = - i \sqrt{r(s;x,t)} \tilde{Y}_{2}(s)$ for $s \in (\eta _{1}, \eta _{2})$, one verifies that $\tilde{f}$ must satisfy the integral equation (A.5), which obviously possesses a unique solution.

Returning to the Riemann–Hilbert problem (A.1)–(A.3), if we take $(a,b)=(1,1)$ we have established the existence and uniqueness of the solution to the soliton gas Riemann–Hilbert problem. But more importantly, if we separately consider $(a,b) = (1,2)$, we find a second independent solution of the Riemann–Hilbert problem, which combined yields a $2 \times 2$ matrix solution to the following Riemann–Hilbert problem:

$$\begin{aligned}&\mathbf{Y}(\lambda ) \text { is analytic for } \lambda \in \mathbb {C}\backslash \left\rbrace \Sigma _1 \cup \Sigma _2 \right\lbrace \nonumber \\&\mathbf{Y}_+ (\lambda ) = \mathbf{Y}_-(\lambda ) {\left\{ \begin{array}{ll} \displaystyle \begin{bmatrix} 1 &{}0 \\ -i r(\lambda ;x,t) &{} 1 \end{bmatrix} &{}\quad \lambda \in \Sigma _1\\ \displaystyle \begin{bmatrix} 1 &{} i r(\lambda ;x,t)\\ 0 &{} 1 \end{bmatrix}&\quad \lambda \in \Sigma _2 \end{array}\right. } \end{aligned}$$

(A.6)

$$\begin{aligned}&\mathbf{Y}(\lambda ) = \begin{bmatrix}1&{}1 \\ 1 &{} 2 \end{bmatrix} + \mathcal {O}\left(\frac{1}{\lambda }\right) \qquad \lambda \rightarrow \infty \ . \end{aligned}$$

(A.7)

This solution is invertible for all $\lambda \in \mathbb {C}$ since $\text{ det } \mathbf{Y} \equiv 1$.