Tropical limit of matrix solitons and entwining Yang-Baxter maps

We consider a matrix refactorization problem, i.e., a"Lax representation", for the Yang-Baxter map that originated as the map of polarizations from the"pure"2-soliton solution of a matrix KP equation. Using the Lax matrix and its inverse, a related refactorization problem determines another map, which is not a solution of the Yang-Baxter equation, but satisfies a mixed version of the Yang-Baxter equation together with the Yang-Baxter map. Such maps have been called"entwining Yang-Baxter maps"in recent work. In fact, the map of polarizations obtained from a pure 2-soliton solution of a matrix KP equation, and already for the matrix KdV reduction, is NOT in general a Yang-Baxter map, but it is described by one of the two maps or their inverses. We clarify why the weaker version of the Yang-Baxter equation holds, by exploring the pure 3-soliton solution in the"tropical limit", where the 3-soliton interaction decomposes into 2-soliton interactions. Here this is elaborated for pure soliton solutions, generated via a binary Darboux transformation, of matrix generalizations of the two-dimensional Toda lattice equation, where we meet the same entwining Yang-Baxter maps as in the KP case, indicating a kind of universality.


Introduction
The quantum Yang-Baxter equation is known to be a crucial structure underlying two-dimensional integrable QFT models. A typical feature of the latter is the factorization of the scattering matrix into contributions from 2-particle interactions (see [1] and references therein). This factorization is also typical for the scattering of solitons of classical nonlinear integrable field equations. Indeed, we tend to think of a multi-soliton solution of some vector or matrix version of an integrable nonlinear partial differential of difference equation as being composed of 2-soliton interactions. Matrix solitons carry "internal degrees of freedom", called "polarization". In some cases, like matrix KdV [2] or vector NLS [3,4,5], the map from incoming to outgoing polarizations has been found to satisfy the Yang-Baxter equation. But why should we expect the Yang-Baxter property? The latter is a statement about three particles, here solitons, and may be thought of as expressing independence of the different ways in which a three-particle interaction can be decomposed into 2particle interactions. First of all, how to decompose a 3-soliton solution into 2-soliton interactions? Because of the wave nature of solitons, there are no definite events at which the interaction takes place, and (with the exception of asymptotically incoming and outgoing solitons) there are no definite values of the dependent variable which could define a corresponding map. However, there is a certain limit, called "tropical limit", that takes soliton waves to "point particles" and then indeed determines events at which an interaction occurs. This has been a decisive tool in our previous work [6,7,8,9,10,11] and this will be so also in this work, which continues an exploration started Let S be the set of rank one m × n K-projection matrices (XKX = X), and R(1, 2) := R(p 1 , q 1 ; p 2 , q 2 ) : S × S → S × S (X 1 , X 2 ) → (X 1 , X 2 ) be given by where 1 m denotes the m × m identity matrix and α 12 := α(p 1 , q 1 , X 1 ; p 2 , q 2 , X 2 ) : This is a parameter-dependent Yang-Baxter map, which means that it satisfies the Yang-Baxter equation on S × S × S. The indices of R ij specify on which two of the three factors the map R acts. Here we have to assume that the constants p i , i = 1, 2, 3, and also q i , i = 1, 2, 3, are pairwise distinct, and that the expressions for α ij make sense.
If q i = −p i , this is the Yang-Baxter map obtained from the 2-soliton solution of the KdV K equation For the matrix KdV equation (where m = n and K = 1 n ), the Yang-Baxter map has first been derived in [2]. 2 In the particular case where n = 1 (and correspondingly for m = 1), the above (generically highly nonlinear) Yang-Baxter map becomes linear : (X 1 , X 2 ) = (X 1 , X 2 ) R(i, j) , R(i, j) := (1.5) 1 In this work, we only consider matrices over the real or complex numbers. 2 The factor α12 is missing in the latter work, but it is necessary for the Yang-Baxter property.
Here we used the fact that, for X ∈ S, KX is now a scalar, so that the K-projection property requires KX = 1. The R-matrix, which ermerges here, solves the Yang-Baxter equation on a threefold direct sum of an m-dimensional vector space, which extends the set S.
In the tropical limit of a pure N -soliton solution of the KP K equation, the dependent variable u has support on a piecewise linear structure in R 3 (with coordinates x, y, t), a configuration of pieces of planes, and the dependent variable takes a constant value on each plane segment. This piecewise linear structure is obtained as the boundary of "dominating phase regions". In the KdV reduction, the support of the dependent variable in the tropical limit is a piecewise linear graph in 2-dimensional space-time. For the KdV 2-soliton solution, we have four dominating phase regions, numbered by 11, 12, 21, respectively 22. 3 Fig. 1 shows an example. Using the general 2-soliton solution to compute the values of the dependent variable u along the boundary line segments of the tropical limit graph, after normalization toû such that tr(Kû) = 1, the map (û 11,21 ,û 21,22 ) → (û 12,22 ,û 11,12 ) yields the above Yang-Baxter map (with the KdV reduction q i = −p i ). Here, e.g.,û 11,21 is the polarization along the boundary line between the dominating phase regions numbered by 11 and 21. For a rank one matrix, the above normalization condition is equivalent to the K-projection property.
But what about a phase constellation different from the one shown in Fig. 1 ? Indeed, Fig. 2 displays alternatives. In the same way as for the phase constellation in Fig. 1, one finds that the  Figure 2: Other dominating phase region constellations and tropical limit graphs in 2dimensional space-time, for 2-soliton solutions of the KdV K equation. In each case, the Yang-Baxter map R is recovered if we consider the graph as defining a map of polarizations not from bottom to top, but in a different space-time direction (here from the soliton lines marked by red numbers to those marked by blue numbers).
first alternative in Fig. 2 leads to the inverse of the above KdV Yang-Baxter map. The remaining 3 The parameters p k , q k belong to the k-th soliton. The first digit of the phase number ab refers to soliton 1, the second to soliton 2. Since we write p k =: p k,1 and q k =: p k,2 , we have a, b ∈ {1, 2}. possibilities, however, determine maps that are not Yang-Baxter. Nevertheless, they are realized in matrix KdV 2-soliton interactions (see Appendix A), regarding them as a process evolving in time t and by choosing the parameters appropriately. We learn that there are matrix KdV 2-soliton solutions for which the map of polarizations in t-direction is not Yang-Baxter! 4 How to understand this, in view of our different expectation?
In all cases of phase constellations, shown in Fig. 2, the Yang-Baxter map R is present, however. It is recovered by regarding the plot not as a process in t-direction, but in a different direction in space-time. 5 First of all, this means that we overlooked something in our analysis of the KP K case in [9]. As in the KdV reduction, of course also the general pure 2-soliton solution of KP K contains constellations, for certain parameter values, where incoming and outgoing polarizations (with respect to a chosen direction) are related by a map that is not a Yang-Baxter map. Besides the above Yang-Baxter map, this map and the inverses of both maps are needed to describe the propagation of polarizations along the support of pure multi-soliton solutions in the tropical limit.
We were actually led to the new insights by exploring a matrix version of the two-dimensional Toda lattice equation (see, e.g., [13,14] for the scalar equation). This is the subject of Section 4. In particular, with the restriction to "pure solitons", it turns out that the same Yang-Baxter map is here at work as in the KP K case, indicating a kind of universality.
The tropical limit associates with a soliton solution a configuration of plane segments, together with values of the dependent variable on the segments. 6 It is found that, at intersections, these polarizations are related by one of two maps (and their inverses), of which only one is a Yang-Baxter map, but the two maps satisfy a mixed version of the Yang-Baxter equation (see (3.6) below). They are "entwining Yang-Baxter maps" in the sense of [16]. Section 2 presents a "Lax representation" for the above map R. This is a matrix refactorization problem. The basic argument 7 is the same as in [2] for the matrix KdV case (also see [18,19]), but we prove more directly, as compared with [2], that the refactorization problem determines the map R.
In Section 3, we show that this refactorization problem, written in a different way, also determines the abovementioned mixed version of the Yang-Baxter equation. It implies further relations which in particular lead to solutions of the "WXZ system" in [20], called "Yang-Baxter system" in [21]. To our knowledge, such a system first appeared in [22].
In Section 4 we explore soliton solutions of the abovementioned matrix 2-dimensional Toda lattice equation. Section 4.1 presents a binary Darboux transformation for the matrix potential two-dimensional Toda lattice equation. Its origin from a general result in bidifferential calculus is explained in Appendix C. We then concentrate on the case of vanishing seed solution. Section 4.2 further restricts to a subclass of soliton solutions, which we call "pure", and we define the tropical limit of such solitons. In Section 4.3 we derive the Yang-Baxter map R from the pure 2-soliton solution. The relevance of the aforementioned additional non-Yang-Baxter map is explained in Section 4.4, and Section 4.5, which treats the case of three pure solitons, shows explicitly how the 4 [2] uses results of [12], where a restriction has been imposed on the parameters of the matrix KdV 2-soliton solution. As a consequence of this, the non-Yang-Baxter cases are excluded in [2]. 5 As a process in t-direction, the first plot in Fig. 2 corresponds to an application of the inverse of the Yang-Baxter map R. 6 Here we think of the discrete independent variable k in the Toda lattice equation as being continuously extended (also see [15]). Such a smoothing of the discrete variable is actually done in the plots presented in Section 4 of this work. But, of course, it is not assumed in any of our computations. 7 It is actually more generally based on the relation between neighboring simplex equations, see [17] and references cited there.
Yang-Baxter map and the non-Yang-Baxter map (and their inverses) are at work, and why they have to be "entwining".
Finally, Section 5 contains some concluding remarks.

