Matrix Boussinesq solitons and their tropical limit

We study soliton solutions of matrix"good"Boussinesq equations, generated via a binary Darboux transformation. Essential features of these solutions are revealed via their"tropical limit", as exploited in previous work about the KP equation. This limit associates a point particle interaction picture with a soliton (wave) solution.


Introduction
The (scalar) Boussinesq equation originated from the study of water waves propagating in a canal [1]. In this work we study the "good" Boussinesq equation, in which the highest derivative term has the opposite sign, as compared to the original Boussinesq equation. It also appears as a continuum limit of certain nonlinear atomic chains, see [2], for example. From a mathematical point of view, it belongs to the main examples of completely integrable PDEs. It is the second Gelfand-Dickey reduction [3] of the KP-II equation, the famous Korteweg-deVries (KdV) equation being the first. Compared with the latter, the behavior of its soliton solutions is considerably more diverse, in particular they can move in both directions, experience head-on collision, and solitons can split or merge (cf. [4,5], also see the references cited in these papers). The scalar good Boussinesq equation is also known as "nonlinear string equation", see e.g. [6,7].
In the KdV and Boussinesq case, a contour plot of a soliton solution displays a localization along a piecewise linear graph in two-dimensional space-time. The definition of the "tropical limit" makes this precise and yields a method to compute it. Restricting the dependent variable to this graph, completes the tropical limit of the soliton solution, which displays the essentials of soliton interactions in a clear way. It is a very convenient tool to describe and to classify soliton solutions. The tropical limit associates with a soliton solution a point particle picture, also in the interaction region of solitons, revealing "virtual solitons".
In this work we address more generally the following m×n matrix potential Boussinesq equation, where K is a constant n × m matrix, β > 0 (chosen as β = 1/4 in all plots in this work), and a subscript indicates a partial derivative of φ with respect to the respective variable x or t. We will refer to this equation as potential Bsq K . It is the second Gelfand-Dickey reduction of the following KP-II equation, in a moving frame (x → x − 3βt 3 , t 2 = t), In terms of u = 2φ x , we obtain the m × n matrix Boussinesq equation, or rather Bsq K , 1 For a solution u of the vector Boussinesq equation, where n = 1 or m = 1, Ku, respectively uK, is a solution of the scalar Boussinesq equation. If m, n > 1, we define the tropical limit graph via that of the scalar tr(Ku), which is not in general a solution of the scalar Boussinesq equation. Our explorations and results fully substantiate this approach. Our analysis largely parallels that in [8], where we explored line soliton solutions of the matrix KP-II equation in the tropical limit. There we concentrated on a class of solutions which we called "pure solitons". Another class has been treated in [9]. Whereas pure solitons exhaust the solitons of the KdV reduction, this is not so for the Boussinesq reduction.
Matrix versions of scalar integrable equations are natural extensions and of interest as models of coupled systems. A further motivation to explore them originated from the fact that in 2-soliton scattering, in particular in case of the matrix KdV [10,11] and the vector Nonlinear Schrödinger (NLS) equation [12,13], in-and outgoing polarizations (values of the dependent variable attached to in-and outgoing solitons) are related by a Yang-Baxter map. For a matrix Boussinesq equation, a Yang-Baxter map is not sufficient to fully describe the situation.
In Section 2 we present a binary Darboux transformation for the potential Bsq K equation (1.1). Appendix A explains its origin from a general result in bidifferential calculus [14,15]. In this work we concentrate on the case of vanishing seed solution. This leads to two cases, treated in Sections 3 and 4. The soliton solutions, obtained via the binary Darboux transformation, depend on parameters that have to be roots of a cubic equation. We introduce a convenient parametrization of these roots (see Section 3.1) that greatly facilitates the further analysis. Section 5 contains some concluding remarks.

