Long-Time Asymptotics of a Three-Component Coupled mKdV System

: We present an application of the nonlinear steepest descent method to a three-component coupled mKdV system associated with a 4 × 4 matrix spectral problem. An integrable coupled mKdV hierarchy with three potentials is ﬁrst generated. Based on the corresponding oscillatory Riemann-Hilbert problem, the leading asympototics of the three-component mKdV system is then evaluated by using the nonlinear steepest descent method.


Introduction
Asymptotic analysis has become an important area in soliton theory and its basic approaches are used to determine limiting behaviors of Cauchy problems of integrable systems. Significant work on long-time asympototics of integrable systems was carried out by Manakov [1] and Ablowitz and Newell [2]. Zakharov and Manakov [3] wrote down precise expressions, depending explicitly on initial data, for the leading asymptotics for the nonlinear Schrödinger equation in the physically interesting region x = O(t). A complete description was presented for the leading asymptotics of the Cauchy problem of the Korteweg-de Vries (KdV) equation by Ablowitz and Segur [4]. The method of Zakharov and Manakov starts an ansatz for the asymptotic form of the solution and utilizes some techniques which are removed from the classical framework of Riemann-Hilbert (RH) problems. Segur and Ablowitz [5] began with the similarity solution form to derive the leading two terms in each of the asymptotic expansions for the amplitude and phase for the nonlinear Schrödinger equation, based on conservation laws. However, all those results were presented, without applying the method of stationary phase.
Its [6] used the stationary phase idea to conjugate the RH problem associated with the nonlinear Schrödinger equation, up to small errors which decay as t → ∞, by an appropriate parametrix to a model RH problem solvable by the technique from the theory of isomonodromic deformations. Deift and Zhou [7] determined the long-time asymptotics of the mKdV equation, by manipulating an associated oscillatory RH problem systematically and rigorously, in the spirit of the stationary phase method. Their technique, further developed in [8,9] and also in [10], opens a nonlinear steepest descent method to explore asymptotics of integrable systems through analyzing oscillatory RH problems associated with matrix spectral problems. A crucial ingredient of the Deift-Zhou approach is the asymptotic analysis of singular integrals on contours by deformations. Applications have been made for a few of integrable equations, including the KdV equation, the nonlinear Schrödinger equation, the sine-Gordon equation, the derivative nonlinear Schrödinger equation, the Camassa-Holm equation and the Kundu-Eckhaus equation (see, e.g., [11][12][13][14][15][16][17][18]). McLaughlin and Miller [19] generalized the steepest descent method to the case when the jump matrix fails to be analytic and their method is now called a nonlinear∂ steepest descent method. One important factor in the nonlinear steepest descent method is the order of involved spectral matrices in RH problems or equivalently the inverse scattering theory. However, only 2 × 2 spectral matrices and their RH problems have been systematically considered (see, e.g., [20]), which lead to algebro-geometric solutions to integrable systems expressed by hyperelliptic functions [21]. There are very few 3 × 3 spectral matrices, whose long-time asymptotics or RH problems are considered (see, e.g., [22,23]) and whose associated inverse scattering transforms are solved (see, e.g., [24,25]). Associated trigonal curves exhibit much more diverse asymptotic behaviors and algebro-geometric solutions than hyperelliptic curves [26,27]. There has been no application of the nonlinear steepest descent method to the 4th-order or higher-order matrix spectral problems so far.
In this paper, we would like to present an application to a 4 × 4 matrix spectral problem where i denotes the unit imaginary number, k is the spectral parameter and the superscript * denotes the complex conjugate. This spectral problem generates a three-component coupled mKdV system p j,t = −p j,xxx + 3(|p 1 | 2 + |p 2 | 2 + |p 3 | 2 )p j,x + 3(p * 1 p 1,x + p * 2 p 2,x + p * 3 p 3,x )p j , 1 ≤ j ≤ 3, where |p i | 2 = p i p * i , 1 ≤ j ≤ 3.
We will compute the leading asymptotics of this mKdV system by analyzing an associated oscillatory RH problem. Symmetry constraints and RH problems have been made for an unreduced combined mKdV system in [28,29]. It is worth noting that this mKdV system (2) has a slightly different matrix spectral problem from the one for the multiple wave interaction equations [30]. Let |M| denote an equivalent matrix norm for a matrix M (may not be square): where M † stands for the Hermitian transpose of M, and where S denotes the Schwartz space. The primary result of the paper can be stated as follows.
Theorem 1. Let p(x, t) = (p 1 (x, t), p 2 (x, t), p 3 (x, t)) solves the Cauchy problem of the three-component coupled mKdV system (2) with initial data in S 3 . Suppose that γ(k) = (γ 1 (k), γ 2 (k), γ 3 (k)) is the reflection coefficient vector in S 3 associated with the initial data and it satisfies Then, in the physically interesting region x = O(t), the leading asymptotics of the solution for x < 0 is given by where Γ(·) being the Gamma function.
The rest of the paper is organized as follows. In Section 2, within the zero-curvature formulation, we derive an integrable coupled hierarchy with three potentials and furnish its bi-Hamiltonian structure, based on the 4 × 4 matrix spectral problem (1) suited for the RH theory. In Section 3, taking the three-component coupled mKdV system (2) as an example, we present an associated oscillatory RH problem. In Section 4, we explore long-time asymptotics for the three-component coupled mKdV system through manipulating the oscillatory RH problem by the nonlinear steepest descent method. In the last section, we give a summary of our conclusions, together with some discussions.