A Lax representation for the Yang-Baxter map
Let K be an n × m matrix with maximal rank, and where X is an m × n matrix and λ a parameter. Then we have If X is a K-projection matrix, which means XKX = X, then Theorem 2.1. Let p 1 , p 2 , q 1 , q 2 be pairwise distinct and X i , i = 1, 2, rank one K-projections, hence X i ∈ S. Then the refactorization equations 8 imply the map R(1, 2), defined in the introduction.
A proof is given in Appendix B. Recalling a well-known argument (see [17] and references cited there), exploiting associativity in different ways, we obtain where we abbreviated R ij (i, j) to R ij and set, for example, R(1, 3)(Y 1 , X 3 ) =: (Z 1 , Y 3 ), and also There are corresponding chains with A i replaced byÃ i . If These are local 1-simplex equations, see [17], for example. determines a unique map (X 1 , X 2 , X 3 ) → (Z 1 , Z 2 , Z 3 ), which means that 9 we can conclude the statement of the following theorem. But it can also be verified directly, using computer algebra.
Theorem 2.2. Let X i ∈ S. Then R, given by (1.1), is a Yang-Baxter map.
(2.2) is called a "Lax representation" for the map R.
Using (2.1), (1.1) can be expressed as In particular, we have which in particular means that α 12 is an invariant of the map R.