A binary Darboux transformation for the matrix Boussinesq equation
The following binary Darboux transformation is a special case of a general result in bidifferential calculus, see Appendix A. Let N ∈ N. The integrability condition of the linear system where θ is an m × N and C a constant N × N matrix, is the potential Bsq K equation for φ 0 . The same holds for the adjoint linear system where χ is an N × n and C a constant N × N matrix, in general different from C. The Darboux potential Ω satisfies the consistent N × N system of equations is then a new solution of the potential Bsq K equation.
Remark 2.1. Taking the transpose of the above equations, and applying the substitution t → −t, we see that, besides K → K T , we also obtain C ↔ C T and θ ↔ χ T . We also note that with any invertible constant N × N matrices A and B.
Using (2.4) and the first of (2.3), we find Such a formula is familiar in the scalar case, where det Ω is the Hirota τ -function. But we will see that, also in the matrix case, det Ω plays a crucial role. In the following we will still call it τ , after multiplication by a convenient factor, which preserves the relation (2.5).

Solutions for vanishing seed
The linear system with φ 0 = 0 reads It possesses solutions of the form where ϑ(P ) = P x + P 2 t , (2.6) and each P a is a solution of the cubic equation The index a runs over any number of distinct roots. Correspondingly, the adjoint linear system takes the form which is solved by if Q b solves the cubic equation The equations for the Darboux potential Ω are reduced to Writing it follows that W ba has to satisfy the Sylvester equation 10) and the constant matrix Ω 0 is subject to As a consequence of the last condition, there are two major cases. In Section 3 we will address the case where C = −C. Section 4 then deals with the complement. We note that (2.5) reduces to tr(Ku) = 2 (log det Ω) xx . (2.12) 3 The case C = −C If Ω 0 is invertible, Remark 2.2 shows that without restriction of generality we can choose Ω 0 = I N , the N × N identity matrix. (2.11) then implies C = −C. The remaining freedom of transformations, according to Remark 2.2, is then given by transformations with B = A −1 . The similarity transformation C → A −1 CA now allows us to assume that C has Jordan normal form. P a and Q b are now solutions of the same cubic equation, so we can set Q a = P a , and the Sylvester equation takes the form If a = b and P a and P b have disjoint spectrum, it is well-known that there is a solution and it is unique. In this case, the sum in the expression for θ or χ is over a disjoint set of solutions of the cubic equation. We will restrict our considerations to diagonal matrices 2 P a = diag(p 1,a , . . . , p N,a ) , C = diag(c 1 , . . . , c N ) , 2 A treatment of the case where C contains larger than size 1 Jordan blocks is left aside in this work.
we have We will only consider real roots. Writing and W ab = (W ab,ij ), we find and thus

A parametrization of the roots of the cubic equation
We are only interested in real soliton solutions, hence we need real roots of the cubic equation, and at least two different ones. This requires |c i | ≤ 2β 3/2 . We can then express the constants c i as follows, where λ i are real parameters. The roots of the cubic equation (3.1) are then given by All the roots satisfy p 2 ≤ 4β. We have lim |λ i |→∞ p ia = √ β for a = 1, 2, and lim |λ i |→∞ p i,3 = −2 √ β. Also see Fig. 1. The two involutive transformations generate the permutation group of the three roots.

