Rarefaction Waves for the Toda Equation via Nonlinear Steepest Descent

We apply the method of nonlinear steepest descent to compute the long-time asymptotics of the Toda lattice with steplike initial data corresponding to a rarefaction wave.

In the case when a = 1 2 , the initial value problem (1.1)-(1.3) is called rarefaction problem. We keep this name for an arbitrary a > 0 and refer to the case a = 1 2 as the classical rarefaction (CR) problem. The long-time asymptotics of the CR problem were studied rigorously by Deift et al. [5] in 1996 in the transitional region where ξ := n t ≈ 0 as t → +∞. To this end the authors applied the nonlinear steepest descent approach for vector Riemann-Hilbert (RH) problems. Using the same approach, our aim is to study the region n t ∈ (−2a + ε, −ε) ∪ (ε, 1 − ε), where ε > 0 is a sufficiently small number. Note that the regions n t ∈ (−∞, −2a − ε) and n t ∈ (1 + ε, +∞), which are called the soliton regions, can also be studied by the vector RH approach (see [15] for decaying initial data a = 1 2 , b = 0). Although the considerations for the soliton regions in the rarefaction case are more technical than in the decaying case, they are essentially the same and lead to a sum of solitons on the respective constant background. In our opinion, the classical inverse scattering transform with the analysis of the Marchenko equation provides this result easier ( [2,3,4,20]), and consequently, we will not study the soliton regions in this paper. Moreover, the transitional regions n t ≈ 1, n t ≈ 0, and n t ≈ −2a require further analysis and are also not the subject of the present paper.
For related results on the KdV equation using an ansatz based approach see [16]. For results on the corresponding shock problem we refer to [11,17,21] and the references therein.
In summary, we will show that there are four principal sectors with the following asymptotic behavior: • In the region n > t, the solution {a(n, t), b(n, t)} is asymptotically close to the constant right background solution { 1 2 , 0} plus a sum of solitons corresponding to the eigenvalues λ j < −1.
• In the region 0 < n < t, as t → ∞ we have (1.4) a(n, t) = n 2t • In the region −2at < n < 0, as t → ∞ we have • In the region n < −2at, the solution of (1.  (1.5), where ε > 0 is an arbitrary small value. Moreover, the terms O(t −1 ) are differentiable with respect to t, and the first derivatives are of order O( n t 3 ). In the two middle regions we derive a precise formula for these error terms (see Theorem 5.1 and Proposition 6.1 below).
The following picture demonstrates the expected behavior of the Toda lattice solution in the middle regions. The numerically computed solution in Fig. 1 corresponds to "pure" steplike initial data a(n, 0) = 1 2 , b(n, 0) = 0 for n ≥ 0 and a(n, 0) = 0.4, b(n, 0) = 2 for n < 0. The apparent line is due to the fact that neighboring points are very close due to the scaling. We observe that the analytically obtained asymptotics (1.4), (1.5) and the numerically computed asymptotics match well. In particular, in the analytic case the coefficient b(n, t) has a jump in the transition region n t ≈ 0 as well. An overview of the asymptotic solutions for (1.1)-(1.2) with arbitrary constant a > 0, b ∈ R can be found in [18].
To simplify considerations we assume in addition to (1.2) that the initial data decay to their backgrounds exponentially fast where ν > 0 is an arbitrary small number. This condition allows to continue the right reflection coefficient analytically to a small vicinity of the respective spectrum.

Statement of the Riemann-Hilbert problem
Let us first recall some elementary facts from scattering theory for Jacobi operators with steplike backgrounds from [7,8,9,10] (see also [19,Chapter 10] for general background). The spectrum of the Jacobi operator H(t) associated with the equation . In addition to the spectral parameter λ we will use two other parameters z and ζ, connected with λ by the Joukowsky transformation Introduce the Jost solutions ψ(z, n, t), ψ (z, n, t) of (2.1) with asymptotic behavior lim n→∞ z −n ψ(z, n, t) = 1, |z| ≤ 1; lim n→−∞ ζ n ψ (z, n, t) = 1, |ζ| ≤ 1. Denote The points z = −1 and z = 1 correspond to the edges of the spectrum of the right background Jacobi operator, and q 2 and q 1 correspond to the respective edges of the left background operator. We will call the points z j discrete spectrum. Denote T = {z : |z| = 1} and D = {z : |z| < 1}. The map z → λ is one-to-one between the closed domains D := clos(D \ [q 2 , q 1 ]) and D : We treat closure as adding to the boundary the points of the upper and lower sides along the cuts, while considering them as distinct points. Since the function ζ n ψ (z) is in fact an analytical function of ζ, it takes complex conjugated values on the sides of the cut along the interval I := [q 2 , q 1 ], which we denote as I ± i0. Note that z ∈ I − i0 corresponds to the arc ζ ∈ {|ζ| = 1, Im ζ < 0}. The Jost solution ψ (z) takes equal real values at z, z −1 ∈ T, which yields the respective properties of the scattering matrix as a function of z. This matrix consist of the right (resp., left) reflection coefficient R(z, t) (resp. R (z, t)), defined for |z| = 1 (resp., |ζ| = 1), and the transmission coefficients T (z, t) and T (z, t) defined on D. They are connected by the scattering relations T (z, t)ψ (z, n, t) = ψ(z, n, t) + R(z, t)ψ(z, n, t), |z| = 1, Moreover, let z j , 1 ≤ j ≤ N , (note |z j | < 1) be the eigenvalues and set The set of right scattering data defines the solution of the Toda lattice uniquely. Under condition (1.6), it has the following properties (we list only those relevant for the present paper, see [5,8]): • The function R(z, t) is continuous and can be continued analytically in the annulus e −ν < |z| < 1. • The right transmission coefficient T (z, t) can be restored uniquely from (2.4) for z ∈ D. It is a meromorphic function with simple poles at z j . • The function χ(z, t) is continuous for z ∈ (q 2 , q 1 ) and vanishes at The transmission coefficient has the same behavior at q i as χ(z, t). On D we define a vector-valued function m(z) = (m 1 (z, n, t), m 2 (z, n, t)), (2.5) m(z, n, t) = T (z, t)ψ (z, n, t)z n , ψ(z, n, t)z −n .
Moreover, m 2 (z) is a meromorphic function in {z : |z| > 1} \ I * with poles at z −1 j . By definition, the vector function m(z), z ∈ C, has jumps along the unit circle and along the intervals I and I * . The statement of the respective Riemann-Hilbert problem with pole conditions is given in [5]. We will not formulate this problem here, but instead give an equivalent statement which is valid in the domain ξ := n t ∈ (−2a, 0) ∪ (0, 1) we are interested in. In this domain we can reformulate the initial meromorphic RH problem as a holomorphic RH problem as in [5,15]. We skip the details and only provide a brief outline below.
Throughout this paper, m + (z) (resp. m − (z)) will denote the limit of m(p) as p → z from the positive (resp. negative) side of an oriented contour Σ. Here the positive (resp. negative) side is the one which lies to the left (resp. right) as one traverses the contour in the direction of its orientation. Using this notation implicitly assumes that these limits exist in the sense that m(z) extends to a continuous function on the boundary. Moreover, all contours are symmetric with respect to the map z → z −1 , i.e., they contain with each point z also z −1 . The orientation on these contours should be chosen in such a way that the following symmetry is preserved for the jump matrix of the vector RH problem and for its solution.
Symmetry condition. LetΣ be a symmetric oriented contour. Then the jump matrix of the vector problem m Moreover, Most of our transformations are conjugations with diagonal matrices, so it is convenient to use the following Lemma 2.2 (Conjugation, [14]). Let m be the solution on C of the RH problem m + (z) = m − (z)v(z), z ∈Σ, which satisfies the symmetry condition as above. Let d : C \ Σ → C be a sectionally analytic function. Set then the jump matrix of the problemm + =m −ṽ is given bỹ If d satisfies d(z −1 ) = d(z) −1 for z ∈ C \ Σ, then the transformation (2.9) respects the symmetry condition.
Recall that the behavior of the solution of the RH problem is determined mostly by the behavior of the phase function Let ξ ∈ (0, 1). Part of the eigenvalues lie in the domain Re Φ(z) > 0 (namely z j ∈ (−1, 0)), while the remaining eigenvalues belong to the set Re Φ(z) < 0. The pole conditions at the eigenvalues are given by ( [15]) The pole conditions can be replaced by jump conditions on small curves around the eigenvalues as in [14]. Let be the Blaschke product corresponding to z j in the domain Re Φ(z) > 0 (if any). Note that P (z) satisfies P (z −1 ) = P −1 (z). Let δ be sufficiently small such that the circles T j = {z : |z − z j | = δ} around the eigenvalues do not intersect and lie away from T ∪ I (the precise value of δ will be chosen later). Set We consider the circles T j as contours with counterclockwise orientation. Denote their images under the map z → z −1 by T * j and orient them clockwise. These curves are not circles, but they surround z −1 j with minimal distance from the curve to z j given by δ zj (zj −δ) . We redefine the vector m by Then m ini (z) is a holomorphic function in C \ {T ∪ I ∪ I * ∪ T δ }, and solves the jump problem Note that the matrix B(z) satisfies the symmetry condition (2.7) and Here the matrix norm is to be understood as the maximum of the absolute value of its elements.
Consider the contour Γ = T ∪ I ∪ I * ∪ T δ , where the unit circle T is oriented counterclockwise and the intervals I, I * are oriented towards the center of the circle. Continue the function (2.3) to I * by χ(z) = −χ(z −1 ). Then the following proposition is valid (cf. [5,15]).
is the unique solution of the following vector Riemann-Hilbert problem: Find a vector-valued function m ini (z) which is holomorphic away from Γ, continuous up to the boundary, and satisfies: The phase function Φ(z) = Φ(z, n/t) is given by (2.10), the matrix B(z) by (2.12), the function P (z) by (2.11), and χ(z) by (2.3). II. The symmetry condition m ini (z −1 ) = m ini (z)σ 1 . III. The normalization condition Moreover, in Proposition 2.3 we assume that the points b − 2a and b + 2a are non-resonant. This means that the initial vector function m has continuous limits on I, I * and that χ(z) is bounded there (otherwise T (z) ∼ (z − q j ) −1/2 and both m and the jump matrix have singularities). For ξ ∈ (0, 1), (ii). The normalization condition m ini 1 > 0 holds since it holds for the initial function m and by definition, P (0) > 0.
We omit the proof of Proposition 2.3 which is essentially the same as in [5,15]. Uniqueness of the solution can be proven as in [1].

Reduction to the model problem
In this section we perform three conjugation-deformation steps which reduce the RH problem I-III to a simple jump problem on an arc of the circle with a constant jump matrix, plus jump matrices which are small with respect to t. The jump problem with the constant matrix can be solved explicitly. Note that since |R(z)| = 1 for z ∈ T, we cannot apply the standard lower-upper factorization of the jump matrices on an arc of T and the subsequent "lens" machinery near this arc (see [15]). For this reason we first have to find a suitable g-function ([6]).
Step 1. In this step we replace the phase function by a function with "better" properties. The g-function has the same asymptotics (up to a constant term) as Φ(z) for z → 0 and z → ∞ and the same oddness property, g(z −1 ) = −g(z). In addition, it has the convenient property that the curves separating the domains with different signs of Re g(z) cross at z 0 (ξ) ∈ T and z 0 (ξ), where g(z 0 ) = g(z 0 ) = 0. A second helpful property is that the g-function has a jump along the arc connecting z 0 and z 0 which satisfies g + (z) = −g − (z) > 0. This simplifies further transformations, because with this property and Lemma 2.2 we do not need the lens machinery around this arc. Note that z 0 does not coincide with the stationary phase point of Φ. Recall that for ξ ∈ (0, 1), the curves Re Φ(z) = 0 intersect at the symmetric points −ξ ± ξ 2 − 1 ∈ T, which are the stationary points of Φ. That is, the stationary phase point corresponds to the angle φ 0 ∈ (0, π) where cos φ 0 = −ξ. Set and introduce where √ s > 0 for s > 0. The cut of the square root in (3.2) is taken between the points z 0 and z 0 along the arc We orient Σ in the same way as T, i.e., from z 0 to z 0 .
Lemma 3.1. The function g(z) defined in (3.2) satisfies the following properties: Proof. We first prove that our choice of z 0 yields Φ(z) − g(z) = O(1) as z → ∞. To match the asymptotics of g(z) and Φ(z) in (2.10) we compute which implies that cos θ 0 = 1 − 2ξ as desired. To prove property (b), substitute s = e iθ and z 0 = e iθ0 in (3.2), then Now it is straightforward to see (c), because the integrand can be represented as For property (d), note that due to their equal asymptotic behavior, the signature table for g as z → 0 or z → ∞ is the same as the signature table for Φ (see Fig. 2). The line T \ Σ corresponds to Re g = 0. Indeed, if z 1 = e iθ1 and Re z 1 > Re z 0 , then The function g(z) has a jump along the contour Σ, but the limiting values are real (compare with the proof of (b)). Thus Re g + = g + > 0 (because this limit is taken from the domain where Re g > 0). Respectively, Re g − = g − < 0 and g + = −g − . Now we are ready to finish the proof of property (a). Since g(z) satisfies (c) we have g(1) = 0. By use of (3.3) and (2.2), and taking into account that we obtain for z → +0, that is when √ λ 2 − 1 > 0 as λ > 1, the asymptotic behavior Differentiating (3.5) we also get To prove (e) we decompose (3.2) in a vicinity of z 0 and integrate. Since arg(z − z 0 ) = (− 3π 2 + θ 0 )(1 + o(1)) for z ∈ Σ in a small vicinity of z 0 , then arg(z − z 0 ) 3/2 + π 4 − 3θ0 2 = −2π(1 + o(1)), which implies the second claim of (e). The signature table for g(z) is given in Fig. 2. With this description of the g- where the following notations have been introduced: (3.7) The matrix B(z) was defined in (2.12). Note that in the non-resonant case for z = q 1 and z = q 2 , To obtain an analogous estimate on T δ , we have to adjust the value for δ. Denote inf j inf{|Φ(z j )|, |g(z j )|} = J > 0, j = 1, . . . , N.
Choose δ > 0 so small that Then sup Step 2: On T \ Σ, the jump matrix can be factorized using the standard upperlower factorization ( [5,15]). Let C be a contour close to the complementary arc T\Σ with endpoints z 0 and z 0 and clockwise orientation. Let C * be its image under the map z → z −1 , oriented clockwise as well. We denote the regions adjacent to      Figure 3. Contour deformation of Step 2.
these contours by Ω and Ω * as in Fig. 3, and set Redefine m (1) inside Ω and Ω * by m (2) (z) = m (1) (z)W (z). Then the new vector does not have a jump along T \ Σ, and satisfies m The symmetry and normalization conditions are preserved in this deformation step.
Before we perform the next conjugation step, we have to study in detail the solution of the following scalar conjugation problem: find a holomorphic functioñ d(z) on C \ Σ, such that Remark 3.2. As for any multiplicative scalar jump problem with non-vanishing jump function on a contour in C, one can find its solution via the associated additive jump problem and the usual Cauchy integral. However, the representation via the usual Cauchy integral requires a prescribed (and known) behavior of the solutions at ∞. On the other hand, this representation cannot provide (3.9), unless the jump functions are not even on Σ. For example, condition (3.9), (i), implies logd(1) = 0, which cannot be obtained by the usual Cauchy integral. That is why we use the Cauchy integral with kernel vanishing at z = 1, In order to solve the conjugation problem (3.8)-(3.9), we first solve an auxiliary conjugation problem: find a holomorphic function F (z) in C\Σ, bounded as z → ∞ and such that Since there are two associated additive jump problems, log F + (z) = log F − (z) ± iπ, z ∈ Σ, and since iπ 2πi , the solution of the jump problem in (3.11) can be given by The symmetry F (z −1 ) = −F (z) is evident from here. It turns out that this symmetry implies that F + (s) is an even function for s ∈ Σ. Moreover, F + (s) is real-valued as s = e iθ and θ = θ 0 , (3.12) because cos θ < cos θ 0 on Σ. On the other hand, and in particular, F (1) = 0 and F (∞) = −F (0) = cos θ0 2 ∈ R. Now assume thatd(z) solves (3.8)-(3.9) and introduce the function f (z) := F (z) logd(z). Then (3.8), (3.9), and (3.11) imply that this function should solve the jump problem (3.14) f . Moreover, f (z −1 ) = f (z) should hold for z ∈ C \ Σ and f (∞) ∈ R. The jump (3.14) can be satisfied by (3.15) f One has to check that the other two conditions are satisfied too. The function log R(z) R(−1) is an odd function. Indeed, recall that the reflection coefficient R(z) is a continuous function which satisfies R(z −1 ) = R −1 (z), |R(z)| = 1, and so does P 2 (z), which implies R(z −1 ) = R −1 (z). Since F + (z) is an even function, h(z) is odd on Σ, and the required evenness of f can be directly verified from (3.15).
To check that f (∞) ∈ R recall that |R(z)| = 1 yields h(z) ∈ iR. When z → ∞ the first integrand in (3.15) vanishes due to oddness of h. For the second we obtain Since 1 / ∈ Σ, the function f (z) is holomorphic in a vicinity of z = 1 and has a zero at least of first order there. Hence F −1 (z)f (z) is well defined, it satisfies the required additive jump problem for logd(z), the symmetry Thus,d(z) = e F −1 (z)f (z) solves the conjugation problem (3.8)-(3.9) and combining (3.10), (3.12)-(3.15) we obtain (3.16) Moreover, taking into account that , and using the first equality in (3.13) we obtain another representation ford(z), .
Remark 3.4. The function R, which is in fact a function of the spectral parameter λ, depends on θ via R(cos θ), therefore the derivative of u(θ) is a real valued function.
Proof. The function u(θ) is an odd function of θ on the interval of integration in the sense that u( π 2 − α) = −u( π 2 + α). In the same sense, sin θ 2 (cos θ 0 − cos θ) −1/2 is an even function. Taking this into account as well as (3.10), we have for z → ∞ Integration by parts then yields plus the term of order z −2 . On the other hand, Combining this with (3.21) yields (3.19) and (3.20).
Since r(s) − r(z) ∼ (s − z), the integral in J 1 (z) is Hölder continuous in a vicinity of z 0 . Therefore, On the other hand, since and (q(z, z 0 )) −1 → 0 as z → ∞, we have , and (3.22) follows from (3.23) in a straightforward manner.
Theorem 3.6. The vector function m (3) (z) is the unique solution of the following RH problem: find a holomorphic vector functionm(z) in C \ (Σ ∪ C ∪ C * ∪ Ξ), which is continuous up to the boundary and has the following properties: • It solves the jump problemm Here E(z) is given by (3.7), (2.12),d(z) by (3.16), andȆ(z) satisfies where J is defined by (3.8); •m(z) satisfies symmetry and normalization conditions as (2.8) and (2.13); moreover,v(z) has the symmetry property (2.7); • For small z, the vector function m in (2.5) andm(z) are connected by   The proof of this Lemma is analogous to the uniqueness proof of the vector model problem in [1].

The model problem solution
For our further investigation we will also need a matrix solution M mod (z) = M mod (z, n, t) of the matrix RHP: find a holomorphic matrix function M mod (z) on C \ Σ satisfying the following jump and symmetry conditions, M mod We find a solution of the matrix problem following [13,1]. Consider the nonresonant case, that is, we first look for a holomorphic solution of the jump problem M 0 + = iM 0 − σ 3 satisfying the symmetry condition as above. Using (3.10) we get where the branch of the fourth root is chosen with a cut along the negative half axis and 1 1/4 = 1. Since β(z −1 ) = β −1 (z), we have the required symmetry. Concerning the original matrix solution M mod (z), this begs for the representation Thus, In the resonant case the matrix solution is represented by the same formula (4.2), but with β −1 (z) instead of β(z) and vice versa. In summary, we have the following  (z) as z → 0 is given by in the non-resonant case. In the resonant case, Proof. Consider first the non-resonant case. Since since cos θ 0 = 1−2ξ. Respectively, sin θ 0 = 2 ξ − ξ 2 , from which (4.4) follows. The resonant case follows from the same computation using β −1 (z) instead of β(z).
The structure of the matrix solution (4.2) and the jump matrix (3.24) as well as the results of Lemmas 3.5, 3.1, (e), allow us to conclude that the solution of the parametrix problem (which has a local character) can be constructed as in [1] using [15,Appendix B] (in particular, the solution can be given in terms of Airy functions). Consequently, as z → 0 (cf. [12]) m (3) (z, ξ) =m(z) = m mod (z, ξ) + F (ξ) t + z H(ξ) t + O(z 2 ) + G(z, ξ, t) t 4/3 , where the vector functions F (ξ) and H(ξ) are Hölder continuous as ξ ∈ I := [ε, 1−ε] with an exponent α ≥ 1/2. The vector-function G(z, ξ, t) is uniformly bounded with respect to z ∈ D δ , ξ ∈ I and t ∈ [T, ∞), where T > 0 is a sufficiently large number and D δ is a small circle centered at 0. By (3.25), we have (5.1) Here the term O(t −4/3 ) is uniformly bounded with respect to z ∈ D δ and ξ ∈ I, and the term O(z 2 ) is uniformly bounded with respect to t ∈ [T, ∞) and ξ ∈ I. 6. Discussion of the region −2at < n < 0.
To obtain the asymptotic of the solution in the region n t ∈ (−2a, 0) we could study an analogous RH problem connected with the left scattering data and the variable ζ (cf. (2.2)). Instead of this extensive analysis let us consider the Toda lattice associated with the functions