Zero Curvature Formulation
Let us first recall the zero curvature formulation to construct integrable hierarchies [31]. Let u be a vector potential, k be a spectral parameter, and I n stand for the n-th order identity matrix. Choose a square spectral matrix U = U(u, k) from a given matrix loop algebra. Suppose that presents a solution to the corresponding stationary zero curvature equation By using the solution W, we define an infinite sequence of Lax matrices where the subscript + denotes the operation of taking a polynomial part in k, and ∆ r , r ≥ 0, are appropriate modification terms, and then construct an integrable hierarchy from an infinite sequence of zero curvature equations The two matrices U and V [r] here are called a Lax pair [32] of the r-th integrable system in the soliton hierarchy (9). The zero curvature equations in (10) present the compatibility requirements of the spatial and temporal matrix spectral problems with φ being the matrix eigenfunction.
To show the Liouville integrability of the soliton hierarchy (9), we usually furnish a bi-Hamiltonian structure [33]: where J 1 and J 2 constitute a Hamiltonian pair and δ δu is the variational derivative. The Hamiltonian structures are normally achieved through the trace identity [31]: or more generally, the variational identity [34]: with ·, · being a symmetric and ad-invariant non-degenerate bilinear form on the underlying matrix loop algebra [35]. The bi-Hamiltonian structure often guarantees the existence of infinitely many commuting Lie symmetries {K n } ∞ n=0 and conserved quantities {H n } ∞ n=0 : and where n 1 , n 2 ≥ 0, N = J 1 or J 2 , and K denotes the Gateaux derivative of K with respect to the potential u: Such Abelian algebras of symmetries and conserved quantities can be generated directly from Lax pairs (see, e.g., [36,37]). We also know that for a system of evolution equations, is a conserved functional if and only if δH δu is an adjoint symmetry [38], and thus, Hamiltonian structures connect conserved functionals with adjoint symmetries and further symmetries. Pairs of adjoint symmetries and symmetries actually correspond to conservation laws [39].
When the underlying matrix loop algebra in the zero curvature formulation is simple, semisimple and non-semisimple, the associated zero curvature equations generate classical integrable hierarchies [40,41], a collection of different integrable hierarchies, and hierarchies of integrable couplings [42], respectively. Integrable couplings require extra care in presenting lumps and solitons.