Pure solitons
This is the case with only a single root of each of the cubic equations (2.7) and (2.8). These equations coincide in the case under consideration, so we have to choose two different roots, P 1 = P and P 2 = Q. We will assume that these matrices are diagonal, hence P = diag(p 1 , . . . , p N ) =: diag(p 1,1 , . . . , p N,1 ) , Q = diag(q 1 , . . . , q N ) =: diag(p 1,2 , . . . , p N,2 ) . It will be convenient to allow both notations for the (diagonal) entries. The constants p i , q i have to solve the cubic equation z 3 = 3βz + c i . We will require q i = p i , i = 1, . . . , N . Moreover, we will mostly also assume p i = p j and q i = q j for i = j. Writing we then have This is exactly the expression for Ω that we found in [8] for the KP K equation. The only difference is that, for i = 1, . . . , N , q i = p i now have to satisfy the cubic equations, so that p 3 i −3βp i = q 3 i −3βq i holds, which is But, apart from this, all formulae derived in [8] for the potential KP K equation also apply to the case under consideration. Next we summarize those that are needed in this work. Introducing where adj(Ω) is the adjugate of the matrix Ω, we find that they have expansions where µ I are constants, M I constant matrices, and Recall that p k,1 = p k and p k,2 = q k , k = 1, . . . , N . We have Since τ I,x = p I τ I , p I = p 1,a 1 + · · · + p N,a N , where I = (a 1 , . . . , a N ), we obtain Note that (cf. [8]) in accordance with (2.12).
If τ = 0 and if µ I > 0 for all I with µ I = 0 in the expression for τ , regularity of the solution is guaranteed. Let We call this the region where ϑ I dominates. As intersection of half-spaces, it is convex.
If µ I > 0, the tropical limit of φ in U I is 3 On such a (visible) segment the value of u is given by We find tr(Ku IJ ) = 1 2 (p I − p J ) 2 and introduce normalized valueŝ which satisfy tr(Kû IJ ) = 1. Using the notation the k-th soliton appears in space-time on segments of the straight lines determined by log τ I k (1) = log τ I k (2) , i.e., This also determines the asymptotic structure of a tropical limit graph of a pure N -soliton solution.
Without restriction of generality, we can order the parameters such that p 1 + q 1 < p 2 + q 2 < · · · < p N + q N . 4 If we represent p i and q i by the first two roots in (3.4), then p i + q i is a strictly increasing function of |λ|, hence this order is obtained by choosing 0 < λ 1 < λ 2 < · · · < λ N . Now it follows from (3.7) that, for t 0, the solitons appear along the x-axis according to their numbering, and for t 0 they appear in reverse order. To the left of the center is the dominating phase region U 1,...,1 . Then follows counterclockwise U 2,1,...,1 , U 2,2,1,...,1 , ..., unless we get to the region U 2,...,2 on the right hand side. Correspondingly, starting again from U 1,...,1 , we get clockwise to the regions U 1,...,2 , U 1,...,2,2 , ..., until we arrive at U 2,...,2 . These regions always appear in a tropical limit graph of a pure N -soliton solution. The remaining U I can only appear as bounded regions. But some may be empty. This depends on the value of higher Boussinesq hierarchy variables, see Fig. 5 below.
Remark 3.1. It can happen that µ I = 0 for some multi-index I, so that e ϑ I is absent in the expression for τ , but that this exponential appears in the numerator of the expression for u, i.e., M I = 0. Then some components of u will exhibit exponential growth in the region where this phase dominates, and a tropical limit does not exist.

Single soliton solution
For the 1-soliton solution we find assuming ηKξ > 0, and This yields Setting it takes the form Via the symmetries (3.5), 1-soliton solutions with other choices of the roots are obtained from the above solution. If λ 2 < 1/3, the soliton moves from left to right. If λ 2 > 1/3, it moves from right to left. For λ 2 = 1/3, it is stationary. In all cases, the absolute value of the velocity is less than 2 √ β. We also note that 0 ≤ tr(Ku) ≤ 6β. The tropical limit graph of the 1-soliton solution is the boundary between the two dominating phase regions U 1 and U 2 . It is the straight line in space-time (xt-plane), determined by with slope −1/(p + q). We have