Further aspects of the Lax representation
According to Section 2, is a Lax representation for the Yang-Baxter map R. More precisely, we have to supplement this equation byÃ 1 (λ, A corresponding extension is also necessary for the other versions of (3.1) considered below, but for simplicity we will suppress it.

A Lax representation for the inverse of R is given by
As a consequence, we have Comparison with (2.3) shows that this is obtained from the latter by exchanging the two indices 1 and 2. This means that R is a reversible Yang-Baxter map,

Let us write
instead of (3.1). As in Section 2 and Appendix B, it can be shown that this equation uniquely determines the map The denominators of the final expressions are both equal to α 12 . This map is invariant under exchange of the two indices 1 and 2, hence Although (3.5) resembles (1.1), in contrast to the latter it does not yield a Yang-Baxter map. This can be checked using computer algebra. As a consequence of associativity, we have and also where we used (3.1) and set, for example, . One can argue that this implies Z i = Z i , i = 1, 2, 3, in which case we can conclude that This can also be verified using computer algebra. Hence, writing the Lax representation (3.1) in the form (3.4), we are led to a "mixed Yang-Baxter equation" for the two maps R and T .

The inverse
A corresponding Lax representation is the version of (3.1).
Remark 3.1. In the special case where n = 1, besides R also T becomes linear : It is easily verified that the matrix T does not satisfy the Yang-Baxter equation.

Further consequences of the Lax representation
There are actually further consequences of the fact that (3.1), (3.2), (3.4) and (3.8) uniquely determine maps.

The two ways to rewrite
by using (3.1) and (3.4), allow us to deduce that in the two possible ways, with certain Z i , using (3.2) and (3.4), leads to But this equivalent to (3.9).
Remark 3.2. A system of equations like that given by the Yang-Baxter equation (1.3), supplemented by (3.9) or (3.12), appeared in [23] under the name "braided Yang-Baxter equations". The same holds for the Yang-Baxter equation for R −1 , supplemented by (3.10) or (3.11). Also see [22]. Via (3.9) and (3.10), as well as via (3.11) and (3.12), we have examples of what has been called "WXZ system" in [20], later also named "Yang-Baxter system" [21]. This system apparently first appeared in [22]. Here we obtained solutions of these systems.

The p2DTL K equation
In this section, we address the following matrix version of the potential two-dimensional Toda lattice equation, where ϕ is an m × n matrix of (real or complex) functions and K a constant n × m matrix of maximal rank. A subscript indicates a partial derivative with respect to the respective variable, here x or y. A superscript + or − means a shift, respectively inverse shift, in a discrete variable, which we will denote by k. We refer to this equation as p2DTL K .
In terms of new independent variables If ϕ is independent of z, the last equation reduces to We will refer to this equation as p1DTL K .
Correspondingly, we may regard (4.3) as a matrix version of the potential one-dimensional Toda lattice equation.

A binary Darboux transformation for the p2DTL K equation
The following binary Darboux transformation is a special case of a general result in bidifferential calculus, see Appendix C. Let N ∈ N. The integrability condition of the linear system where θ is an m × N matrix, is the p2DTL K equation for ϕ 0 . The same holds for the adjoint linear system where χ is an N × n matrix. So let ϕ 0 be a given solution of (4.1). Let the Darboux potential Ω satisfy the consistent system of N × N matrix equations Where Ω is invertible, is then a new solution of the p2DTL K equation (4.1).
This observation is helpful in order to reduce the set of parameters, on which a generated solution depends.
Using so that det Ω plays a role similar to the (Hirota) τ -function of the (scalar) 2DTL equation.

Solutions for vanishing seed
The linear system with ϕ 0 = 0 reads It possesses solutions of the form Here θ a , a = 1, . . . , A, are constant m×N matrices, k denotes the discrete variable, P a , a = 1, . . . , A, are constant N × N matrices, and ϑ(P ) = (P − I) x + (I − P −1 ) y . (4.10) Correspondingly, the adjoint linear system takes the form which is solved by The equations for the Darboux potential Ω are reduced to Writing with a constant N × N matrix Ω 0 , it follows that W ba has to satisfy the Sylvester equation and if p i,a = q j,b for all i, j = 1, . . . , N and a = 1, . . . , A, b = 1, . . . , B, then the unique solution is known to be given by the Cauchy-like N × N matrices Assuming that Ω 0 is invertible, Remark 4.3 shows that we can set Ω 0 = 1 N without loss of generality. The remaining transformations, according to Remark 4.3, can be used to reduce the parameters in θ or χ.

Pure solitons
Here we further restrict the class of solutions specified in Section 4.1.1 by setting A = B = 1. Then we set P = P 1 , Q = Q 1 , and assume these matrices to be diagonal, so that (4.13) holds. Furthermore, we write P = diag(p 1,1 , . . . , p N,1 ) =: diag(p 1 , . . . , p N ) , Q =: diag(p 1,2 , . . . , p N,2 ) =: diag(q 1 , . . . , q N ) , where ξ i , i = 1, . . . , N , are constant m-component column vectors and η i , i = 1, . . . , N , are constant n-component row vectors. We assume that p i > 0 and q i > 0, i = 1, . . . , N , since otherwise the generated solution of (4.1) will be singular. With the further assumption that p i = q j for i, j = 1, . . . , N , we have provisionally 10 assuming p > 0, we obtain Let us introduce Instead of using (a 1 , . . . , a N ) as a subscript, we simply write a 1 . . . a N in the following. For example, ϑ a 1 ...a N = ϑ (a 1 ,...,a N ) . From (4.8) we find that the pure soliton solutions of the p2DTL K equation are given by with τ := e ϑ 2 det Ω ,  Besides conditions imposed on p i and q i such that all the ϑ I appearing in (4.14) are real, regularity of a pure N -soliton solution requires µ I ≥ 0 for all I ∈ {1, 2} N , and µ J > 0 for at least one J ∈ {1, 2} N . We will impose the stronger condition µ I > 0 for all I ∈ {1, 2} N . These conditions will be appended to our definition of pure solitons.

It follows that
Example 4.5. For N = 1, writing p 1 = p, q 1 = q, ξ 1 = ξ, η 1 = η and κ = ηKξ, we have the single soliton solution which leads to In terms of the variables t = x + y and z = x − y, it reads We have to restrict the parameters such that p/q > 0 and qκ > 0. The solution becomes independent of z if we choose q = p −1 , in which case the above ϕ reduces to a single soliton solution of the p1DTL K equation (4.3), and we have It is obvious from (4.8) and the sizes of its matrix constituents that, for N = 1, the binary Darboux transformation with zero seed can only yield a rank one solution.

Tropical limit of pure solitons
We define the tropical limit of a matrix soliton solution via the tropical limit of the scalar function τ (cf. [6,7,8]). Let In the region of R 2 × Z, where a phase ϑ I dominates all others, in the sense that log(µ I e ϑ I ) > log(µ J e ϑ J ) for all participating J = I, the tropical limit of the potential ϕ is given by (4.19). 11 These expressions do not depend on the variables x, y, k (respectively, z, t, k). The boundary between the regions associated with the phases ϑ I and ϑ J is determined by the condition (4.20) Not all parts of such a boundary are "visible", in general, since some of them may lie in a region where a third phase dominates the two phases. The tropical limit of a soliton solution, more precisely, of the variable u, has support on the visible parts of the boundaries between the regions associated with phases appearing in τ . For I = (a 1 , . . . , a N ) we set (a 1 , . . . , a j−1 , a, a j+1 , . . . , a N ) .
The j-th soliton (having parameters p j and q j ) lives, in the tropical limit, on the set of twodimensional plane segments determined, via (4.20), by for all I ∈ {1, 2} N . More explicitly, the last equation reads All these plane segments are parallel. In general there are relative shifts between the segments, they do not constitute together a single plane. This gives rise to the familiar (asymptotic) "phase shift" of solitons caused by their interaction. Fig. 3 below shows this for a 2-soliton example, considered at constant time, so that the configuration of planes is projected to a graph in two dimensions. For j = 1, . . . , N , the regularity conditions (4.21) will be assumed in the following. On a (visible) boundary segment, the value of u is given by Instead of the above expressions for the tropical values of u, we will rather consider  As a consequence, we have the normalization tr(Kû IJ ) = 1 . Remark 4.6. We note that tr(Ku I j (1),I j (2) ) = − 1 4 (p I j (1) −p I j (2) )(p I j (1) − p I j (2) ) = − 1 4 which shows that its value is the same everywhere (i.e., for all I) on the tropical support of the j-th soliton.
The above expressions for τ and F coincide with those derived in the KP K case [9]. The only difference is in the expressions for the phases, but the latter do not enter the expressions for the polarizations.
Let us now consider the graph in Fig. 3 as a scattering process evolving from bottom to top. Defining we find We observe that u i and u i all have rank one. They are K-projections, i.e., Furthermore, they satisfy The equations (4.26) imply 3) shows that (u 1 , u 2 ) → (u 1 , u 2 ) provides us with a realization of the Yang-Baxter map R. 13 (1, 2) on a threefold direct sum. Also see [9] for the case of the vector KP K equation. Introducing we have R(i, j) = S(i, j)R(i, j)S(i, j) −1 , andR(i, j) also satisfies the Yang-Baxter equation. We further note that RA(i, j) If we drop the normalization condition (4.24) in the computation of the Yang-Baxter map, the resulting R-matrix turns out to be of the latter form. The same holds if we consider v = ϕ + − ϕ instead ofû.