Three-Component mKdV Hierarchy
Let us consider a 4 × 4 matrix spectral problem where k is a spectral parameter and u = p T = (p 1 , p 2 , p 3 ) T is a three-component potential. A special case of p 3 = 0 transforms (17) into the Manakov type spectral problem [43]. Since the leading matrix diag(−1, 1, 1, 1) has a multiple eigenvalue, the spectral problem (17) is degenerate, which is different from the case of multiple wave interaction equations [30].
To derive an associated integrable coupled mKdV hierarchy, we first solve the stationary zero curvature Equation (7) corresponding to (17). We assume that a solution W is determined by where a is a real scalar, b is a three-dimensional row, and d is a 3 × 3 Hermitian matrix. It is direct to observe that the stationary zero curvature Equation (7) becomes We search for a formal series solution as follows: with b [m] and d [m] being assumed to be Obviously, the system (19) is equivalent to the following recursion relations: The three-component mKdV system (2) can correspond to the specific initial values: Besides those initial values, we set all constants of integration in (22c) to be zero, that is, demand This way, with a [0] and d [0] given by (23), all matrices W m , m ≥ 1, will be uniquely determined. For instance, by a direct computation, (22) generates b [1] j = 4p j , a [1] = 0, d [1] where 1 ≤ j, l ≤ 3. Based on (22b) and (22c), we can obtain a recursion relation for b [m] and where Ψ is a 6 × 6 matrix operator with q = −p † . Obviously, this tells To generate an integrable coupled mKdV hierarchy with three components, we introduce, for all integers r ≥ 0, the following Lax matrices where the modification terms are chosen as zero. The compatibility requirements of (11), i.e., the zero curvature equations (10), lead to the integrable coupled mKdV hierarchy with three components: The first two nonlinear systems in the above integrable coupled mKdV hierarchy (30) are the three-component coupled nonlinear Schrödinger system and the three-component coupled mKdV system (2). Under a reduction p 3 = 0, the three-component coupled system (31) is reduced to the Manakov system [43], for which a decomposition into finite-dimensional integrable Hamiltonian systems was made in [44], whileas the three-component coupled system (2) contains various examples of mKdV equations, for which there are diverse kinds of integrable decompositions through symmetry constraints (see, e.g., [45,46]). The three-component coupled mKdV hierarchy (30) with an extended six-component potential u = (p, q T ) T = (p, −p * ) T possesses a Hamiltonian structure [29,38], which can be computed through applying the trace identity [31], or more generally, the variational identity [34]. A direct computation tells Plugging these into the corresponding trace identity and considering the case of m = 2 tell γ = 0, and thus A bi-Hamiltonian structure for the extended six-component coupled mKdV systems, consisting of (30) and its conjugate compartment, then follows: where the Hamiltonian pair (J 1 , J 2 = J 1 Ψ) is defined by where q = −p † again. Adjoint symmetry constraints (or equivalently symmetry constraints) decompose a four-component nonlinear Schrödinger system with p 3 = q 3 = 0 into two commuting finite-dimensional Liouville integrable Hamiltonian systems in [38]. In the next section, we will concentrate on the three-component coupled mKdV system (2).

An Associated Oscillatory Riemann-Hilbert Problem
For a matrix M = M(k; x, t), we define where |M| is the matrix norm of M given by (3) and f p is the L p -norm of a function f of k. We assume that as t → ∞, If the constant C depends on a few parameters α 1 , α 2 , · · · , α n , then we write A(t) α 1 ,α 2 ,··· ,α n B(t), but on some situations where the dependence of C on some parameters is not important, we will often suppress those parameters for brevity.
For an oriented contour, we denote the left-hand side by + and the right-hand side by −, as one travels on the contour in the direction of the arrow. A Riemann-Hilbert (RH) problem (Γ, J) on an oriented contour Γ ⊂ C (open or closed) with a jump matrix J defined on Γ reads where J 0 is a given matrix determining a boundary condition at infinity, and M ± are analytic in the ± side regions and continuous to Γ from the ± sides, and M = M ± in the ± side regions, respectively.

An Equivalent Matrix Spectral Problem
In the previous section, we have seen that the spectral problems of the three-component coupled mKdV system (2) with the Lax pair being of the form where u = p T = (p 1 , p 2 , p 3 ) T and The other two matrices P and Q are given by where (25), and thus As normal, for the spectral problems in (38), upon making the variable transformation we can impose the canonical normalization condition: The equivalent pair of matrix spectral problems reads by a generalized Liouville's formula [47].
Applying the method of variation in parameters, and using the canonical normalization condition (44), we can transform the x-part of (38) into the following Volterra integral equations for ψ ± [48]: Similarly, we can turn the t-part of (38) into the following Volterra integral equations: Observing the structures of Λ and Ω, we can know that the first column of ψ − and the last three columns of ψ + consist of analytical functions in the upper half plane C + , and the first column of ψ + and the last three columns of ψ − consist of analytical functions in the lower half plane C − (also see [29,49]). All this will help us formulate an associated RH problem for the three-component coupled mKdV system (2).

