1. Introduction
A soliton gas can be thought of as a very large number of solitons, with possibly random amplitudes and phases, interacting weakly (rarefied gas) or strongly (dense gas) in a medium modeled by an integrable nonlinear wave model. This concept goes back to the work of Zakharov [Reference Zakharov17]. There has been a growing interest in the study of such solutions at the analytical, numerical and experimental fronts in the past decade or so, and the literature on soliton gases is vast. We refer the reader to the comprehensive review article [Reference Suret, Randoux, Gelash, Agafontsev, Doyon and El12] on the subject and the references therein. The notion behind the construction of such solutions from an integrable systems point of view is the so-called ‘dressing’ method [Reference Zakharov and Manakov16]. This idea was employed in [Reference Dyachenko, Zakharov and Zakharov7] to construct the so-called primitive potentials, which can be taken as initial data for a solution modeling a soliton gas.
The current state of the art for computing soliton gas solutions of integrable nonlinear wave models is largely limited to either computing a suitable nonlinear superposition of N solitons for N large or numerically implementing an available asymptotic formula in a suitable asymptotic region of the space-time domain. The former approach is equivalent to computing an iterated (or N-fold) Darboux transformation, and it requires high-precision arithmetic because the size of the underlying linear algebra system is proportional to N (the number of solitons), and that linear system becomes numerically ill-conditioned as N grows. The latter approach has its difficulties as the accurate evaluation of asymptotic formulæ involving Riemann theta functions is a challenging task, if at all possible, especially when the genus of the underlying Riemann surface is not small. Moreover, the availability of spacetime asymptotic formulæ is limited to soliton gases obtained by the accumulation of eigenvalues on one pair of bands in the spectral plane. Regardless, such formulæ cease to be accurate for intermediate values of (x, t) outside of asymptotic regimes. An approach that is more closely related to our work was taken in [Reference Dyachenko, Zakharov and Zakharov7], where soliton gas (primitive) potentials were computed via the solution of a system of singular integral equations arising from a dressing construction. The numerical solution of this system necessitated the use of high-precision (i.e., quadruple or higher) arithmetic. In the current work, we avoid the use of high-precision arithmetic by using Riemann–Hilbert (RH) steepest-descent deformations [Reference Deift and Zhou6] to, in effect, precondition the singular integral equations. Despite the growing body of work on both theoretical and experimental fronts, an efficient framework for computing soliton gas solutions of integrable nonlinear wave models has been elusive.
Specifically, with this article, we introduce a fast and accurate method to compute a class of soliton gas primitive potentials in the context of the Korteweg–de Vries (KdV) equation
\begin{equation}
u_t + 6 u u_x + u_{xxx} = 0,\quad -\infty \lt x \lt +\infty,
\end{equation}based on the numerical solution of their Riemann–Hilbert Problem (RHP) representations. An advantage of this approach is that there is no time-stepping involved, as (x, t) enter the problem as (explicit) parameters. The method is easily parallelizable over values of (x, t).
The method presented here can compute primitive potentials of soliton gases (at t = 0, for all x) and their time evolution under the KdV flow for small values of t, still on the entire x-axis. In this regime, the method is asymptotically accurate as 
$x\to\pm \infty$ and does not require high-precision arithmetic. Importantly, it allows us to compute soliton gas solutions of the KdV equation with associated RH jump conditions supported on many pairs of bands. We also compute such soliton gas solutions nonlinearly superposed with a number of solitons. See Figure 1 for a soliton gas solution of (1.1) with associated RH jump conditions supported on five (pairs of) bands and nonlinearly superposed with five solitons, computed with our method. Pointwise evaluation of the solution in Figure 1 does not require higher precision arithmetic and takes about 0.7 s in the unmodulated region (see Section 2.2) and 0.08 s in the quiescent region (see Section 2.1) on a standard laptop. Solution animations and Julia code used to generate the plots in this paper can be found at [Reference Ballew, Bilman and Trogdon1].
Remark 1 [Hardware used]
All computations in this paper can be performed on a Lenovo laptop running Ubuntu version 20.04 with 8 cores and 16 GB of RAM with an Intel® Core™ i7-11800H processor running at 2.30 GHz. However, due to the ease of parallelization, a cluster can greatly reduce the time it takes to generate the figures.
The method presented here can also compute soliton gas solutions of the KdV equation in the entire (x, t)-plane outside of an unbounded wedge-shaped region emanating from 
$(x,t)=(0,0)$. See Figure 2 for a pure soliton gas (no solitons) supported on five pairs of bands. The bottom panel presents the computed large-time evolution in the tail 
$x \lt -K t$ for some suitable choice of constant K > 0. Figure 3 provides the density plot of a soliton gas with associated RH jump conditions supported on a single pair of bands, nonlinearly superposed with two solitons, in the complement of the aforementioned wedge-shaped region. As is apparent from this figure, the validity of our method extends a little bit into the wedge (see the boundary curves in black). The extension stops once the exponential factors supported on the suitable deformed jump contours become too large as (x, t) penetrates into the wedge.
Figure 1. KdV soliton gas with 
$r_1(\lambda)$ supported on five pairs of bands, nonlinearly superposed with five solitons. In the notation of (1.5), 
$I_1 = (0.25,0.5)$, 
$I_2 = (0.8,1.2)$, 
$I_3=(1.5,2)$, 
$I_4=(2.5,3)$, 
$I_5=(4,5)$, with 
$f_1(z)=1$, 
$f_2(z)=1/2$, 
$f_3(z)=1/4$, 
$f_4(z)=1/8$, 
$f_5(z)=1/16$ and 
$\alpha_j=\beta_j=\frac{1}{2}$ for 
$j=1,\ldots,5$. The solitons are associated with (see Riemann--Hilbert Problem 2) 
$\kappa_1=0.1$, 
$\kappa_2 = 0.7$, 
$\kappa_3=2.25$, 
$\kappa_4=3.5$, 
$\kappa_5 = 5.5$ with the norming constants 
$\chi_1=10^5$, 
$\chi_2=1000$, 
$\chi_3=100$, 
$\chi_4=10$ and 
$\chi_5=10^{-6}$.
Figure 2. A pure KdV soliton gas with 
$r_1(\lambda)$ supported on five pairs of bands. In the notation of (1.5), 
$I_1=(0.25,0.5)$, 
$I_2=(0.8,1.2)$, 
$I_3=(1.5,2)$, 
$I_4=(2.5,3)$, 
$I_5=(4,5)$ with 
$f_1(z)=1$, 
$f_2(z)=1/2$, 
$f_3(z)=1/4$, 
$f_4(z)=1/8$, 
$f_5(z)=1/16$ and 
$\alpha_j=\beta_j=\frac{1}{2}$ for 
$j=1,\ldots,5$.
Figure 3. Density plot of the computed soliton gas with 
$r_1(\lambda)$ supported on a single pair of bands with two solitons. In the notation of (1.5), 
$I_1=(1.5,2.5)$, 
$f_1(z)=1$ and 
$\alpha_1=\beta_1=\frac{1}{2}$. The solitons are associated with the eigenvalue parameters 
$\kappa_1=1$, 
$\kappa_2=4$ and the norming constants 
$\chi_1=10$, 
$\chi_2=10^{-10}$. Outside of the wedge region, the numerical method presented here is seen to be uniformly accurate with a computational cost that is independent of (x, t). Inside the wedge, the numerical method begins to break down, and additional RH deformations will need to be incorporated.
In Figure 4, we present a plot of the computed large-time evolution of a soliton gas with associated RH jump conditions supported on two pairs of bands (instead of five as in Figures 1 and 2), and in the notation of (1.5), we choose 
$f_1(z)$ to be much larger than 
$f_2(z)$.
Figure 4. A pure KdV soliton gas with 
$r_1(\lambda)$ supported on two pairs of bands. In the notation of (1.5), 
$I_1=(1,2)$, 
$I_2=(2.5,3)$ with 
$f_1(z)=100$, 
$f_2(z)=1$ and 
$\alpha_{j}=\beta_{j}=\frac{1}{2}$, 
$j=1,2$. Bottom panel: The same solution plotted along the ray 
$x/t = -32$.
This work is a culmination of recent advances made in computing orthogonal polynomials that are orthogonal on multiple disjoint intervals [Reference Ballew and Trogdon3] and work in computing large-genus solutions of the KdV equation [Reference Bilman, Nabelek and Trogdon4]. The key ingredient behind these developments is the choice of a basis that encodes the behaviour (e.g., the singularity structure) of the solution of an RHP at the endpoints of the jump contour. We also make use of the machinery available in the OperatorApproximation.jl software package [Reference Trogdon15].
This work focuses mainly on the computation of soliton gas potentials supported on many bands and completing the groundwork for nonlinear superpositions of such potentials with many solitons. Extensions of our method to the entire (x, t)-plane, along with routines to compute a broader class of soliton gas solutions, will appear in a forthcoming article.
1.1. RH definition of a KdV soliton gas
Let 
$n\in\mathbb{Z}_+$ be fixed and consider a collection of disjoint intervals 
$\mathrm{i}(a_j, b_j)$, 
$1\leq j \leq n$, on the imaginary axis with 
$0 \lt a_j \lt b_j \lt a_{j+1}$ for each 
$1\leq j \leq n-1$. We denote the union of these intervals and of their reflections by
\begin{equation}
\Sigma_+ := \bigcup_{j=1}^n (a_j, b_j)\qquad \text{and} \qquad \Sigma_- := \bigcup_{j=1}^n (-b_j, -a_j),
\end{equation} and take all of these intervals to be oriented upwards. We consider a function 
$r_1 \colon (\mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-) \to \mathbb{C}$ such that 
$r_1(\lambda)$ is positive for λ on each interval 
$\mathrm{i} (a_j, b_j)$, extending analytically to a neighbourhood of each 
$\mathrm{i} (a_j, b_j)$ with the local behaviour
\begin{equation}
r_1(\mathrm{i} z) = f_{a_j}(z) (z-a_j)^{\alpha_j},\quad \alpha_j\in\left\{-\tfrac{1}{2},\tfrac{1}{2}\right\},
\end{equation} for z in a neighbourhood of aj, where 
$f_{a_j}(z)$ is analytic and satisfies 
$f_{a_j}(z) \gt 0$ for 
$z\geq a_j$. Similarly,
\begin{equation}
r_1(\mathrm{i} z) = f_{b_j}(z) (b_j-z)^{\beta_j},\quad \beta_j\in\left\{-\tfrac{1}{2},\tfrac{1}{2}\right\},
\end{equation} for z in a neighbourhood of bj, where 
$f_{b_j}(z)$ is analytic and satisfies 
$f_{b_j}(z) \gt 0$ for 
$z\leq b_j$. Here, all of the roots are taken to be the principal branch. Therefore, in the computations performed for this paper, we take
\begin{equation}
r_1(\mathrm{i} z) = f_j(z)(z-a_j)^{\alpha_j}(b_j-z)^{\beta_j},\quad \text{for}~z\in I_j:=(a_j,b_j), \quad \Sigma_+=\bigcup_{j=1}^n I_j,
\end{equation} for analytic functions fj and constants 
$\alpha_j,\beta_j\in\left\{-\tfrac{1}{2},\tfrac{1}{2}\right\}$. For 
$\lambda\in\mathrm{i}\Sigma_-$, r 1 is defined by symmetry: 
$r_1(\lambda)=r_1(-\lambda)$. In what follows 
$\mathrm{cl}(\Sigma)$ denotes the closure of a set Σ. With these ingredients, we consider the following problem, which we take as the definition of a soliton gas:
Riemann–Hilbert Problem 1 (Pure KdV soliton gas)
Let 
$(x,t)\in\mathbb{R}\times\mathbb{R}_{+}$ be fixed. Find a 
$1\times 2$ row vector-valued function 
$\mathbf{M}(\lambda)\equiv \mathbf{M}(\lambda;x,t)$ with the following properties:
• Analyticity:
$\mathbf{M}(\lambda)$ is analytic for
$\lambda\in\mathbb{C}\setminus\mathrm{cl}( \mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-)$, admitting continuous boundary values on
$\mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-$.• Jump conditions: The boundary values
$\mathbf{M}^+(\lambda)$ (resp.
$\mathbf{M}^-(\lambda)$) from the left (resp. right) are related via
(1.6)\begin{align}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\mathbf{M}^+(\lambda)=\mathbf{M}^-(\lambda)\begin{bmatrix} 1 & 0 \\ -2\mathrm{i} r_1(\lambda)\mathrm{e}^{2\mathrm{i}(x \lambda + 4 t \lambda^3)} & 1\end{bmatrix},\quad \lambda\in \mathrm{i}(a_j,b_j),
\end{align}(1.7)\begin{align}
\mathbf{M}^+(\lambda)&=\mathbf{M}^-(\lambda)\begin{bmatrix} 1 & 2\mathrm{i} r_1(\lambda)\mathrm{e}^{-2\mathrm{i}(x\lambda + 4 t \lambda^3)} \\ 0 & 1\end{bmatrix},\quad \lambda\in \mathrm{i}(-b_j,-a_j), \quad j=1,2,\ldots,n.
\end{align}• Normalization:
$\mathbf{M}(\lambda)=\begin{bmatrix} 1 & 1\end{bmatrix} + \mathrm{O}(\lambda^{-1})$ as
$\lambda\to\infty$.• Symmetry condition:
$\mathbf{M}(-\lambda)=\mathbf{M}(\lambda) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$.
This problem has appeared in [Reference Dyachenko, Zakharov and Zakharov7, Reference Girotti, Grava, Jenkins and McLaughlin9] with different assumptions on the behaviour of 
$r_1(\lambda)$ at the endpoints. A derivation of Riemann–Hilbert Problem 1 via a limiting procedure involving the accumulation of eigenvalues of the Schrödinger operator on a single pair of bands with respect to a suitable density can be found in [Reference Girotti, Grava, Jenkins and McLaughlin9]. The setting of [Reference Girotti, Grava, Jenkins and McLaughlin9] assumes that 
$r_1(\lambda)$ is nonvanishing and bounded at the endpoints, in contrast with the assumptions in our work. As the reader will see, the endpoint behaviour that we assume for 
$r_1(\lambda)$ enables a particularly fast computational method. Riemann–Hilbert Problem 1 has a unique solution [Reference Girotti, Grava, Jenkins and McLaughlin9] (see Section 1.2 for the notion of a solution when αj or βj is negative), and we can recover the soliton gas solution of (1.1) by the second-order term at infinity [Reference Egorova, Piorkowski and Teschl8]
\begin{equation}
u(x,t) = -\lim_{\lambda\to\infty} 2\lambda^2(m_1(z;x,t) m_2(z;x,t) - 1),
\end{equation} where 
$m_j(z;x,t)$ is the jth component of the row vector 
$\mathbf{M}(\lambda;x,t)$.
We proceed with rotating the λ-plane clockwise by introducing 
$r(z) :=2 r_1(\mathrm{i} z)$ and 
$\mathbf{Y}(z):=\mathbf{M}(\mathrm{i} z)$, which satisfies the jump conditions
\begin{equation}
\mathbf{Y}^+(z)=\mathbf{Y}^-(z) \begin{bmatrix} 1 & 0 \\ -\mathrm{i} r(z)\mathrm{e}^{-2\theta(z;x,t)} & 1\end{bmatrix},\quad z\in \Sigma_+,\qquad \theta(z;x,t):= xz - 4 tz^3,
\end{equation}
\begin{equation}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\mathbf{Y}^+(z)=\mathbf{Y}^-(z)\begin{bmatrix} 1 & \mathrm{i} r(z)\mathrm{e}^{2\theta(z;x,t)} \\ 0 & 1\end{bmatrix},\quad\quad\,\,\! z\in \Sigma_-,
\end{equation} and is analytic elsewhere. Note that all (real) intervals are oriented in the increasing direction. We consider the numerical solution of the RHP satisfied by 
$\mathbf{Y}(z)$, which inherits the normalization and symmetry conditions in Riemann--Hilbert Problem 1.
Regarding the behaviour of 
$\mathbf{Y}(z)$ at the endpoints, let cj denote either of the endpoints aj or bj, and let γj denote the corresponding power αj or βj as in (1.3) and (1.4). As 
$z\to c_j$, we have
\begin{equation}
\mathbf{Y}(z) =
\begin{cases}
\begin{bmatrix} \mathrm{O}(1) & \mathrm{O}(1) \end{bmatrix}, &\quad \text{if} \quad \gamma_j = \frac{1}{2},\\[0.5em]
\begin{bmatrix} \mathrm{O}(|z-c_j|^{-\frac{1}{2}}) & \mathrm{O}(1) \end{bmatrix}, &\quad \text{if} \quad\gamma_j = -\frac{1}{2}.
\end{cases}
\end{equation} This situation is mirrored on the negative real axis. As 
$z\to-c_j$, we have
\begin{equation}
\mathbf{Y}(z) =
\begin{cases}
\begin{bmatrix} \mathrm{O}(1) & \mathrm{O}(1) \end{bmatrix}, &\quad \text{if} \quad \gamma_j = \frac{1}{2},\\[0.5em]
\begin{bmatrix} \mathrm{O}(1) &\mathrm{O}(|z+c_j|^{-\frac{1}{2}})\end{bmatrix}, &\quad \text{if} \quad\gamma_j = -\frac{1}{2}.
\end{cases}
\end{equation} Therefore, there is no ambiguity in the statement of the RHP for 
$\mathbf{Y}(z)$ (or in Riemann--Hilbert Problem 1, for that matter) arising from unspecified behaviour of the solution at the endpoints of the jump contour — it is implied by the behaviour of the jump matrix at the endpoints.
Writing 
$\lambda = \mathrm{i} z$ in the large-λ expansions that led to (1.8), the soliton gas solution 
$u(x,t)$ of (1.1) is recovered from 
$\mathbf{Y}(z)\equiv \mathbf{Y}(z;x,t)$ via the formula
\begin{equation}
u(x,t) = \lim_{z\to\infty} 2 z^2 (y_1(z;x,t)y_2(z;x,t) - 1).
\end{equation}1.2. From local solutions to consistent numerical approximations
A typical approach to the numerical solution of Riemann--Hilbert Problem 1, for example, is to seek a singular integral equation that the Cauchy density solves [Reference Olver11, Reference Trogdon and Olver14]. More specifically, if
\begin{align}
\mathbf{Y}^+(z) = \mathbf{Y}^-(z) \mathbf{J}(z), \quad z \in \Gamma,
\end{align}then we write
\begin{align}
\mathbf{Y}(z) = \begin{bmatrix} 1 & 1 \end{bmatrix} + \mathcal C_{\Gamma} \mathbf{ U}(z), \quad \mathcal C_{\Gamma}\mathbf{ U}(z) = \frac{1}{2 \pi \mathrm{i}} \int_{\Gamma} \frac{\mathbf{ U}(\zeta)}{\zeta - z} \mathrm{d}\zeta.
\end{align} If we denote the boundary values of the Cauchy transform by 
$\mathcal C_{\Gamma}^\pm$, 
$\mathbf{ U}$ satisfies the singular integral equation
\begin{align}
\mathcal C_{\Gamma}^+\mathbf{ U}(z) - \mathcal C_{\Gamma}^- \mathbf{ U}(z) \mathbf{J}(z) = \begin{bmatrix} 1 & 1 \end{bmatrix}(\mathbf{J}(z) - \mathbb{I}), \quad z \in \Gamma.
\end{align} To construct an efficient numerical method for this equation, one must choose a basis in which 
$\mathbf{ U}$ can be accurately represented, something that will typically result in the consistency of the numerical method. To illustrate this, we consider a toy problem where 
$\Gamma = \Gamma_1 \cup \Gamma_2$ are disjoint, 
$\Gamma_1 = [-1,1]$, oriented from 
$z=-1$ to z = 1, and
\begin{align}
\mathbf{J}|_{\Gamma_1}(z) = \begin{bmatrix} 0 & \mathrm{i} h(z) \\ \mathrm{i} /h(z) & 0 \end{bmatrix}, \quad h(z) = \sqrt{1 -z^2 }.
\end{align}We will build a local solution to this problem near z = 1, and it will tease out the singularity structure. Technically speaking, we will also need to verify that our local solution satisfies the requisite endpoint conditions, see (2.9), for example. We use this example because it represents the structure encountered in the unmodulated region of Section 2.2. For that encountered in the quiescent region of Section 2.1, the same analysis can be performed, but it is simpler because the jump matrices are triangular and local solutions are found directly by Cauchy integrals.
Consider 
$\ell(z) = (z-1)^{1/4} (z + 1)^{1/4}$, where we use the principal branch cut along 
$(-\infty,0]$ for the quarter root function. It follows that this function satisfies 
$\ell^+(z) \ell^-(z) = h(z)$, 
$z \in [-1,1]$. Then, define
\begin{align}
\mathbf{ L}(z) = \ell(z)^{-\sigma_3}, \quad \sigma_3 = \mathrm{diag}(1,-1).
\end{align}We compute, in a neighbourhood of z = 1,
\begin{align}
\mathbf{ L}^-(z)\mathbf{J}|_{\Gamma_1}(z) \mathbf{ L}^+(z)^{-1} = \begin{bmatrix} 0 & \mathrm{i} \frac{h(z)}{\ell^+(z) \ell^-(z)} \\ \mathrm{i} \frac{\ell^+(z) \ell^-(z)}{h(z)} & 0 \end{bmatrix} = \begin{bmatrix} 0 & \mathrm{i} \\ \mathrm{i} & 0 \end{bmatrix}.
\end{align}Then, we factorize
\begin{align}
\begin{bmatrix} 0 & \mathrm{i} \\ \mathrm{i} & 0 \end{bmatrix} = \mathbf{Q} a^-(z)^{-\sigma_3} a^+(z)^{\sigma_3} \mathbf{Q}^{-1}, \quad \mathbf{Q} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}, \quad a(z) = \left( \frac{z -1}{z+ 1}\right)^{1/4},
\end{align}giving
\begin{align}
a^-(z)^{\sigma_3}\mathbf{Q}^{-1} \mathbf{ L}^-(z)\mathbf{J}|_{\Gamma_1}(z) \mathbf{ L}^+(z)^{-1}\mathbf{Q} a^+(z)^{-\sigma_3} = \mathbb{I}.
\end{align}All of this is to say that
\begin{align}
\boldsymbol{\Psi}(z) = a(z)^{\sigma_3} \mathbf{Q}^{-1} \mathbf{ L}(z) = \frac{1}{\sqrt{2}}\begin{bmatrix}
\frac{1}{\sqrt{z+1}} & \sqrt{z-1} \\
-\frac{1}{\sqrt{z-1}} & \sqrt{z+1} \\
\end{bmatrix},
\end{align} is a local solution to the RHP near z = 1, where the square roots denote the principal branch cut along 
$(-\infty,0]$. Then define, for ϵ small,
\begin{align}
\widehat{\mathbf{Y}}(z) = \begin{cases} \mathbf{Y}(z) & |z-1| \gt \epsilon,\\
\mathbf{Y}(z) \boldsymbol{\Psi}({z})^{-1} & |z-1| \lt \epsilon,\end{cases}
\end{align}which will be an analytic function near z = 1 (using the requisite endpoint conditions). We find the representation
\begin{align}
\mathbf{Y}(z) = \widehat{\mathbf{Y}}(z) \boldsymbol{\Psi}(z),
\end{align} for z near z = 1. To back out the behaviour of 
$\mathbf{ U}$, we use the Plemelj lemma,
\begin{align}
\mathbf{ U}(z) = \mathbf{Y}^+(z) - \mathbf{Y}^-(z) = \widehat{\mathbf{Y}}(z) \left(\boldsymbol{\Psi}^+(z) - \boldsymbol{\Psi}^-(z)\right) = \widehat{\mathbf{Y}}(z) \begin{bmatrix} 0& \sqrt{2(z-1)} \\ -\sqrt{\frac{2}{z-1}} & 0 \end{bmatrix}, \quad |z-1| \lt \epsilon.
\end{align} A similar analysis can be performed near 
$z = -1$, and patching together these arguments, it follows that
\begin{align}
\mathbf{ U}|_{\Gamma_1}(z) = \mathbf{A}(z)\begin{bmatrix} 0& \sqrt{1-z^2}\\ \frac{1}{\sqrt{1-z^2}} &0\end{bmatrix},
\end{align} where 
$\mathbf{A}(z)$ is an analytic vector-valued function that can therefore be approximated well with polynomials of low degree. This local analysis informs the ansatz (2.2) below.
2. Computing soliton gases
2.1. Quiescent region
When t = 0, the jump matrices in (1.9)–(1.10) are 
$\mathrm{O}(1)$ (away from the endpoints) for 
$x\geq 0$ and become exponentially close to identity as 
$x\to+\infty$ since 
$\textrm{Re}(-2 z x) \lt 0$ on 
$\Sigma_+$ and 
$\textrm{Re}(2 z x) \lt 0$ on 
$\Sigma_-$ in this case. This configuration extends to a region of the (x, t)-plane for t > 0 as long as 
$\textrm{Re}(-2\theta(z;x,t)) \lt 0$ is maintained on 
$\Sigma_+$ (the situation is mirrored on 
$\Sigma_-$ automatically). The axis 
$\textrm{Re}(z)=0$ is always a part of the locus 
$\textrm{Re}(-2\theta(z;x,t))=0$. The other branch of this locus is given by 
$\textrm{Re}(z)^2 = 3\,\textrm{Im}(z)^2 + \frac{x}{4t}$. Thus, the aforementioned boundedness or exponential decay of the jump matrices to the identity is preserved for all (x, t) with 
$t \geq 0$ satisfying
\begin{equation}
x \geq 4 K t,
\end{equation} for an arbitrary constant 
$K \gt b_n^2$. In this setting, the RHP at hand can be treated numerically as is, without implementing any contour deformations, by using an appropriate basis with weights encoding the endpoint behaviour of 
$\mathbf{Y}(z)$. We employ the approach put forth in [Reference Ballew and Trogdon3, Reference Bilman, Nabelek and Trogdon4] and consider the ansatz given by the superposition of weighted Cauchy transforms
\begin{equation}
\mathbf{Y}(z) = \begin{bmatrix} 1 & 1 \end{bmatrix} + \sum_{j=1}^n\left[ \int_{a_j}^{b_j} \frac{\mathbf{F}_{+j} (\zeta) {\begin{bmatrix}w_{+j}^{[1]}(\zeta)& 0 \\ 0 & w_{+j}^{[2]}(\zeta)\end{bmatrix}}\mathrm{d}\zeta}{\zeta-z} + \int_{-b_j}^{-a_j} \frac{\mathbf{F}_{-j} (\zeta) {\begin{bmatrix}w_{-j}^{[1]}(\zeta)& 0 \\ 0 & w_{-j}^{[2]}(\zeta)\end{bmatrix}}\mathrm{d}\zeta}{\zeta-z}\right],
\end{equation} for unknown row vector-valued densities 
$\mathbf{F}_{\pm j}(\zeta)$ on each interval. Here, each of the weights 
$w_{\pm j}^{[1]}(\zeta)$ and 
$w_{\pm j}^{[2]}(\zeta)$ (used for the first and second columns, respectively) is equal to one of the Chebyshev weights
\begin{align*}
\rho_1(\zeta)&:= \frac{1}{\sqrt{(\zeta-a)(b-\zeta)}},
\quad \rho_2(\zeta):=\rho_1(\zeta)^{-1}, \\
\rho_3(\zeta)&:=\sqrt{(\zeta-a)(b-\zeta)^{-1}},
\quad \rho_4(\zeta):=\rho_3(\zeta)^{-1},
\end{align*} for the four kinds of Chebyshev orthogonal polynomials on an interval 
$[a,b]$ on 
$[a_j, b_j]$ for 
$w_{+j}^{[1],[2]}(\zeta)$ and on 
$[-b_j,-a_j]$ for 
$w_{-j}^{[1],[2]}(\zeta)$. The choice is made so that the weight captures the behaviour of the given r(z) (and hence that of 
$\mathbf{Y}(z)$) at the endpoints, and the corresponding Chebyshev polynomials are used as the basis on the relevant interval. The associated singular integral equation is discretized and solved by collocation, as described in [Reference Olver11, Reference Trogdon and Olver14], by enforcing that the singular integral equation should be satisfied exactly at mapped Chebyshev first-kind zeros.
Remark 2.1. The triangular nature of the jump matrices (1.9)–(1.10) implies that the second column of 
$\mathbf{F}_{+j}(\zeta)$ and the first column of 
$\mathbf{F}_{-j}(\zeta)$ are identically equal to 0. Therefore, the weights 
$w^{[2]}_{+j}(\zeta)$ and 
$w^{[1]}_{-j}(\zeta)$ are immaterial. In practice, we do not encode this structure and allow the linear system arising from collocation to force those entries to be equal to 0. This structure is lost in other asymptotic regions upon contour deformations; therefore, we present the ansatz (2.2) in its general form for future reference.
Finally, note that since 
$r(z)=r(-z)$, 
$w_{-j}^{[2]}(\zeta) = w_{+j}^{[1]}(-\zeta)$ for 
$\zeta\in[-b_j, -a_j]$. The region described by (2.1) is the analogue of the ‘constant region’ as 
$t\to+\infty$ in [Reference Girotti, Grava, Jenkins and McLaughlin9] for the case n = 1.
2.2. Unmodulated region
When t = 0 and x < 0, the situation is dramatically different: All of the exponentials in (1.9)–(1.10) grow exponentially and unboundedly as 
$x\to-\infty$. This behaviour cannot be put under control solely by employing one of the four canonical matrix factorizations for the jump matrices since the exponents are real-valued. In this setting, the RHP is stabilized by modifying the exponent via introducing a g-function. This is a generalization of the approach taken in [Reference Girotti, Grava, Jenkins and McLaughlin9, Reference Girotti, Grava, Jenkins, McLaughlin and Minakov10] for the case n = 1. Let 
$\Gamma_{\pm j}$ denote the gaps: 
$\Gamma_0=[-a_1, a_1]$, 
$\Gamma_{j}=[b_{j},a_{j+1}]$ and 
$\Gamma_{-j}=[-a_{j+1},-b_{j}]$ for 
$j=1,2,\ldots,n-1$. To have a uniform treatment that also captures a region when t > 0, we proceed by allowing for 
$t\geq 0$ with x < 0 and do not designate a particular large parameter. We seek a function 
$g(z)\equiv g(z;x,t)$ analytic for 
$z\notin [-b_n , b_n]$, with the following additional properties: g(z) admits continuous boundary values 
$g^{\pm}(z)$ on 
$[-b_n,b_n]$, taken from 
$\mathbb{C}^{\pm}$, satisfying
\begin{align}
g^+(z)+g^-(z) &= 2 \theta(z;x,t), \quad\! z\in \Sigma_+ \cup \Sigma_-, \nonumber \\
g^+(z)-g^-(z) &= \mathrm{i} \Omega_{\pm j}(x,t), \quad z\in \Gamma_{\pm j},\quad j=0,1,2,\ldots,n-1,
\end{align} where 
$\Omega_{\pm j}(x,t)$ are real-valued constants. Finally, 
$g(z)=\mathrm{O}(z^{-1})$ as 
$z\to\infty$. Let R(z) be the function analytic for 
$z\notin \mathrm{cl}(\Sigma_+ \cup \Sigma_-)$ satisfying 
$
R(z)^2 = \prod_{j=1}^n (z^2-a_j^2)(z^2-b_j^2)
$ along with 
$R(z)=z^{2n} + \mathrm{O}(z^{2n-2})$ as 
$z\to\infty$. Then,
\begin{equation}
g(z) = \frac{R(z)}{2\pi \mathrm{i}} \int\limits_{\Sigma_{+} \cup \Sigma_{-}} \frac{2\theta(\zeta;x,t) \mathrm{d}\zeta}{R^+(\zeta)(\zeta-z)} + \sum_{\ell=-(n-1)}^{n-1}\frac{\Omega_\ell(x,t)R(z)}{2\pi} \int\limits_{\Gamma_\ell} \frac{\mathrm{d} \zeta}{R(\zeta)(\zeta - z)}\,,
\end{equation} where the 
$2n-1$ constants 
$\Omega_{\ell}(x,t)$ are chosen to ensure the desired decay 
$g(z)=\mathrm{O}(z^{-1})$ as 
$z\to\infty$. These moment conditions result in a linear system:
\begin{equation}
\sum_{\ell=-(n-1)}^{n-1}\mathrm{i}\Omega_{\ell}(x,t) \int\limits_{\Gamma_\ell} \frac{\zeta^k \mathrm{d} \zeta}{R(\zeta)} = -2\int\limits_{\Sigma_-\cup\Sigma_+} \frac{\zeta^k \theta(\zeta;x,t)}{R^+(\zeta)} \mathrm{d}\zeta,\quad k =0,1,2,\ldots, 2n-2.
\end{equation} While this linear system can be shown to have a unique solution, it becomes ill-conditioned for large values of n due to the presence of the monomials ζ k. We solve an equivalent but empirically well-conditioned linear system obtained by replacing the monomials ζ k with a suitable basis of polynomials of degree 
$2n-2$, chosen to vanish at the midpoints of the 
$2n-2$ gaps. That is, we compute the constants 
$\Omega_\ell(x,t)$ using the basis
\begin{equation*}
p_{k}(x)=\prod_{j=-n+1, j\neq k}^{n-1}(x-\mu_j), \quad k=-n+1,\ldots,n-1,
\end{equation*} where µj denotes the midpoint of 
$\Gamma_j$. We compute the integrals in the linear system (2.6) via Gauss–Chebyshev quadrature.
Now, introduce open and disjoint disks Dj, 
$j=1,2,\ldots,n$, chosen so that Dj encloses 
$[a_j,b_j]$ in its interior, and let D −j be the analogue of Dj for 
$[-b_j, -a_j]$. We take the disk boundaries to be oriented counter-clockwise and make the following definitions (with correlated signs on each line):
\begin{equation}
\ \ \!\!\!\!\!\!\mathbf{S}(z) := \mathbf{Y}(z)
\begin{bmatrix}
1 & \frac{\mathrm{i}}{r(z)}\mathrm{e}^{2\theta(z;x,t)} \\ 0 & 1
\end{bmatrix}^{\mp 1} \mathrm{e}^{g(z;x,t)\sigma_3},\qquad\ \ z\in D_j \cap \mathbb{C}^{\pm},
\end{equation}
\begin{equation}
\mathbf{S}(z) := \mathbf{Y}(z)
\begin{bmatrix}
1 &0 \\ \frac{-\mathrm{i}}{r(z)}\mathrm{e}^{-2\theta(z;x,t)} & 1
\end{bmatrix}^{\mp 1} \mathrm{e}^{g(z;x,t)\sigma_3},\qquad z\in D_{-j} \cap \mathbb{C}^{\pm},
\end{equation} for 
$j=1,2,\ldots,n$. Recall that 
$\sigma_3 = \mathrm{diag}(1,-1)$. We set 
$\mathbf{S}(z):= \mathbf{Y}(z)\mathrm{e}^{g(z;x,t)\sigma_3}$ elsewhere. Using 
$r^+(z)=-r^-(z)$ on 
$\Sigma_+\cup\Sigma_-$, we find that the jump conditions satisfied by 
$\mathbf{S}(z)$ are as given in Figure 6. We denote the modified exponent by 
$\varphi(z;x,t):=g(z;x,t)-\theta(z;x,t)$.
A contour plot of the sign of 
$\mathrm{Re}(\varphi(z;x,t))$ is given in Figure 5, and one can see that the jump matrices on the circles are exponentially close to the identity matrix if t > 0 becomes large for 
$x/t$ in this region. This behaviour together with the fact that the circular jump contours are detached from 
$\Sigma_+ \cup \Sigma_-$ (no self-intersection points with jump matrices involving r(z)) is what enables an efficient numerical method. Therefore, the assumptions on the behaviour of r(z) at the endpoints are absolutely essential for our approach.
Again, let cj denote either of the endpoints aj or bj, and let γj denote the corresponding power αj or βj as in (1.3) and (1.4). From (1.11), (2.7) and (2.8), we find that as 
$z\to c_j$,
\begin{equation}
\mathbf{S}(z) =
\begin{cases}
\begin{bmatrix} \mathrm{O}(1) & \mathrm{O}(|z-c_j|^{-\frac{1}{2}}) \end{bmatrix}, &\quad \text{if} \quad \gamma_j = \frac{1}{2},\\[0.5em]
\begin{bmatrix} \mathrm{O}(|z-c_j|^{-\frac{1}{2}}) & \mathrm{O}(1) \end{bmatrix}, &\quad \text{if} \quad\gamma_j = -\frac{1}{2}.
\end{cases}
\end{equation}
Figure 5. Contour plot showing the sign of 
$\mathrm{Re}(\varphi(z;x,t))$ in the complex plane for t > 0 in the unmodulated region treated in Section 2.2. Here, 
$x=-2$, t = 0.01, and in the notation of (1.5), 
$I_1=(1,2)$, 
$I_2=(2.5,3)$.
Figure 6. Jump contours and jump matrices associated with 
$\mathbf{S}(z;x,t)$ near each interval. 
$\varphi(z;x,t) = g(z;x,t) - \theta(z;x,t)$.
The situation is again mirrored as z approaches an endpoint of 
$[-b_j,-a_j]$, but the behaviour of the columns is flipped.
As in [Reference Ballew and Trogdon3, Section 3.2.3], we now introduce a correction to the g-function to eliminate the constant jump conditions on the gaps 
$\Gamma_{\pm j}$, 
$j=0,1,\ldots n-1$ at the expense of introducing constants to the jump matrices on the intervals 
$\Sigma_+\cup\Sigma_-$. We seek a function 
$h(z)\equiv h(z;x,t)$ analytic for 
$z\notin [-b_n,b_n]$ with the properties
\begin{equation} \begin{array}{lll} h^+(z) - h^-(z) = \log(\exp(-\mathrm{i} \Omega_{\pm k}(x,t))),& z\in \Gamma_{\pm k}, & k=0,1,2,\ldots, n-1,\\ h^+(z) + h^-(z) = A_{+j}(x,t), & z\in(a_j, b_j), & j=1,2,\ldots, n,\\ h^+(z) + h^-(z) = A_{-j}(x,t), & z\in(-b_j, -a_j),& j=1,2,\ldots, n, \end{array} \end{equation} along with 
$h(z)=\mathrm{O}(z^{-1})$ as 
$z\to\infty$. Then,
\begin{align}
h(z)=&\frac{R(z)}{2\pi\mathrm{i}}
\sum_{j=1}^{n}\left[\int_{a_j}^{b_j}\frac{A_{j+}(x,t)\mathrm{d} \zeta}{R^+(\zeta)(\zeta-z)}+\int_{-b_j}^{-a_j}\frac{A_{j-}(x,t)\mathrm{d} \zeta}{R^+(\zeta)(\zeta-z)}\right]\\\notag &+\frac{R(z)}{2\pi\mathrm{i}}\sum_{\ell=-(n-1)}^{n-1}\int\limits_{\Gamma_\ell} \frac{\log(\exp(-\mathrm{i}\Omega_\ell(x,t)))\mathrm{d}\zeta}{R(\zeta)(\zeta-z)},
\end{align}
where the constants 
$A_{\pm j}(x,t)$ are determined to ensure the desired decay at infinity. This results in the following linear system of 2n equations:
\begin{equation}
\sum_{j=1}^{n} \left[
A_{+j}(x,t)\int_{a_j}^{b_j}\frac{\zeta^k \mathrm{d}\zeta}{R^+(\zeta)} + A_{-j}(x,t)\int_{-b_j}^{-a_j}\frac{\zeta^k \mathrm{d}\zeta}{R^+(\zeta)} \right]
= - \sum_{\ell=-(n-1)}^{n-1}\log(\exp(-\mathrm{i}\Omega_\ell(x,t)))\int\limits_{\Gamma_\ell} \frac{\zeta^k \mathrm{d}\zeta}{R(\zeta)},
\end{equation} indexed by 
$k=0,1,2,\ldots,2n-1$. Note that the right-hand side of (2.11) and the coefficients of the unknowns 
$A_{\pm j}(x,t)$ are all purely imaginary for each k. Therefore, 
$A_{\pm j}(x,t)$ are real-valued. This is alarming at first; however, the coefficient matrix for this linear system is independent of (x, t), and the right-hand side is bounded in (x, t) (even though 
$\Omega_\ell(x,t) \in\mathbb{R}$ may, in principle, grow unboundedly). This structure ensures that 
$A_{\pm j}(x,t)\in\mathbb{R}$ are bounded in (x, t). The integrals that make up the linear system (2.11) are again computed numerically via Gauss–Chebyshev quadrature and by employing a suitable basis of polynomials. In particular, we now use
\begin{equation*}
q_{k}(x)=\prod_{j=-n, j\neq 0,k}^{n}(x-\nu_j), \quad k=\pm1,\ldots,\pm n,
\end{equation*} where νj denotes the midpoint of the jth band (
$(a_j,b_j)$ if j > 0 or 
$(-b_{-j},-a_{-j})$ if j < 0).
We now make a global definition and introduce an exponent:
\begin{equation}
\mathbf{R}(z):= \mathbf{S}(z)\mathrm{e}^{h(z;x,t)\sigma_3}\quad \text{and}\quad \phi(z;x,t):=h(z;x,t) + g(z;x,t) - \theta(z;x,t).
\end{equation} Note that 
$\mathbf{R}(z) = \begin{bmatrix}1& 1 \end{bmatrix} + \mathrm{O}(z^{-1})$ as 
$z\to\infty$. Since the transformations 
$\mathbf{Y}(z)\mapsto\mathbf{S}(z)\mapsto\mathbf{R}(z)$ involve right-multiplications by diagonal unimodular matrices near 
$z=\infty$, it follows from (1.13) that
\begin{equation}
u(x,t) = \lim_{z\to\infty} 2z^2(r_1(z;x,t) r_2(z;x,t) - 1).
\end{equation} If 
$z\in\Gamma_\ell$, then 
$h^+(z;x,t) - h^-(z;x,t) + \mathrm{i}\Omega_{\ell}(x,t) = \log(\exp(-\mathrm{i} \Omega_{\ell}(x,t)) +\mathrm{i}\Omega_{\ell}(x,t) = \mathrm{i} 2 k \pi$ for some 
$k\in\mathbb{Z}$, so 
$\mathbf{R}(z)$ has no jumps across the gaps 
$\Gamma_\ell$. The jump conditions satisfied by 
$\mathbf{R}(z)$ are described in Figure 7. This final RHP satisfied by 
$\mathbf{R}(z)$ is again treated numerically by using appropriate basis functions and an ansatz involving appropriate weighted Cauchy transforms, similar to (2.2). More specifically, we use Laurent polynomials for the components on the circular contours and appropriate Chebyshev polynomials and their weight functions for the components on 
$\Sigma_+ \cup \Sigma_-$, capturing the endpoint behaviour of 
$\mathbf{R}(z)$, which is exactly the same as in (2.9). This approach closely follows the computational framework in [Reference Ballew and Trogdon3].
Remark 2.2. This method could, in principle, be sped up significantly by considering r with a common analytic extension to a region containing 
$[-b_n,b_n]$. In this case, the 2n circular jump contours can be replaced by a single large one. See [Reference Ballew and Trogdon2, Appendix A].
Note that the gap in the (x, t)-plane for t > 0 between the regions described in Section 2.1 and Section 2.2 disappears at t = 0, allowing us to compute the soliton gas primitive potentials quickly and accurately for all 
$x\in\mathbb{R}$.
Figure 7. Jump contours and jump matrices associated with 
$\mathbf{R}(z;x,t)$ near each interval. 
$\phi(z;x,t) = \varphi(z;x,t) + h(z;x,t)$.
3. Nonlinear superposition of a pure soliton gas with solitons
It is possible to insert a multi-soliton into a soliton gas by supplying suitable residue conditions for the RHP satisfied by 
$\mathbf{Y}(z)$ for a number of simple poles. In this case, we consider solving the following RHP:
Riemann–Hilbert Problem 2 (KdV soliton gas with solitons)
Let 
$(x,t)\in\mathbb{R}\times\mathbb{R}_{+}$ be fixed. Find a 
$1\times 2$ (row vector-valued) function 
$\mathbf{N}(\lambda)\equiv \mathbf{N}(\lambda;x,t)$ with the following properties:
• Analyticity:
$\mathbf{N}(\lambda)$ is analytic for
$\lambda\in\mathbb{C}\setminus\mathrm{cl}( \mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-)$, admitting continuous boundary values on
$\mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-$, with the exception of simple poles at
$\lambda=\pm\mathrm{i} \kappa_j$,
$\kappa_j \gt 0$,
$j=1,2,\ldots,N$, in the complement of the collection of intervals
$\mathrm{cl}( \mathrm{i}\Sigma_+ \cup \mathrm{i} \Sigma_-)$.• Jump conditions: The boundary values
$\mathbf{N}^+(\lambda)$ (resp.
$\mathbf{N}^-(\lambda)$) from left (resp. from right) are related via
(3.1)\begin{equation}
\!\!\!\!\!\!\!\ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\mathbf{N}^+(\lambda)=\mathbf{N}^-(\lambda)\begin{bmatrix} 1 & 0 \\ -2\mathrm{i} r_1(\lambda)\mathrm{e}^{2\mathrm{i}(x \lambda + 4 t \lambda^3)} & 1\end{bmatrix},\quad \lambda\in \mathrm{i}(a_j,b_j),
\end{equation}(3.2)\begin{equation}
\mathbf{N}^+(\lambda)=\mathbf{N}^-(\lambda)\begin{bmatrix} 1 & 2\mathrm{i} r_1(\lambda)\mathrm{e}^{-2\mathrm{i}(\lambda x+ 4 t \lambda^3)} \\ 0 & 1\end{bmatrix},\quad \lambda\in \mathrm{i}(-b_j,-a_j), \quad j=1,2,\ldots,n.
\end{equation}• Residue conditions: At each simple pole
$\lambda=\pm\mathrm{i}\kappa_j$,
$j=1,2,\ldots,N$,
$\mathbf N(\lambda)$ satisfies
(3.3)\begin{equation}\ \ \ \ \ \ \ \underset{\lambda=\mathrm i\kappa_j}{\mathrm{Res}}\ \mathbf N(\lambda)=\lim_{\lambda\rightarrow\mathrm i\kappa_j}\mathbf N(\lambda)\begin{bmatrix}0&0\\\mathrm i\chi_j\mathrm e^{2\mathrm i(\mathrm i\kappa_jx-4\mathrm i\kappa_j^3t)}&0\end{bmatrix},\end{equation}(3.4)\begin{equation}\underset{\lambda=-\mathrm{i}\kappa_j}{\mathrm{Res}}\mathbf{N}(\lambda)=\lim_{\lambda\to-\mathrm{i}\kappa_j} \mathbf{N}(\lambda)\begin{bmatrix}0&-\mathrm{i}\chi_j\mathrm{e}^{-2\mathrm{i}(\mathrm{i}\kappa_j x-4\mathrm{i}\kappa_j^3 t)}\\0& 0 \end{bmatrix},\end{equation}with norming constants
$\chi_j\in\mathbb{R_+}\setminus\{0\}$.• Normalization:
$\mathbf{N}(\lambda)=\begin{bmatrix} 1 & 1\end{bmatrix} + \mathrm{O}(\lambda^{-1})$ as
$\lambda\to\infty$.• Symmetry condition:
$\mathbf{N}(-\lambda)=\mathbf{N}(\lambda) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$.
As with Riemann--Hilbert Problem 1, we rotate the λ-plane and consider the related function 
$\mathbf{Y}(z):=\mathbf{N}(\mathrm{i} z)$. The jump contours and pole locations 
$z=\pm\kappa_j$ for 
$\mathbf{Y}(z)$ are then on the real axis. The solution 
$u(x,t)$ of the KdV equation (1.1) is obtained from (1.13). In Figure 8, we present the computed large-time evolution of a soliton gas supported on two pairs of bands, nonlinearly superposed with three solitons. In the notation of (1.5), we choose 
$f_1(z)$ to be much larger than 
$f_2(z)$. We also note here that our method is by no means limited to constant functions 
$f_j(z)$ in (1.5) (see Figure 11).
Figure 8. A KdV soliton gas with 
$r_1(\lambda)$ supported on two pairs of bands. In the notation of (1.5), 
$I_1=(1,2)$, 
$I_2=(2.5,3)$ with 
$f_1(z)=100$, 
$f_2(z)=1$ and 
$\alpha_{j}=\beta_{j}=\frac{1}{2}$, 
$j=1,2$. The three solitons superposed are associated with 
$\kappa_1=0.8$, 
$\kappa_2=2.25$ and 
$\kappa_3=3.5$ and the norming constants 
$\chi_1=10^6$, 
$\chi_2=10^5$, 
$\chi_3=10^{-12}$. Bottom panel: The same solution plotted along the ray 
$x/t = -32$.
Figure 9. A KdV soliton gas with 
$r_1(\lambda)$ supported on one pair of bands. In the notation of (1.5), 
$I_1=(1.5,2.5)$ with 
$f_1(z)=1$ and 
$\alpha_{1}=\beta_{1}=\frac{1}{2}$. The superposed soliton is associated with 
$\kappa_1=3$ and the norming constant 
$\chi_1=10^{-4}$.
Figure 10. A KdV soliton gas with 
$r_1(\lambda)$ supported on one pair of bands. In the notation of (1.5), 
$I_1=(1.5,2.5)$ with 
$f_1(z)=1$ and 
$\alpha_{1}=\beta_{1}=\frac{1}{2}$. The superposed soliton is associated with 
$\kappa_1=1$ and the norming constant 
$\chi_1=10$.
Figure 11. A KdV soliton gas with 
$r_1(\lambda)$ supported on five pairs of bands, nonlinearly superposed with five solitons. In the notation of (1.5), 
$I_1=(0.25,0.5)$, 
$I_2=(0.8,1.2)$, 
$I_3=(1.5,2)$, 
$I_4=(2.5,3)$ and 
$I_5=(4,5)$ with 
$f_1(z)=(z-0.375)^2+1$, 
$f_2(z)=(z-1)^4+1$, 
$f_3(z)=(z-1.75)^6+1$, 
$f_4(z)=\exp(z-2.75)+1$, 
$f_5(z)=\exp(-(z-4.5)^2)+1$ and 
$\alpha_1=\beta_1=-\frac{1}{2}$, 
$\alpha_2=\beta_2=\frac{1}{2}$, 
$\alpha_3=\frac{1}{2}$, 
$\beta_3=-\frac{1}{2}$, 
$\alpha_{4}=-\frac{1}{2}$, 
$\beta_4=\frac{1}{2}$, 
$\alpha_5=\beta_5=-\frac{1}{2}$. The five solitons are associated with the eigenvalue parameters 
$\kappa_1=0.1$, 
$\kappa_2=0.7$, 
$\kappa_3=2.25$, 
$\kappa_4=3.5$, 
$\kappa_5=5.5$ and the norming constants 
$\chi_1=10^5$, 
$\chi_2=1000$, 
$\chi_3=100$, 
$\chi_4=10$, 
$\chi_5=10^{-6}$.
3.1. The solution procedure
An overview of our procedure is as follows: We first ignore the poles in Riemann--Hilbert Problem 2 (in the rotated plane of 
$z=-\mathrm{i} \lambda$) and compute the matrix solution 
$\mathbf{P}(z)$, normalized so that 
$\mathbf{P}(z)\to \mathbb{I}$ as 
$z\to\infty$, and use 
$\mathbf{P}(z)$ as a global parametrix for the RHP with poles, reducing it to a discrete finite-dimensional linear algebra problem. We note that a matrix solution ceases to exist at a set of isolated points 
$(x, t)$, but we ignore this behaviour for our numerical purposes just as in [Reference Bilman and Trogdon5] and [Reference Trogdon, Olver and Deconinck13]. We solve the discrete problem, computing the solution 
$\mathbf{S}(z)$, and superpose its solution via 
$\mathbf{S}(z) \mathbf{P}(z)$ to compute the nonlinear superposition of solitons with soliton gases. We now describe the steps in this procedure in more detail.
To handle potential exponential growth in the residue conditions, we flip their triangularities as needed via the transformation
\begin{equation*}
{\hat{\mathbf{Y}}}(z)={\mathbf{Y}}(z)v(z)^{\sigma_3},\quad v(z)=\prod_{j\in K}\frac{z-\kappa_j}{z+\kappa_j},
\end{equation*}where K is an index set corresponding to ‘large’ residue conditions. In particular, for some constant c > 0, we define
\begin{equation*}
K = K(x,t) =\left\{j:\left|\chi_j\mathrm{e}^{8\kappa_j^3t-2\kappa_j x}\right| \gt c\right\},
\end{equation*}in the quiescent region and
\begin{equation*}
K=\left\{j:\left|\chi_j\mathrm{e}^{8\kappa_j^3t-2\kappa_j x+2g(z)+2h(z)}\right| \gt c\right\},
\end{equation*}in the unmodulated region. Our implementation takes c = 10 as the default value.
This transformation has the effect of modifying the residue conditions to
\begin{align*}&\underset{z=\kappa_j}{\operatorname{Res}}\,\hat{\mathbf{Y}}(z)=\lim_{z\to\kappa_j}\hat{\mathbf{Y}}(z)\begin{bmatrix}0&0 \\
\chi_j\mathrm{e}^{8\kappa_j^3t-2\kappa_j x}v(\kappa_j)^2 & 0
\end{bmatrix},\quad j\notin K,\\
&\underset{z=\kappa_j}{\operatorname{Res}}\,\hat{\mathbf{Y}}(z)=\lim_{z\to\kappa_j}\hat{\mathbf{Y}}(z)\begin{bmatrix}0 & \frac{\mathrm{e}^{-8\kappa_j^3t+2\kappa_j x}}{\chi_jv'(\kappa_j)^2}\\
0 & 0
\end{bmatrix},\quad j\in K.\end{align*}The residue conditions on the negative real axis follow from the symmetry condition.
In the unmodulated region, the disjoint disks Dj are assumed to be sufficiently small so as not to contain the poles κj. The steepest descent deformations then result in the residue conditions
\begin{align*}
&\underset{z=\kappa_j}{\operatorname{Res}}\ \mathbf{R}(z)=\lim_{z\rightarrow\kappa_j}\mathbf{R}(z)\begin{bmatrix}0 & 0 \\ \chi_j\mathrm e^{8\kappa_j^3t-2\kappa_jx+2g(\kappa_j)+2h(\kappa_j)}v(\kappa_j)^2 & 0\end{bmatrix},\quad j\not\in K,\\
&\underset{z=\kappa_j}{\operatorname{Res}}\ \mathbf{R}(z)=\lim_{z\rightarrow\kappa_j}\mathbf{R}(z)\begin{bmatrix}0 & \frac{\mathrm e^{-8\kappa_j^3t+2\kappa_jx-2g(\kappa_j)-2h(\kappa_j)}}{\chi_jv'(\kappa_j)^2}\\0 & 0 \end{bmatrix},\quad j\in K.
\end{align*} Once the matrix solution 
$\mathbf{P}$ is obtained, the residue conditions that 
$\mathbf{S}(z)=\mathbf{R}(z)\mathbf{P}(z)^{-1}$ satisfies are computed through the following formula: Given a residue condition
\begin{equation*}\underset{z=\pm{\kappa_j}}{\operatorname{Res}}\,{\mathbf{R}}(z)=\lim_{z\to\pm\kappa_j}{\mathbf{R}}(z){\boldsymbol{\Xi}_{\pm j}},\end{equation*} of 
$\mathbf{R}$, the corresponding residue condition of 
$\mathbf{S}$ is given by
\begin{equation*}\underset{z=\pm\kappa_j}{\operatorname{Res}}\,\mathbf{S}(z)=\lim_{z\to \pm\kappa_j}\mathbf{S}(z)\mathbf{P}(\pm\kappa_j){\boldsymbol{\Xi}_{\pm j}}\mathbf{P}(\pm\kappa_j)^{-1}\left(\mathbb{I}-\mathbf{P}^{\prime}(\pm\kappa_j){\boldsymbol{\Xi}_{\pm j}}\mathbf{P}(\pm\kappa_j)^{-1}\right)^{-1},\end{equation*} where 
$\mathbf{P}^{\prime}(z)$ denotes the componentwise derivative of 
$\mathbf{P}(z)$. We compute 
$\mathbf{P}^{\prime}(\pm\kappa_j)$ by numerically expanding 
$\mathbf{P}(z)$ in a Laurent series in a small circle centered at 
$\pm\kappa_j$.
4. Examples, demonstration of convergence and performance
4.1. Examples
In this subsection, we enumerate the examples found throughout this paper.
4.1.1. Five pairs of bands with five solitons
In Figures 1 and 11, we plot a soliton gas corresponding to 
$r_1(\lambda)$ supported on five pairs of bands, nonlinearly superposed with five solitons. In the notation of (1.5), Figure 11 considers variable fj and various endpoint behaviours 
$\alpha_j,\beta_j$.
4.1.2. Five pairs of bands
In Figure 2, we plot a pure soliton gas corresponding to 
$r_1(\lambda)$ supported on five pairs of bands, including for large negative x at t = 100.
4.1.3. One pair of bands with two solitons
In Figure 3, we consider a soliton gas corresponding to 
$r_1(\lambda)$ supported on one pair of bands, nonlinearly superposed with two solitons. We include a density plot of this solution in the (x, t)-plane outside of the wedge-shaped region where our method begins to break down.
4.1.4. Two pairs of bands
In Figure 4, we consider a pure soliton gas corresponding to 
$r_1(\lambda)$ supported on two pairs of bands. We plot its tail at t = 100 and its behaviour along the ray 
$x/t=-32$.
4.1.5. Two pairs of bands with three solitons
In Figure 8, we include the same plots as Figure 4, but with the solution nonlinearly superposed with three solitons.
4.1.6. One pair of bands with a taller soliton
In Figure 9, we plot the time evolution of a soliton gas corresponding to 
$r_1(\lambda)$ supported on one pair of bands, nonlinearly superposed with one taller soliton, showing the soliton escape the soliton gas.@@
4.1.7. One pair of bands with a shorter soliton
In Figure 10, we plot the time evolution of a soliton gas corresponding to 
$r_1(\lambda)$ supported on one pair of bands, nonlinearly superposed with one shorter soliton, showing the soliton becoming trapped within the soliton gas.
4.2. Convergence
In Figure 12, we demonstrate the convergence of our computational method as we increase the number of collocation points used on each component of the jump contour associated with the RHP that is treated numerically. We plot the absolute pointwise error made in computing a soliton gas potential
\begin{equation}
E_\mathrm{p}(x):=|u_{\mathrm{p}}(x,0) - u_{\mathrm{true}}(x,0)|,
\end{equation} where 
$\mathrm{p}$ refers to the number of collocation points used on each interval, denoted by ‘PPI’ in Figure 12. Here, 
$u_{\mathrm{true}}(x,0)$ denotes the solution computed with a large value of 
$\mathrm{p}$. The soliton gas potential computed for Figure 12 is associated with 
$r_1(\lambda)$ supported on two pairs of bands (in the notation of (1.5)) 
$I_1=(1.2,2)$, 
$I_2=(2.5,3)$ with 
$f_1(z)=100$ and 
$f_2(z)=1$ and 
$\alpha_{j}=\beta_{j}=\frac{1}{2}$, 
$j=1,2$. This soliton gas potential is also nonlinearly superposed with two solitons associated with the eigenvalue parameters 
$\kappa_1=1$, 
$\kappa_2=4$ and the norming constants 
$\chi_1=10$, 
$\chi_2=10^{-10}$. While we observe that our method is quite accurate, some precision is lost due to the large condition number of the collocation matrix (empirically on the order of 107). Future work will aim to improve the conditioning of this linear system.
Figure 12. Pointwise errors of various numbers for collocation points. PPI (Points Per Interval) denotes the number of collocation points used per interval Ij, and 10 times that number of points are used on each circle. 50 PPI is used as the exact solution in comparisons.
Finally, in all of the plots presented in this paper, 20 collocation points are used on each interval and 120 points are used on each circle in the jump contour for the RHP treated numerically.
4.3. Notes on performance
Perhaps the most attractive feature of our method is that it is trivially parallelizable. Since x and t are only parameters in the RHP, all computations in this paper can, in principle, be done in the time it takes to solve a single RHP, i.e., a pointwise evaluation of 
$u(x,t)$. With this perspective in mind, we give rough pointwise evaluation runtimes on a standard laptop for the solution plots above in Table 1.
Table 1. Pointwise evaluation runtimes in seconds (s) or milliseconds (ms) to compute the solutions 
$u(x,t)$ presented in each figure listed in Section 4.1
Acknowledgement
This work was supported by the National Science Foundation under Grant No. DMS-2108029 (DB) and Grant No. DMS-2306438 (TT). Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. Parts of this work were completed on Hyak, UW’s high performance computing cluster.
Open access