Yang-Baxter and non-Yang-Baxter maps at work
In Section 4.3 we looked at the relation between the polarizations associated with the boundary segments of dominant phase regions of a pure 2-soliton solution, selecting a "propagation direction". But there is actually no preferred direction. It is therefore more adequate to regard (4.26) just as determining a relation between four polarizations, and there are several ways in which this determines a map from two "incoming" to two "outgoing" polarizations. Example 4.9. We keep the choices for K, ξ i and η i , made in Example 4.7. Then κ ij = 1 and α 12 = (q 2 − q 1 )/(q 2 − p 1 ), so that the regularity conditions (4.21) read p 2 , q 2 > 0, p 1 (q 2 − q 1 )/(q 2 − p 1 ) > 0 and p 1 /q 1 > 0. Furthermore, without loss of generality, we can choose the parameters such that, for large enough negative value of k, soliton 1 appears in z-direction to the left of soliton 2.
1. q 1 < p 1 < p 2 < q 2 or q 1 < p 2 < q 2 < p 1 . In this case the phase constellation is that shown in Fig. 3. Regarding it as a process in k-direction, the map of polarizations is the Yang-Baxter map R.
2. p 1 < q 1 < q 2 < p 2 or p 1 < q 2 < p 2 < q 1 . The phase constellation is that shown in the first plot of Fig. 4 (which is generated with p 1 = 1/4, q 1 = 3/2, q 2 = 2 and p 2 = 3). The map of polarizations, in k-direction, leads us to the inverse of the Yang-Baxter map R, which is also a Yang-Baxter map.
4. p 1 < q 1 < p 2 < q 2 or p 1 < p 2 < q 2 < q 1 . The phase constellation is that shown in the third plot of Fig. 4 (which is generated with p 1 = 1/4, q 1 = 3/2, p 2 = 2 and q 2 = 3). In this case, the map of polarizations, from bottom to top, is the inverse of T .
For any one of the plots in Figs. 3 or 4, we obtain realizations of all the maps by choosing different directions.
The naive expectation that 2-soliton scattering yields a Yang-Baxter map is therefore wrong. But we have to keep in mind that the Yang-Baxter property is a statement about three solitons. A 3-soliton solution involves three 2-particle interactions. In the tropical limit, this means that a composition of three of the above maps carries the polarizations along the tropical limit support. We will see in the next subsection that the Yang-Baxter equation is indeed only required to hold for a mixture of the maps, but not for each map separately. Remark 4.10. As seen above, the constellation of dominant phase regions depends on the concrete values of the parameters. If, for a certain constellation, we select a direction and obtain a Yang-Baxter map, then the latter has the Yang-Baxter property for all choices of parameters. From the above, we conclude that the map, relating the (relativ to our choice of direction) incoming and outgoing polarizations, is a Yang-Baxter map if and only if the phase region that lies between the two incoming solitons is that of ϑ 12 or ϑ 21 .  The regularity condition (4.21) is then (p 2 − p 1 )/(p 2 − 1) > 0.
1. 0 < p 2 < p 1 < 1. The map of incoming to outgoing polarizations, in t-direction, is the reduced Yang-Baxter map. The tropical limit graph, for parameters p 1 = 1/2 and p 2 = 1/10 is shown in the first plot of Fig. 5.
2. 1 < p 1 < p 2 . In this case, the polarization map is the inverse of the (reduced) Yang-Baxter map. For p 1 = 2 and p 2 = 10, we obtain the second plot of Fig. 5. 3. 0 < p 1 < 1 < p 2 . In this case, the (reduced) map T is at work. For p 1 = 1/2 and p 2 = 10, the third plot of Fig. 5 is obtained.