An Oscillatory Riemann-Hilbert Problem
The scattering matrix S is determined by where e αM X = e αM Xe −αM for a scalar α and two same order square matrices M and X. We also adopt the simple expressions we have and Notice that det S(λ) = 1 due to det ψ ± = 1. It then follows from det S = 1 and (55) that where adj(M) is the adjoint matrix of M, and S is the block matrix with S 22 being a 3 × 3 matrix. Thus, the scattering matrix S can be written as where a is a 3 × 3 matrix. Introduce where ψ ±,L and ψ ±,R denote the first column and the rest columns of ψ ± . This matrix M solves an associated oscillatory RH problem where with In the above RH problem, we assume that the matrix a(k) is invertible, and so we will only analyze the solitonless case. The behavior of the oscillatory factor e 2itθ(k) depends on the increasing and decreasing of θ(k) and the signature of Re iθ(k) (see Figures 1 and 2).  In what follows, we assume that γ lies in S 3 and satisfies The analysis on RH problems in [50] shows that the above RH problem (58) is equivalent to a Fredholm integral equation of the second kind. For such Fredholm equations, the existences and uniqueness of solutions can be guaranteed by the vanishing lemma [51]. A direct computation also shows [29] that one can compute the potential p(x, t) through the solution M(k; x, t) to the RH problem (58) as follows.
Theorem 2. Assume that γ lies in S 3 and satisfies the conditions in (61). Then there exists a unique solution M(k; x, t) to the Riemann-Hilbert problem (58), and the solution of the three-component coupled mKdV system (2) is recovered via where we partition M into a block matrix M = (M jl ) 2×2 with M 22 being a 3 × 3 matrix.

Long-Time Asymptotic Behavior
We will first deal with the Riemann-Hilbert (RH) problem (58) and then compute the long-time aymptotics of the three-component coupled mKdV system (2) with the leading term, by using the nonlinear steepest descent method [7]. We will be concentrated on the physically interesting region | x t | ≤ C, where C is a positive constant.

Transformation of the RH Problem
Obviously, the jump matrix J has the following upper-lower and lower-upper factorizations: and We would like to introduce another RH problem to replace the upper-lower factorization with the lower-upper factorization to make the analytic and decay properties of the two factorizations consistent.
We easily see that the stationary points of θ are ±k 0 , where which implies a scalar RH problem Due to (61), the jump matrix Thus, it follows from the vanishing lemma (see, e.g., [51]) that the RH problem (65) has a unique solution δ(k). In addition, the Plemelj formula [51] gives the unique solution of the above scalar RH problem (66): where By the uniqueness, we obtain It then follows that and Therefore, by the maximum principle for analytic functions, we can have from the canonical normalization condition of the above two RH problems.
Let us now introduce and a vector-valued spectral induced function and reverse the orientation for |k| > k 0 as shown in Figure 3. Then, we see that M ∆ presents the solution of the RH problem on R oriented as in Figure 3: where the jump matrix is defined by We will deform this RH problem to evaluate the asymptotics of the three-component coupled mKdV system (2).