2-soliton solution
In this case (N = 2), we find Hence, if α, κ 1,1 , κ 2,2 = 0, Choosing we obtain via φ = F/τ a 2-soliton solution of the potential Bsq K equation. 5 Fig. 2 shows examples of corresponding tropical limit graphs. Applying the symmetries (3.5), 2-soliton solutions with other choices of the roots are obtained from the above solution. Furthermore, we find The map R : (û 11,21 ,û 21,22 ) → (û 12,22 ,û 11,12 ) is a nonlinear Yang-Baxter map, it satisfies the (quantum) Yang-Baxter equation  [17,7]), here denoted as a "composite" 12, which finally splits. This has no counterpart in the KdV case. In the other cases, the numbering of the virtual soliton draws on a formal analogy with particles and anti-particles (the latter indicated by a bar over the respective number).
This map is also obtained from that of the KP theory, see Section 5 in [8], with the parameters restricted by (3.9). The fact that R satisfies the Yang-Baxter equation follows from the Yang-Baxter property of the more general KP K Yang-Baxter map. Alternatively, it follows from the 3-soliton solution and the fact that the polarizations do not depend on the independent variables, including higher hierarchy variables, see Section 3.2.5.
In the vector case, we obtain a linear map, where Remark 3.2. Dropping the exponential factor in (3.6), which has only been introduced to achieve a convenient numbering of phases, in the N = 2 case we obtaiñ Comparison with a known expression for the τ -function of the 2-soliton solution of the scalar Boussinesq (or KP) equation shows that this determines a solution of the scalar Boussinesq equation if κ 1,1 κ 2,2 = κ 1,2 κ 2,1 . We also note that, if κ 1,2 κ 2,1 = 0, the above expression factorizes toτ = (1 + e ζ 1 )(1 + e ζ 2 ). In this case, the tropical limit graph is simply the superimposition of those of the factors 6 , hence there is no phase shift.
The corresponding tropical limit graphs are Y-shaped (a soliton splits into two), respectively reverse Y-shaped (two solitons merge), see Fig. 3. Approaching such a solution by letting q 2 → q 1 , respectively p 2 → p 1 , in the 2-soliton solution in Section 3.2.2, we see that the edge representing the virtual soliton (cf. the second and third graph in Fig. 2) gets longer and longer, in such a way that the dominating phase region U 1,1 finally disappears at infinity.
In these cases the R-matrix reduces to
I ij (2,2) i i j j i j Figure 4: A 2-soliton part of an N -soliton tropical limit graph.
Also see, e.g., [18] (Section 5 therein) for such approximations in the case of KP solitons. The two parallel line segments corresponding to the path of the i-th soliton are determined by Their shift along the x-axis, caused by the interaction with the j-th soliton, is The parallel line segments of the j-th soliton are given by Their shift along the x-axis, caused by the interaction with the i-th soliton, is If N = 2, we have A I 1,2 = α/(κ 1,1 κ 2,2 ) (also see [17] for the scalar case).
The first and the third graph correspond to a large negative, respectively large positive value of s. The sequences of 2-soliton interactions are according to the left, respectively right hand side of the Yang-Baxter equation. Since the polarizations along edges of a tropical limit graph do not depend on the variables (x, t, s), we conclude that the Yang-Baxter equation holds.

Other soliton configurations
Now we consider solutions involving three roots of the cubic equation. There are then two cases, either θ = θ 1 e ϑ(P 1 ) + θ 2 e ϑ(P 2 ) , The second, "dual" choice can be obtained from the first by applying the symmetry φ → φ T , t → −t, and using χ a ↔ θ T a , P a → −P T a . In the following we use again the decomposition (3.2), but with the rescaling ξ ia → (p i,3 − p ia ) ξ ia , a = 1, 2. and We set The solution describes the merging of two solitons into a single one, also see Fig. 6. If λ ∈ {0, ±1/3, ±1}, two of the p i are equal and the solution reduces to a 1-soliton solution. The dual case describes the splitting of a single soliton into two.
Remark 3.4. In contrast to the KdV reduction of KP, where a quadratic equation rules the game, Bsq K thus admits solutions with a tropical limit graph in space-time having the most elementary rooted binary tree shape: three edges meeting at a vertex. We may speculate that in higher Gelfand-Dickey reductions there are corresponding limits for the number of edges forming a rooted binary tree.  Example 3.6. According to Proposition 3.7 below, in the scalar and vector case there is no regular solution in the class given by Example 3.5, with all possible phases present in the expression for τ . But corresponding regular solutions exist, for example, in the 2 × 2 matrix case. Fig. 7 shows two examples of tropical limit graphs. Here we chose K = diag(1, 1) and λ 1 = 7/10, λ 2 = 4. Proof. In the vector case, where η i are scalars, we have κ 2,1,a κ 1,2,b = κ 1,1,a κ 2,2,b , hence
The last proposition tells us that, however small (with respect to a suitable norm) we choose a neighborhood of such a regular solution, in the set of solutions, it contains singular solutions.
An N > 2 solution locally consists approximately of N = 2 solutions. We should thus expect the above proposition to extend to N > 2. But we will not attempt to provide a rigorous proof here. In the following we show that solutions with N = 3 exist, where the singularities only appear in a compact region of space-time. N = 3. We set P a = diag(p a , q a , r a ), a = 1, 2, 3, and obtain τ = 3 a,b,c=1 µ abc eθ abc (3.13) withθ abc = ϑ(p a ) + ϑ(q a ) + ϑ(r a ) and where δ i a is the Kronecker symbol. Furthermore, If τ has also negative summands, strictly its tropical limit is not defined. Notwithstanding this, we can determine (and plot) the regions where the logarithm of the absolute value of a summand of τ dominates. On a boundary segment between such (generalized) dominating phase regions, where the corresponding summand in τ is positive for one and negative for the other, the soliton solution is singular. This is so because close to such a boundary segment contributions from all other summands are exponentially suppressed, hence negligible. Furthermore, singularities can only appear at such a boundary segment, since everywhere else either a single summand of τ dominates and all others are negligible, or we have a boundary segment along which two summands with the same sign have equal values and all other summands are negligible. The latter two cases exclude a singularity. Fig. 9 shows an example of such a modified tropical limit graph. The white regions are dominated by phases having negative contributions to the τ -function. In these regions the solution is still regular, but on the interface between a "positive" and a "negative" phase region (plotted in red), and only there, the solution becomes singular. A similar solution of the scalar Boussinesq equation appeared in [5] (see Figure 8 therein).