Pure 3-soliton solution
For N = 3 we find τ = κ 11 κ 22 κ 33 β e ϑ 111 + κ 11 κ 22 α 12 e ϑ 112 + κ 11 κ 33 α 13 e ϑ 121 +κ 22 κ 33 α 23 e ϑ 211 + κ 11 e ϑ 122 + κ 22 e ϑ 212 + κ 33 e ϑ 221 + e ϑ 222 , and Again, we set κ ij = η i Kξ j . Furthermore, we have Note that α ij = α ijj . The above expressions coincide with those obtained in the KP K case [9]. The only difference is in the phases entering the exponentials. Here we wrote F in a more compact form.  Regarding it as an evolution of three particles (solitons in the tropical limit) in k-direction, first particle 1 meets particle 2, then particle 1 and particle 3 meet, and finally particle 2 interacts with particle 3. This corresponds to one side of the (mixed) Yang-Baxter equation. The right graph is obtained for t = 50 and the sequence of interactions corresponds to the other side of the (mixed) Yang-Baxter equation.
The polarizations at a node of the tropical limit graph at t = t 0 are completely determined if we know the soliton numbers associated with two "incoming lines", e.g. soliton 1 with the left and soliton 2 with the right line, and the number of the enclosed phase region, for example 222. In this configuration, to the left of incoming soliton line 1 we have the phase with number 122, to the right of incoming soliton line 2 the phase with number 212, and the remaining one has necessarily number 112. Hence there are incoming polarizationsû 122,222 ,û 212,222 and outgoing polarizationŝ u 112,212 ,û 112,122 of solitons 1, 2. These can be computed from the above 3-soliton solution and it can be verified that they are related by one of the maps considered in Section 4.4. This also holds for all possible other nodes.
In the following, we consider the situation shown in the left plot in Fig. 6, regarding it as a process from bottom to top. Accordingly, we set where a subscript m1 indicates that the line segment appears in the middle of the first plot in Fig. 6. For the right plot in Fig. 6, we only have to replace the polarizations of middle segments by u 1,m2 =û 122,222 , u 2,m2 =û 222,212 , u 3,m2 =û 221,222 . Now we can check that either the Yang-Baxter map R, the map T , or one of their inverses acts at each crossing of the tropical limit graph. Proceeding from bottom to top in the left plot of Fig. 6, the first crossing involves only solitons 1 and 2, so we can ignore the last digit of the numbers of the involved phases. The situation is that of the 2-soliton interaction sketched in the first drawing of Fig. 7. Accordingly, T , given by (3.5), should map (u 1,in , u 2,in ) to (u 1,m1 , u 2,m1 ). Indeed, this can be varified using the above data.  Fig. 6, proceeding from bottom to top. In the phase numbers we marked (with red color) the digit corresponding to the soliton that does not take part in the respective interaction. Disregarding this digit, the drawing describes the phase constellation of a 2-soliton interaction. In this way, for example, the first drawing corresponds to an application of T .
The next crossing, where solitons 1 and 3 interact, is sketched as a 2-soliton interaction in the second drawing of Fig. 7. In this situation, the Yang-Baxter map R should yield (u 1,m1 , u 3,in ) → (u 1,out , u 3,m1 ). Indeed, we find (which can be deduced from (D.1)). Similarly, Finally, at the last crossing solitons 2 and 3 interact. The situation is sketched in the last drawing in Fig. 7. Accordingly, we expect the map T −1 to be at work and this can indeed be verified. The 2-soliton subinteractions appearing in the right plot of Fig. 6 are sketched in Fig. 8. For the lowest, corresponding to the first drawing in Fig. 8, we can verify (see Appendix D) that the polarizations are related by the map T −1 , given by (3.7). The second drawing in Fig. 8 corresponds to an application of the Yang-Baxter map R, the third to an application of T .
At least for the parameter range, for which the tropical limit graphs at large negative, respectively positive time t have the structure shown in Fig. 6, we can now conclude that (3.6) holds. This is so because the polarizations associated with line segments of the tropical limit graph do not depend on the variables z, t, k. Since the two situations shown in Fig. 6 belong to the same solution of the 2DTL K equation, the result of the application of the two sides of (3.6) are the same. We know that (3.6) indeed holds for all parameter values (provided that p i , q i are all different).
Corresponding tropical limit graphs are shown in Fig. 10, from which we read off  Figure 10: Tropical limit graphs of yet another pure 3-soliton solution of the p2DTL K equation for a negative (left graph) and a positive (right graph) value of t, using the data of Example 4.14.