Decomposition of the Spectral Induced Function
To determine the required deformation, we first make a decomposition of the spectral induced function ρ(k). By L, we denote the contour and by Σ, denote the contour with the orientation in Figure 4. Further let L ε denote the following part of the contour L: where 0 < ε < √ 2. We will focus on the contour Σ, though any other contour of the same shape as Σ, which locates in the the region where Re iθ(k) is positive, will work. Proposition 1. The vector-valued spectral induced function ρ(k) has the following decomposition on the real axis: where R(k) is a polynomial on |k| < k 0 and a rational function on |k| ≥ k 0 , h 2 (k) has an analytic continuation from R to L in the region where Re iθ(k) > 0, and h 1 , h 2 and R have the estimates as t → ∞: where l is an arbitrary positive integer. Taking the Hermitian conjugate yields the same estimates for e 2itθ(k) h † 1 (k * ), e 2itθ(k) h † 2 (k * ) and e 2itθ(k) R † (k * ) on R, L * and L * ε , respectively.
Proof. Throughout the proof, we use the differential notationds = 1 √ 2π ds for brevity, while dealing with the Fourier transform. We only consider the physically interesting region x = O(t) and so k 0 is bounded. Let r be a fixed positive integer.
(a) First, we consider the case of |k| < k 0 . In this case, we have ρ(k) = −γ(k)(1 − γ(k)γ † (k * )) −1 . Splitting it into even and odd parts leads to where H e and H o are two vector functions in S 3 . By Taylor's theorem with the integral form of the remainder, we have and Set Then, we have and the coefficients decay rapidly as k 0 → ∞. We assume that r is of form with an even number q. Express By the characteristic of R in (89), we have Based on this property, we will try to split h into two parts where h 1 is small and h 2 has an analytic continuation from R to L in the region where Re iθ(k) > 0. This way, we obtain the required splitting of ρ. Let us introduce We consider the Fourier transform with respect to θ. Because k → θ(k) is one-to-one in |k| < k 0 (see Figure 1), we define Based on (94), we have and as dk/dθ = [12(k 2 (θ) − k 2 0 )] −1 , we see that where the H j s are Hilbert spaces. It follows from the Fourier inversion theorem that where h/α is the Fourier transform By the formulae (86), (87), (88) and (93) we have Then, it follows that from which, by the Plancherel theorem, we know Now make a splitting for h as follows: On one hand, based on (104) and noting |k| < k 0 , we have On the other hand, we know that Re iθ(k) is positive in the hatched region in Figure 5, and thus that h 2 (k) has an analytic continuation from R to the line segments in the upper half k-plane: Let k be on the line segment (107). Then, we can have that again based on (104). But Therefore, we have Similarly, we can show that the same estimate holds on the line segment (108). Finally fix 0 < ε < √ 2. Then on the parts of the line segments (107) and (108) with ε < w ≤ √ 2 away from ±k 0 , we have (by, e.g., (110) on the line segment (107) and R(k) is a polynomial), where τ = tk 3 0 . (b) Second, we consider the case of |k| ≥ k 0 . We only consider the sub-case k ≥ k 0 , since the other sub-case k ≤ −k 0 is completely similar. In this case, we have Once more, we use Taylor's theorem to obtain Define and set h(k) = ρ(k) − R(k).
As before, we know that and decay rapidly as k 0 → ∞. Similarly, introduce Following the Fourier inversion theorem, we have where h/β is the Fourier transform Based on the formulae (114), (116) and (119), we see that where from which it follows that d j g(k, k 0 ) dk j 1, k ≥ k 0 , j ≥ 0.
Noting that we have By the Plancherel theorem, Again, we split On one hand, for k ≥ k 0 , similarly as in the previous example (106), we see that On the other hand, h 2 (k) has an analytic continuation from R to the ray in two parts of the lower half k-plane, where Re iθ(k) > 0, as shown in Figure 6. Let k be on the ray (130). Similarly as in the previous case |k| < k 0 , we can show that However, we have Re iθ(k) ≥ 4k 3 0 w 2 ( and thus, we have since k 0 is bounded. This completes the proof of Proposition 1. We are now ready to make a deformation of the RH problem (76).

Deformation of the RH Problem
Let us first explain what a deformation of a RH problem is. Suppose that we have a RH problem (Γ, J) on an oriented contour Γ (see Figure 7): and that on a part (which could be the whole contour Γ) of Γ from k 1 to k 2 in the direction of Γ, denoted by Γ k 1 k 2 , we have where b ± have invertible and analytic continuations to the ± sides of a region D (see Figure 7) supported by k 1 and k 2 , respectively. Define an extended oriented contour Γ (see Figure 7): and an extended jump matrix J : Then the RH problem (133) on Γ is equivalent to the following RH problem on Γ : and the relation between the two solutions reads It is clear that one RH problem is solvable if and only if the other RH problem is solvable, and the solution to the one problem gives the solution to the other problem automatically by (138). We call (Γ , J ) a deformation of the original RH problem (Γ, J). If D is not bounded, then we need to require that b + (k) and b − (k) tend to the identity matrix as k → ∞, in order to keep the same normalization condition. We now deform the RH problem (76) from R to the augmented contour Σ. The first important step is to observe that the jump matrix J ∆ (k; , x, t) can be rewritten as where Moreover, noting the decomposition (81), we have and thus, the decompositions for b ± : We can also state that It finally follows that the jump matrix J ∆ has the following factorization: It is now standard to introduce and then the RH problem (76) on R becomes the following new RH problem on Σ: where The canonical normalization condition in (148) can be verified, indeed. For example, (b a + ) −1 converges to I 4 as k → ∞ in Ω 6 ∪ Ω 8 , because we observe that for fixed x, t, by the definition of h 2 in (128) and the boundedness of det δ and δ in (72), we have and by the definition of R in (115) and the boundedness of det δ and δ in (72), we have both of which converge to 0 as k → ∞ in Ω 6 ∪ Ω 8 .
The above RH problem (148) can be solved by using the Cauchy operators as follows (see [52,53]). Let us denote by C ± the two Cauchy operators on Σ. As is well known, these two operators C ± are bounded from L 2 (Σ) to L 2 (Σ), and satisfy for a 4 × 4 matrix-valued function f , where Assume that µ = µ (k; x, t) ∈ L 2 (Σ) + L ∞ (Σ) solves the fundamental inverse equation Then, we can see that presents the unique solution of the RH problem (148). (2) is expressed by

Theorem 3. The solution p(x, t) of the three-component coupled mKdV system
Proof. This statement can be shown by Theorem 2, (75) and (147).

Contour Truncation
Let Σ be the truncated contour Σ = Σ \ (R ∪ L ε ∪ L * ε ) with the orientation as in Figure 8. We will reduce the RH problem (148) from Σ to Σ , and estimate the difference between the two RH problems, one on Σ and the other on Σ . Let ω e : Σ → M(4, C) be a sum of three terms where ω a = ω R is supported on R and is composed of the contributions to ω from terms of type h 1 (k) and h † 1 (k * ), ω b = ω L ∪ L * is supported on L ∪ L * and is composed of terms of type h 2 (k) and h † 2 (k * ), and ω c = ω L ε ∪ L * ε is supported on L ε ∪ L * ε and is composed of terms of type R(k) and R † (k * ). Define Obviously ω = 0 on Σ \ Σ , and hence ω is supported on Σ from terms of type R(k) and R † (k * ).
Proof. The estimates (158), (159) and (160) can be derived, similarly to Proposition 1. Based on the early definition of R(k), we can see that On this contour, by (110), we have Re iθ(k) ≥ 8k 3 0 α 2 , and then using (72), we find that It then follows from a direct computation that the estimate (161) holds.
Similarly to Proposition 2.23 in [7], we can show the following estimate.

Proposition 2.
In the physically interesting region x = O(t), as τ → ∞, there exists the inverse of the operator 1 − C ω : L 2 (Σ) → L 2 (Σ), which is uniformly bounded: Corollary 1. In the physically interesting region x = O(t), as τ → ∞, there exists the inverse of the operator 1 − C ω : L 2 (Σ) → L 2 (Σ), which is uniformly bounded: Proof. Noting C ω = C ω + C ω e , we have In addition, it follows from Lemma 1 that in the physically interesting region x = O(t) (where k 0 is bounded and thus t = O(τ)). Therefore, we have Second, the second resolvent identity tells Again based on (164), we see from this identity that the estimate in the corollary is just a consequence of Proposition 2.

Theorem 4.
As τ → ∞, we have Proof. First, from (165), we can obtain Second, it follows directly from Lemma 1 that This proves the theorem, together with (167).
Note that as k ∈ Σ \ Σ , ω (k) = 0, we can reduce C ω from L 2 (Σ) to L 2 (Σ ), and for simplicity's sake, we still denote this reduced operator by C ω . Then and thus it follows from Theorems 3 and 4 that the following statement holds.
Let L = L \ L ε , which yields Σ = L ∪ L * , and denote µ = (1 − C ω ) −1 I 4 . Then, we see that solves the following RH problem where the jump matrix is defined by

Disconnecting Contour Components
Denote the contour Σ = Σ A ∪ Σ B (see Figure 8), and set where Further define and thus The two Cauchy operators C ω A and C ω B (see (150) for definition) are bounded operators from L 2 (Σ ) → L 2 (Σ ).
One can prove a similar lemma to Lemma 3.5 in [7].

Lemma 2.
We have the estimates Proof. From the identity In the above computation, we denote for simplicity's sake. Note that if an operator T is bounded and the inverse of 1 − T exists, then (1 − T) −1 is also bounded. Further, based on the two Lemmas, Lemmas 1 and 2, and one proposition, Proposition 2, we can derive the estimate (179). Now, from Theorem 5 and Proposition 3, we see that the following result holds.

Rescaling and Reduction of the Disconnected RH Problem
Extend the two contours Σ A and Σ B to the following two contours respectively, whereω Let Σ A and Σ B denote the contours {k = k 0 αe ± πi 4 : α ∈ R} oriented outward as in Σ A ,Σ A and inward as in Σ B ,Σ B , respectively (see Figure 9). Motivated by the method of stationary phase, define the scaling transformations and denote A direct change-of-variable argument tells that where C ω A (or C ω B ) is a bounded operator from L 2 (Σ A ) (or L 2 (Σ B ) into L 2 (Σ A ) (or L 2 (Σ B )). On the part of the contour Σ A , we have and on the conjugate part L * A , we have where Lemma 3. As t → ∞, we have the estimates Proof. We only prove the estimate (193) and the proof of the other estimate is completely similar. From (65) and (66), we see thatδ solves the following RH problem where The solution of this vector RH problem is given bỹ with Noting that we can have f (k) = O((k 2 − k 2 0 ) l+1 ) when k → ±k 0 , upon noting the definition of h = h 1 + h 2 , and decompose f (k) into two parts: and f 2 (k) has an analytical continuation from R to L t (see Figure 10): and satisfies Let k ∈ L A , and we decompose where with Based on (193), (194) and Lemma 3.35 in [7], we can obtain Therefore, we have Together with a similar argument for B, we can obtain For k ∈ C \ Σ A , we set which solves the following RH problem Particularly, from an asymptotic expansion we get 2πi .
An analogous RH problem for B o on Σ B reads and its solution is given by The above jump matrix is defined by where

The Model RH Problem and Its Solution
To determine (M A o ) 12 explicitly, we solve a model RH problem The solution to this RH problem is given by Since the jump matrix v(k 0 ) does not depend on k along each ray of Σ A , we have Together with (216), this implies that dΦ(k) dk Φ −1 (k) has no jump discontinuity along any of the four rays. Directly from the solution, we obtain It then follows from Liouville's theorem that where Particularly, we have (M A o 1 ) 12 = iδ 2 A η 12 .
We partition Φ(k) into the following form where Φ 11 (k) is a scalar and Φ 22 (k) is a 3 × 3 matrix. From the differential Equation (219) for Φ, we obtain Upon plugging (228) into (227), separating the coefficients of the two independent special functions presents Finally, we conclude that (4) is a consequence of (215), (221) and (229).

Concluding Remarks
The paper is dedicated to determination of long-time asymptotics of the three-component coupled mKdV system, based on an associated Riemann-Hilbert (RH) problem. The crucial step is to deform the associated RH problem to the one which is solvable explicitly, and estimate small errors between solutions to different deformed RH problems. This is an example of applications of the nonlinear steepest descent method to asymptotics in the case of 4 × 4 matrix spectral problems.
There are various other approaches in the field of integrable systems, which include the Hirota direct method [59], the generalized bilinear technique [60], the Wronskian technique [61,62] and the Darboux transformation [63]. Connections between different approaches would be interesting. About coupled mKdV equations, there are many other studies such as integrable couplings [64], super hierarchies [65] and fractional analogous equations [66,67], and an important topic for further study is long-time asymptotics of those generalized integrable counterparts via the nonlinear steepest descent method. It is hoped that our result could be helpful in computing limiting behaviors of solutions incorporating features of other exact solutions, such as lumps [68,69], from the perspective of steepest descent based on RH problems.