The case C = −C
We can use the transformations in Remark 2.2 to achieve that both, C and C , have Jordan normal form. Then (2.11) generically implies Ω 0 = 0. This will be assumed in the following. Choosing diagonal matrices P a = diag (p 1,a , . . . , p N,a ) , Q a = diag(q 1,a , . . . , q N,a ) , C = diag(c 1 , . . . , c N ) , C = diag(c 1 , . . . , c N ) , (2.7) and (2.8) require For each i = 1, . . . , N , we represent the roots p ia , a = 1, 2, 3, of the first cubic equation as in (3.4), using a parameter λ i . In the same way we represent the roots q ia , a = 1, 2, 3, of the second cubic equation using a parameter ν i . We assume that p ia = p jb and q ia = q jb for i = j or a = b. Now we have θ = Assuming q ib = p ja for all combinations of indices, from (2.9) and (2.10) we obtain η ib Kξ ja q ib − p ja e ϑ(p ja )−ϑ(q ib ) i, j = 1, . . . , N .
The Bsq K solution is given by ϑ IJ . If a phase ϑ IJ is present in τ , and if the tropical limit exists, then the tropical value of φ in the corresponding dominating phase region is φ IJ . We decompose θ a and χ a as in (3.2). In the following we restrict our considerations to the case N = 1. Then Ω is a scalar, consisting of at most nine summands. We obtain φ = F/τ with τ = Ω = Hence φ ab = p a − q b κ ba θ a ⊗ χ b .
In the scalar and vector Boussinesq case, one can show that there is no regular solution of the above form, with real parameters and with all coefficients of exponentials in τ different from zero. This means that, however small (using a suitable norm) we choose a neighborhood of such a regular solution, in the set of solutions, it contains singular solutions. The first statement is not true in the matrix case. and λ 1 = 1/5, ν 1 = 3/2. The function τ , given above, is then a positive linear combination of nine independent exponentials, which is the maximal number. All nine phases, appearing in the corresponding function τ , are visible in the tropical limit graph, see Fig. 10. The three interior phase regions are hardly expected from the plot on the left hand side and reveal a complicated interaction pattern. Fig. 11 shows another, though simpler example, where the relation is evident.  Figure 11: Plot of tr(Ku) = 2(log τ ) xx and tropical limit graph of a 3 × 3 matrix solution, as given in Section 4, with N = 1 and θ 1 = (−1, 1, 1) T , θ 2 = (−2, 1, 1), θ 3 = (−1, −1, −1), λ 1 = 1/10, ν 1 = 1/10 + 10 −5 . Only seven of the nine phases, appearing in the corresponding function τ , are visible. the matrix Boussinesq equation We recall a binary Darboux transformation result of bidifferential calculus [14,15]. Let Ω solve the compatible linear system Γ Ω − Ω ∆ = −η K θ , dΩ = (dΩ) ∆ − (dΓ) Ω + (dη) K θ + κ Ω + Ω λ .
In the above theorem, we have to assume that all objects are such that the corresponding products are defined and that d andd can be applied. Next we define a bidifferential calculus via on the algebra A = A 0 [∂ t , ∂ x ], where A 0 is the algebra of smooth functions of two variables, x and t, and ∂ x is the operator of partial differentiation with respect to x. ζ 1 , ζ 2 are a basis of a two-dimensional vector space V , from which we form the Grassmann algebra Λ(V ). d andd extend to Ω = A ⊗ Λ(V ) in a canonical way, and to matrices with entries in Ω.