Conclusions
We presented a "Lax representation" for the Yang-Baxter map R obtained in [9] from the pure 2-soliton solution of the (K-modified) matrix KP equation. For the matrix KdV reduction, such a Lax representation had been found earlier in [2], also see [18]. Using also the inverse of the Lax matrix A, we were led to entwining Yang-Baxter maps, in the sense of [16]. Besides the Yang-Baxter map, this involves another map, which is not Yang-Baxter, but both maps satisfy a "mixed Yang-Baxter equation".
We demonstrated that this structure is indeed realized in 3-soliton interactions. Here we concentrated on an analysis of matrix generalizations of the two-dimensional Toda lattice equation, of which solitons can be generated via a binary Darboux transformation. For the subclass of "pure solitons", which only exhibit elastic scattering (i.e., no merging or splitting of solitons), we elaborated the tropical limit of the general 2-and 3-soliton solution. It turned out that exactly the same Yang-Baxter map as in the KP K case is a work. The crucial new insight is that the flow of polarizations for three solitons is not, in general, described by the Yang-Baxter map alone, but by entwining Yang-Baxter maps. This also holds for KP K and even for KdV K (with K = 1 n , for example), also see Appendix A. This result crucially relies on our tropical limit analysis of solitons.
More precisely, in the preceding section we found that (1.3), (3.6) and also (3.12) are realized by pure 2DTL K solitons. We do not know whether this is also so for the remaining mixed Yang-Baxter equations in Section 3.
Since the soliton with number i is specified via the parameters p i , q i , and the polarization u i , we can associate the matrices A i (λ, u i ) and B i (λ, u i ) with it (as well asÃ i (λ, u i ) andB i (λ, u i )). Comparing the three-fold products of these Lax matrices, which imply a mixed Yang-Baxter equation, with the tropical limit plots for pure 3-solitons, one observes the following rule. Whether A i or B i is at work, depends on the constellation of phases to the left and to the right of the line. For example, if the first soliton has phase 1ab to the left and 2cd to the right (a, b, c, d ∈ {1, 2}), we have to choose A 1 , whereas with 2ab to the left and 1cd to the right, it has to be B 1 . The 3-soliton interaction shown in the first plot in Fig. 6 corresponds to (hiding away the parameters). The Lax matrices, which are subject to the refactorization equations (2.2), generate a Zamolodchikov-Faddeev algebra. In S-matrix theory of an integrable QFT, these matrices play a role as creation and annihilation operators (see, e.g., Section 4.2 in [1], and references cited there).
We also found that the two maps R, T , and their inverses, provide us with solutions of the "WXZ system" [20], called Yang-Baxter system in [21].
We further note that the Yang-Baxter map obtained for the matrix Nonlinear Schrödinger (NLS) equation in [3,4,5] has the same form as the Yang-Baxter map R for matrix KP and matrix two-dimensional Toda lattice.
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 [S, S −1 ], where A 0 is the algebra of smooth functions of two variables, x and y, and also dependent on a discrete variable on which the shift operator S acts. ζ 1 , ζ 2 constitute 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 Ω. Setting φ = ϕ S −1 , the equation (C.1) is equivalent to the p2DTL K equation (4.1). Choosing a solution φ 0 = ϕ 0 S −1 and setting ∆ = Γ = S −1 , κ = λ = 0 , the linear system (C.2) and the adjoint linear system (C.3) lead to (4.5) and (4.6), respectively. Furthermore, via Ω → ΩS, (C.4) implies (4.7). According to the theorem, (C.5) yields a new solution of the p2DTL K equation (4.1).