Physical meaning of the Yang-Baxter equation I'm a graduate student in mathematics, and I have lately been interested in the relation between knot theory and statistical mechanics. As I understood, the Yang-Baxter equation (shown below) is the equivalent of the Reidemeister III move (RIII), and it appeared as a convenient hypothesis to solve lattice models such as the Ising model.

What I can't find anywhere though is any kind of convincing physical interpretation of the Yang-Baxter equation. As I have little physics education, I'm not looking for the whole technical story, but rather some convincing argument that it is a sensible relation to assume. Here is one I came up with:

Thinking of the three strands of RIII as the trajectories of three particles $p_1$, $p_2$ and $p_3$, the RIII move means that It does not matter in which order particles interact: either in the order $(p_1,p_2)$, $(p_1,p_3)$ and $(p_3,p_2)$ or in the order $(p_2,p_3)$, $(p_1,p_3)$ and $(p_1,p_2)$.

This is highly unconvincing: why is there a relation between this orders, and not other orders? How is the topology of the crossings (whether crossings are positive or negative) taken into account? So my question is: do you have a better motivation for assuming the Yang-Baxter equation?
NB: as an example of a simple motivation, I convinced myself that the $q$-Potts model ($E\sim\sum_{(i,j)} \delta(s_i,s_j)$, where $s_i$ is the state of vertex $i$ and the sum is over all edges in the lattice) was sensible by thinking of magnets whose positions are fixed at some points: they lowest energy state is the one where all poles are in opposite direction, justifying the term $\delta(s_i,s_j)$.

Question following Yu-v answer. (Sorry to keep asking, but I feel I haven't get the point yet.) The part I don't get is 

As each of this orderings must have a well-defined amplitude, different ways to relate them to each other must agree.

Here is my understanding of the situation : the system is the superposition of the solutions of each ordering. Moreover, we have a way to relate the solutions of an ordering $S_{12}$ to the solution of the ordering $S_{21}$ obtained from $S_{12}$ by switching the particles 1 and 2. Call this relation $R$ and write $S_{12}\cong_R S_{21}$ (here I just mean that there is a one-to-one correspondence between the solutions). $S_{21}$ is not defined from $S_{21}$ by $R$: all solutions exist on their own, and it just so happened that there is a relation between them.
Now we can define relations $R_1$ and $R_2$ between $S_{123}$ and $S_{321}$ according to the two ways of reordering the particles 1, 2 and 3: $S_{123}\cong_{R_1} S_{321}$ and $S_{123}\cong_{R_2} S_{321}$. Only considering what I just said, there is no reason to think that $R_1=R_2$. For example, you can have different isomorphisms between two vector spaces: having different isomorphisms does not prevent the vector spaces from being well-defined. I deduce that I'm missing an essential point.
 A: I will try to give a reason, based on interacting quantum models in 1d. Basically we have an equation (Schrödinger's equation) that we need to solve, which means that we are looking for a function $\Psi(x_1, x_2, \ldots, x_N)$ that will satisfy it. In the function above, the range of the positions of the particles is some $[-L/2, L/2]$, which mean that we can divide it into many different sectors, each corresponding to a specific ordering of the particles
$$ \Psi(x_1, \ldots, x_N) = \sum_{\mathcal{S}}\theta(x_{i_1} \leq x_{i_2} \leq x_{i_3} \leq \ldots \leq x_{i_N}) F_S(x_1, \ldots, x_N) $$
and $S$ runs over the $N!$ permutations of the particles. Here, the $\theta$ function is a short-hand for multiplication of successive $\theta$ functions ensuring that in this sector $x_{i_j} \leq x_{i_{j+1}}$ for the specific permutation.
The Schrödinger equation now tells us how different sectors are related to each other. When we exchange the place of two particles, that is going from a sector where $x_3 \leq x_4$ to the sector where $x_4 \leq x_3$, the particles interact (this is the case when the interaction is short-ranged between pairs of particles). So we can relate $F_S$ to $F_{S'}$ if they differ by exchanging of two particles, by solving the Schrödinger equation.
The problem is, that now if I look on rearranging three particles I have two ways to relate, say, $F_{S_{123}}$ (where $x_1\leq x_2 \leq x_3$) to $F_{S_{321}}$, (where $x_3\leq x_2 \leq x_1$):
$$ F_{S_{123}} \to F_{S_{213}} \to F_{S_{231}}\to F_{S_{321}}$$
and 
$$ F_{S_{123}} \to F_{S_{132}} \to F_{S_{312}} \to F_{S_{321}}$$
with each arrow must satisfy the condition set by Schrödinger equation. In order for our solution to be consistent, both ways must lead to the same result. This is basically what the Yang-Baxter equation gives us, and this is why we must satisfy it. It is guaranteed that if we satisfy it, the generalization to $N$ particles is also consistent.
So to conclude - the YB equation gives us a consistency condition for our solution. Without it, there is no way to uniquely decide on the wave function in a given sector.

Edit following addition to question
The function we are looking for is not a superposition of different solutions. The superposition itself is the solution. There is a unique solution to the equation, and this solution is derived by matching the different amplitudes between the different orderings.
I will give a more explicit example, by considering the Lieb-Liniger model. We are trying to find an eigenfunction of the Hamiltonian operator
$$ H = -\sum_{i} \partial_i^2 + c \sum_{i<j}\delta(x_i-x_j)$$
for some $N$ variables $x_1, \ldots , x_N$. When $x_i \neq x_j$ for all $i$ and $j$ the solution is immediate $F_S(x_1, \ldots , x_N) = A_S \exp(-i \sum k_j x_j)$ which is just to say that when the particles do not "touch" one another they are free waves. So the solution is of the form
$$ \Psi(x_1, \ldots, x_N) = \sum_{S} A_S \theta(S) e^{-i \sum k_j x_j} $$
where $\theta(S)$ is the product of the theta-functions (that is, it is $1$ if the order of the particles is according to $S$ and zero otherwise).
However, in order to satisfy the $\delta$-function part of the Hamiltonian (that is, the interaction), the set of coefficients $A_S$ cannot be chosen freely, but must be related, and we can derive that when we compare the coefficients $A_S$ and $A_{S'}$ that are related to each other by exchanging of particle $x_i$ and $x_j$ (that is, in $S$ we have $x_i < x_j$ and in $S'$ the opposite), then $A_S = P_{ij} A_{S'}$. In this case $P_{i,j}$ is some phase that depends on $k_i-k_j$. Now in order for the solution to be valid, we must have for example $P_{ji}P_{ij}=1$, since $A_S = P_{ji}A_{S'} = P_{ji}P_{ij}A_{S}$. The Yang-Baxter equation is a similar consistency condition, but on different exchanges of 3 particles
$$ P_{ij}P_{ik}P_{jk}P_{ji}P_{ki}P_{kj}=1$$
Turns out that these two conditions are sufficient and maintaining them ensures that all ordering of $N$ different particles can be connected safely by the appropriate $P$. For the Lieb-Liniger model the Yang-Baxter equation is trivial since $P$ are just $c$-numbers but for more complicated models where different $P$ not necessarily commute finding the correct form of $P$ that maintains Yang-Baxter is quite challenging.
