Integrable couplings [1, 2] are coupled systems of integrable equations, which has been introduced when we study of Virasoro symmetric algebras. It is an important topic to look for integrable couplings because integrable couplings have much richer mathematical structures and better physical meanings. In recent years, many methods of searching for integrable couplings have been developed [3-11], especially authors have used the enlarged matrix spectral problem method [3] to find out integrable couplings. For example, they ever used the following the spectral matrix to obtain integrable couplings

$ \begin{equation} \bar{U}=\begin{bmatrix} U(u)&0 \\ U_{a}(v)& U(u) \\ \end{bmatrix}, \end{equation} $

(1.1)

where the sub-spectral matrix $U$ is associated with a given integrable equation $u_{t}=K(u)$. However, soon afterwards, researcher find out that obtained integrable couplings is a relatively simple. So in order to get better integrable couplings to the known integrable system, we need to introduce an enlarged relatively complex spectral matrix

$ \begin{equation} \bar{U}=\begin{bmatrix} U(u)&0 \\ U_{a}(v)&U(u)+U_{a}(v) \\ \end{bmatrix}. \end{equation} $

(1.2)

Therefore, from zero curvature equation

$ \begin{equation} \bar{U}_{t}-\bar{V}_{x}+[\bar{U}, \bar{V}]=0, \end{equation} $

(1.3)

where

$ \begin{equation} \bar{V}=\bar{V}(\bar{u})=\begin{bmatrix} V(u)&0 \\ V_{a}(\bar{u})&V(u)+V_{a}(\bar{u}) \\ \end{bmatrix}, \end{equation} $

(1.4)

and $\bar{u}$ consist of $u$ and $v$, we can give rise to

$ \begin{eqnarray} \left\{% \begin{array}{ll} U_{t}-V_{x}+[U, V]=0, \\ U_{a, t}-V_{a, x}+[U, V_{a}]+[U_{a}, V]+[U_{a}, V_{a}]=0.\\ \end{array}% \right. \end{eqnarray} $

(1.5)

This is an integrable couplings of (1.1), due to (1.3), and it is a nonlinear integrable coupling because the commutator $[U_{a}, V_{a}]$ can generate nonlinear terms.

Let us further take a solution $\bar{W}$ to the enlarged stationary zero curvature equation

$ \begin{equation} \bar{W}_{x}=[\bar{U}, \bar{W}]. \end{equation} $

(1.6)

Then, we use the quadratic-form identity or variational identity [6, 7]

$ \begin{equation} \frac{\delta}{\delta\bar{u}}\int\langle\bar{W}, \bar{U}_{\lambda}\rangle dx=\lambda^{-\gamma}\frac{\partial}{\partial\lambda}\lambda^{\gamma} \langle\bar{W}, \bar{U}_{\bar{u}}\rangle, \end{equation} $

(1.7)

where $\gamma$ is a constant, to search for the Hamiltonian structures for the integrable couplings [7]. In the variational identity (1.7), $\langle., .\rangle$ is non-degenerate, symmetric and ad-invariant bilinear form over the Lie algebra

$ \begin{equation} \bar{g}=\{\begin{bmatrix} A&0 \\ B&A+B \\ \end{bmatrix}|A, B\in g\}. \end{equation} $

(1.8)

As is well known, the conservation laws play an important roles on discussing the integrability for soliton equation. Since Miura, Gardner, and Kruscal's discovery [13] of an infinite number of conservation laws for KdV equation, many methods have been developed to find them, mainly due to the contribution of Wadati et al [13-15]. Many papers dealing with symmetries and conservation laws were presented, the direct contribution method of multipliers for the conservation laws was presented [16]. Comparatively, the less nonlinear integrable couplings of the soliton equations have been considered for their conservation laws. In what follows, we will make above idea to apply the Guo hierarchy.

For the Guo spectral problem [12]

$ \begin{equation} \begin{split} \begin{array}{ll} \phi_x=U\phi=U(u, \lambda)\phi, \phi_t=V\phi, \lambda_t=0, \\ U=\frac{1}{2}\left(% \begin{array}{ccc} 1/\lambda&q+r \\ q-r&-1/\lambda \\ \end{array}% \right), u=\left(% \begin{array}{cc} q \\ r \\ \end{array}% \right).\\ \end{array}% \end{split} \end{equation} $

(2.1)

Setting

$ \begin{eqnarray*} V=\frac{1}{2}\left(% \begin{array}{ccc} a&b+c \\ b-c&-a \\ \end{array}% \right)=\frac{1}{2}\sum\limits_{n\geq0}\left(% \begin{array}{ccc} a_n&b_n+c_n \\ b_n-c_n&-a_n \\ \end{array}% \right)\lambda^{n}.\\ \end{eqnarray*} $

(2.2)

The stationary zero curvature equation $V_{x}=[U, V]$ yields that

$ \begin{eqnarray} \begin{array}{ll} a_{nx}=rb_{n}-qc_{n}, \\ b_{nx}=c_{n+1}-ra_{n}, \\ c_{nx}=b_{n+1}-qa_{n}.\\ \end{array}% \end{eqnarray} $

(2.3)

Choose the initial data

$ \begin{equation} a_0=1, \ \ \ b_0=c_0=0, \end{equation} $

(2.4)

then we have

$ \begin{eqnarray} \begin{array}{ll} a_{1}=0, b_{1}=q, c_{1}=r, \\ a_{2}=\frac{1}{2}(r^{2}-q^{2}), b_{2}=r_{x}, c_{2}=q_{x}, \\ a_{3}=q_{x}r-qr_{x}, b_{3}=q_{xx}+\frac{1}{2}q(r^{2}-q^{2}), \\ c_{3}=r_{xx}+\frac{1}{2}r(r^{2}-q^{2}), \cdots.\\ \end{array}% \end{eqnarray} $

(2.5)

From the compatibility conditions of the following problems

$ \begin{eqnarray} \begin{array}{ll} \phi_x=U\phi, \quad \phi_t=V^{(n)}\phi, \\ V^{(n)}=(\lambda^{-n}V)_{-}+\Delta_{1n}, \end{array}% \end{eqnarray} $

(2.6)

where $\Delta_{1n}=0$. We can determine the Guo hierarchy of soliton equations

$ \begin{equation} u_{t_n}=\left(% \begin{array}{cc} q \\ r \\ \end{array}% \right)_{t_n}=K_n=\left(% \begin{array}{cc} c_{n+1} \\ b_{n+1}\\ \end{array}% \right)=J\left(% \begin{array}{cc} b_{n} \\ - c_{n} \\ \end{array}% \right)=JL^n\left(% \begin{array}{cc} b_{1} \\ -c_{1} \\ \end{array}% \right), n\geq0, \end{equation} $

(2.7)

with the Hamiltonian operator $J$, the hereditary recursion operator $L$ and the Hamiltonian functions $H_{n}$, respectively, as follows

$ \begin{equation} J=\left(% \begin{array}{cc} \partial+r\partial^{-1}r&r\partial^{-1}q \\ q\partial^{-1}r&-\partial+q\partial^{-1}q \end{array}% \right), L=\left(% \begin{array}{cc} q\partial^{-1}r &-\partial+q\partial^{-1}q \\ -\partial-r\partial^{-1}r&-r\partial^{-1}q \\ \end{array}% \right), \\ \\ H_{n}=\int-\frac{a_{n+1}}{n}dx, ~~~~n\geq 0. \end{equation} $

(2.8)

Let us now begin with an enlarged spectral matrix

$ \begin{equation} \bar{U}=\left(% \begin{array}{ccc} U&0 \\ U_{a}&U+U_{a} \end{array}% \right), U_a=\frac{1}{2}\left(% \begin{array}{cc} 0& p_{1}+p_{2} \\ p_{1}-p_{2}&0 \\ \end{array}% \right), \\ \end{equation} $

(2.9)

$ \begin{equation} \bar{V}=\left(% \begin{array}{cc} V&0 \\ V_{a}& V+V_{a} \end{array}% \right), V_a=\frac{1}{2}\left(% \begin{array}{ccc} e&f+g \\ f-g&-e \\ \end{array}% \right)\\ \end{equation} $

(2.10)

with the help of the corresponding enlarged stationary zero curvature equation $\bar{V}_{x}=[\bar{U}, \bar{V}]$, we have

$ \begin {equation}V_{ax}=[U, V_{a}]+[U_{a}, V]+[U_{a}, V_{a}], \end{equation} $

(2.11)

which equivalently generates

$ \begin{eqnarray} \begin{array}{ll} e_{x}=-qg+rf-p_{1}c+p_{2}b-p_{1}g+p_{2}f, \\ f_{x}=\lambda^{-1}g-re-p_{2}a-p_{2}e, \\ g_{x}=\lambda^{-1}f-qe-p_{1}a-p_{1}e.\\ \end{array}% \end{eqnarray} $

(2.12)

Setting

$ \begin{eqnarray*} V_a=\frac{1}{2}\left(% \begin{array}{ccc} e&f+g \\ f-g&-e \\ \end{array}% \right)=\frac{1}{2}\sum\limits_{n\geq0}\left(% \begin{array}{ccc} e_{n}&f_{n}+g_{n} \\ f_{n}-g_{n}&-e_{n} \\ \end{array}% \right)\lambda^{n}.\\ \end{eqnarray*} $

(2.13)

And then, the (2.12) can be transformed into

$ \begin{eqnarray} \begin{array}{ll} e_{nx}=-qg_n+rf_n-p_{1}c_n+p_{2}b_n-p_{1}g_n+p_{2}f_n, \\ f_{nx}=g_{n+1}-re_n-p_{2}a_n-p_{2}e_n, \\ g_{nx}=f_{n+1}-qe_n-p_{1}a_n-p_{1}e_n.\\ \end{array}% \end{eqnarray} $

(2.14)

We choose the initial data

$ \begin{equation} e_{0}=1, f_{0}=g_{0}=0, \end{equation} $

(2.15)

then we have

$ \begin{eqnarray} \begin{array}{ll} e_{1}=0, f_{1}=q+2p_{1}, g_{1}=r+2p_{2}, \\ e_{2}=\frac{1}{2}(r^{2}-q^{2})-2qp_{1}+2rp_{2}-p_{1}^{2}+p_{2}^{2}, \\ f_{2}=r_{x}+p_{2, x}, g_{2}=q_{x}+p_{1, x}, \\ e_{3}=2(rp_{1, x}-p_{1}r_{x}+p_{2}p_{1, x}-p_{1}p_{2, x}+p_{2}q_{x}-qp_{2, x})+rq_{x}-qr_{x}+qrp_1, \\ f_{3}=q_{xx}+2p_{1, xx}+\frac{1}{2}p_{1}(r^{2}-q^{2})+(q+p_{1})(\frac{1}{2}r^{2}-\\ \frac{1}{2}q^{2}-2qp_{1}+2rp_{2}-p_{1}^{2}+p^{2}_{2}), \\ g_{3}=r_{xx}+2p_{2, xx}+\frac{1}{2}p_{2}(r^{2}-q^{2}) +(r+p_{2})(\frac{1}{2}r^{2}-\frac{1}{2}q^{2}-2qp_{1}+\\ 2rp_{2}-p_{1}^{2}+p^{2}_{2}), \cdots.\\ \end{array}% \end{eqnarray} $

(2.16)

Using the zero curvature equation

$ \begin{eqnarray} \left\{% \begin{array}{ll} \bar{U}_{t_n}-\bar{V}_{x}+[\bar{U}, \bar{V}]=0, \\ U_{a{t_n}}=V^{(n)}_{ax}+\Delta_{2n, x}-[U, V_a^{(n)}]-[U, \Delta_{2n}]-[U_a, V^{(n)}]-[U_a, V^{(n)}_a]-[U_a, \Delta_{2n}], \\ \end{array}% \right. \end{eqnarray} $

(2.17)

where $\Delta_{2n}=0$.

We can rewrite (2.17) as

$ \begin{equation} \left(% \begin{array}{cc} p_{1} \\ p_{2} \\ \end{array}% \right)_{t_n}=\left(% \begin{array}{cc} g_{n+1}\\ f_{n+1}\\ \end{array}% \right).\\ \end{equation} $

(2.18)

Also, we have the following results

$ \begin{equation} \bar{u}_{t_n}=\left(% \begin{array}{ccc} q \\ r \\ p_{1} \\ p_{2} \\ \end{array}% \right)_{t_n}=\left(% \begin{array}{ccc} c_{n+1}\\ b_{n+1}\\ g_{n+1}\\ f_{n+1} \\ \end{array}% \right), n\geq 0. \end{equation} $

(2.19)

Obviously, when $p_1=p_2$=0 in (2.19), the above results become (2.7). So we can say (2.19) is integrable couplings of the Guo hierarchy. When $n=2$, we have

$ \begin{eqnarray} \left\{% \begin{array}{ll} q_{t_2}=r_{xx}+\frac{1}{2}r(r^{2}-q^{2}), \\ r_{t_2}=q_{xx}+\frac{1}{2}q(r^{2}-q^{2}), \\ p_{1, t_2}=r_{xx}+2p_{2, xx}+\frac{1}{2}p_2(r^{2}-q^{2})+\\ (r+p_{2})(\frac{1}{2}r^{2}-\frac{1}{2}q^{2}-2qp_{1}+2rp_{2}-p_{1}^{2}+p^{2}_{2}), \\ p_{2, t_2}=q_{xx}+2p_{1, xx}+\frac{1}{2}p_1(r^{2}-q^{2})+(q+p_{1})(\frac{1}{2}r^{2}-\\ \frac{1}{2}q^{2}-2qp_{1}+2rp_{2}-p_{1}^{2}+p^{2}_{2}).\\ \end{array}% \right. \end{eqnarray} $

(2.20)

So, we can say that the system in (2.19) with $n\geq 2$ provide a hierarchy of nonlinear integrable couplings for the Guo hierarchy of the soliton equation. Now, we proceed to search for the Hamiltonian structure of the equations hierarchy (2.19).

To construct Hamiltonian structures of the integrable couplings obtained, we need to compute non-degenerate, symmetric and invariant bilinear forms on the following Lie algebra

$ \begin{equation}\label{Eq-32} \bar{g}=\{\begin{bmatrix} A &0 \\ B&A+B \\ \end{bmatrix}|A, B\in gl(2)\}. \end{equation} $

(3.1)

For computations convenience, we transform this Lie algebra $\bar{g}$ into a vector from through the mapping

$ \begin{equation}\label{Eq-33} \delta:\bar{g}\rightarrow R^{6}, A\mapsto(a_{1}, a_{2}, \cdots, a_{6})^{T}, A=\begin{bmatrix} a_{1}&a_{2}&0&0 \\ a_{3}&-a_{1}&0&0 \\ a_{4}&a_{5}&a_{1}+a_{4}&a_{2}+a_{5} \\ a_{6} &-a_{4}&a_{3}+a_{6}&-a_{1}-a_{4} \\ \end{bmatrix} \in \bar{g}. \end{equation} $

(3.2)

The mapping $\delta$ induces a Lie algebraic structure on $R^{6}$, isomorphic to the matrix Lie algebra $\bar{g}$ above. It is easy to see that the corresponding commutator [., .] on $R^{6}$ is given by

$ \begin{equation} [a, b]^{T}=a^{T}R(b), a=(a_{1}, a_{2}, \cdots, a_{6})^{T}, b=(b_{1}, b_{2}, \cdots, b_{6})^{T}\in R^{6}, \end{equation} $

(3.3)

where

$ \begin{equation} R(b)=\begin{pmatrix} 0&2b_{2}&-2b_{3}&0&0&0 \\ b_{3}&-2b_{1}&0 &0&0&0 \\ -b_{2}&0&2b_{1}&0& 0&0 \\ 0&2b_{5}& -2b_{6}&0&2b_{2}+2b_{5}&-2b_{3}-2b_{6} \\ b_{6}&-2b_{4}&0&b_{3}+b_{6}&-2b_{1}-2b_{4}&0 \\ -b_{5}& 0&2b_{4}&-b_{2}-b_{5}&0&2b_{1}+2b_{4} \\ \end{pmatrix}. \end{equation} $

(3.4)

Define a bilinear form on $R^{6}$ as follows

$ \begin{equation} \langle a, b\rangle=a^{T}Fb, \end{equation} $

(3.5)

where $F$ is a constant matrix, which is main idea by Zhang and Guo presented in 2005 [6].

Then the symmetric property $\langle a, b\rangle=\langle b, a\rangle$ and the ad-invariance property under the Lie product

$ \begin{equation} \langle a, [b, c]\rangle=\langle [a, b], c\rangle, \end{equation} $

(3.6)

requires that $F^{T}=F$ and

$ \begin{equation} (R(b)F)^{T}=-R(b)F, ~~~~ \forall b\in R^{6}. \end{equation} $

(3.7)

So we can obtain

$ \begin{equation} F=\begin{bmatrix} 2\eta_{1}&0&0&2\eta_{2}&0&0 \\ 0&0&\eta_{1}&0&0&\eta_{2} \\ 0&\eta_{1}&0&0&\eta_{2}&0 \\ 2\eta_{2}&0&0&2\eta_{2}&0&0 \\ 0&0&\eta_{2}&0&0&\eta_{2} \\ 0&\eta_{2}&0&0&\eta_{2}&0 \\ \end{bmatrix}, \end{equation} $

(3.8)

where $\eta_{1}$ and $\eta_{2}$ are arbitrary constants.

Therefore, a bilinear form on the underlying Lie algebra $\bar{g}$ is defined by

$ \begin{equation} \langle a, b\rangle_{\bar{g}}=(2b_{1}\eta_{1}+2b_{4}\eta_{2})a_{1}+(b_{3}\eta_{1}+b_{6}\eta_{2})a_{2}+(b_{2}\eta_{1}+b_{5}\eta_{2})a_{3}\\ +(2b_{1}\eta_{2}+2b_{4\eta_{2}})a_{4}+(b_{3}\eta_{2}+b_{6}\eta_{2})a_{5}+(b_{2}\eta_{2}+b_{5}\eta_{2})a_{6}, \\ \end{equation} $

(3.9)

where

$ \begin{equation} A=\begin{bmatrix} a_{1}&a_{2}&0&0\\ a_{3}&-a_{1}&0 &00 \\ a_{4} &a_{5} &a_{1}+a_{4}&a_{2}+a_{5} \\ a_{6}& -a_{4}&a_{3}+a_{6}&-a_{1}-a_{4} \\ \end{bmatrix}, B=\begin{bmatrix} b_{1}&b_{2} &0& 0 \\ b_{3}&-b_{1}&0&0 \\ b_{4}&b_{5}&b_{1}+b_{4}&b_{2}+b_{5} \\ b_{6}&-b_{4}&b_{3}+b_{6}&-b_{1}-b_{4}. \\ \end{bmatrix}. \end{equation} $

(3.10)

It is non-degenerate if and only if

$ \begin{equation} (\eta_{1}-\eta_{2})\eta_{2}\neq 0. \end{equation} $

(3.11)

Based on (3.9), (2.9) and (2.10), we can easily compute that

$ \begin{eqnarray}\nonumber \begin{array}{ll} \langle \bar{V}, \bar{U}_{\lambda}\rangle=-\frac{1}{2\lambda^{2}}(a\eta_{1}+e\eta_{2}), \\ \langle \bar{V}, \bar{U}_{\bar{u}}\rangle =(\frac{1}{2}(b\eta_{1}+f\eta_{2}), -\frac{1}{2}(c\eta_{1}+g\eta_{2}), \frac{1}{2}\eta_{2}(b+f), -\frac{1}{2}\eta_{2}(c+g), \\ \gamma=-\frac{\lambda}{2}\frac{d}{d\lambda}ln|\langle\bar{V}, \bar{V}\rangle|=0. \end{array}% \end{eqnarray} $

By applying the operator $\Gamma_{n+2}$ to both sides of variational identity (1.7) we deduce that

$ \begin{equation} \frac{\delta}{\delta\bar{u}}\int\frac{\eta_{1}a_{n+2}+\eta_{2}e_{n+2}}{n+1}dx= [\eta_{1}b_{n+1}+\eta_{2}f_{n+1}, -\eta_{1}c_{n+1}-\eta_{2}g_{n+1}, \\ \eta_{2}(b_{n+1}+f_{n+1}), -\eta_{2}(c_{n+1}+g_{n+1})]^{T}, \end{equation} $

(3.12)

So we obtain that equation hierarchy (2.19) possess the Hamiltonian structure

$ \begin{equation} \bar{U}_{t_n}=\bar{K}_{n}(\bar{u})=\bar{J}\frac{\delta H_{n}}{\delta \bar{u}}, ~~~~n\geq 0, \end{equation} $

(3.13)

where the Hamiltonian operator and the Hamiltonian functions are given by

$ \begin{equation} \bar{J}=\frac{1}{\eta_{1}-\eta_{2}}\begin{bmatrix} 0&-1&0&1 \\ 1&0&-1&0 \\ 0&1&0&-\frac{\eta_{1}}{\eta_{2}} \\ -1&0&\frac{\eta_{1}}{\eta_{2}}&0 \\ \end{bmatrix}, \end{equation} $

(3.14)

$ \begin{equation} H_{n}=\int\frac{\eta_{1}a_{n+2}+\eta_{2}e_{n+2}}{n+1}dx, ~~~~n\geq 0. \end{equation} $

(3.15)

With the help of (2.14), we can see a recursion relation

$ \begin{equation} \bar{L}\frac{\delta H_{n}}{\delta\bar{u}}=\frac{\delta H_{n+1}}{\delta\bar{u}}, ~~~~n\geq 0. \end{equation} $

(3.16)

with

$ \begin{equation} \bar{L}=\begin{bmatrix} L&L_{a}\\ 0&L+L_{a} \\ \end{bmatrix}, \end{equation} $

(3.17)

where $L$ is given by (2.8) and

$ \begin{equation} L_{a}=\begin{bmatrix} (p_{1}+q)\partial^{-1}p_2+p_1\partial^{-1}r &(p_{1}+q)\partial^{-1}p_1+p_1\partial^{-1}q \\ -(p_{2}+r)\partial^{-1}p_2-p_2\partial^{-1}r &-(p_{2}+r)\partial^{-1}p_1-p_2\partial^{-1}q \\ \end{bmatrix}. \end{equation} $

(3.18)

Up to now, we have already obtained Hamiltonian structure (42) of integrable couplings of the Guo hierarchy, we must point out that by changing the nonlinear coupling terms of equations, more nonlinear integrable couplings with physical meaning can be obtained. So, with the help of this method, more meaningful results of other integrable hierarchies can be generated.

In what follows, we will construct conservation laws for the nonlinear integrable couplings of the Guo hierarchy. For the coupled spectral problem of Guo hierarchy

$ \begin{equation} \bar U=\begin{bmatrix} 1/\lambda&q+r &0& 0 \\ q-r& -1/\lambda&0&0 \\ 0& p_{1}+p_2& 2/\lambda &q+r+ p_{1}+p_2\\ p_{1}-p_2&0& q-r+ p_{1}-p_2&-2/\lambda \\ \end{bmatrix}, \end{equation} $

(4.1)

we introduce the variables

$ \begin{equation} M=\frac{\phi_{2}}{\phi_{1}}, N=\frac{\phi_{3}}{\phi_{1}}, K=\frac{\phi_{4}}{\phi_{1}}. \end{equation} $

(4.2)

From (4.1), we have

$ \begin{eqnarray} \begin{array}{ll} M_{x}=\frac{1}{2}[q-r-2M\lambda^{-1}-(q+r)M^{2}], \\ N_{x}=\frac{1}{2}[N\lambda^{-1}+(p_1+p_2)M+(q+r+p_1+p_2)K-(q+r)MN], \\ K_{x}=\frac{1}{2}[p_1-p_2+(q-r+p_1-p_2)N-3K\lambda^{-1}-(q+r)MK].\\ \end{array}% \end{eqnarray} $

(4.3)

We expand $M, N, K$ in powers of $\lambda$ as follows

$ \begin{eqnarray*} M=\sum^{\infty}_{j=1}m_{j}\lambda^{j}, N=\sum^{\infty}_{j=1}n_{j}\lambda^{j}, K=\sum^{\infty}_{j=1}k_{j}\lambda^{j}, \end{eqnarray*} $

(4.4)

where $p(m_{j})=0, p(n_{j})=0, p(k_{j})=1$, Substituting (4.4) into (4.3) and comparing the coefficients of the same power of $\lambda$, we obtain

$ \begin{eqnarray}\nonumber \begin{array}{ll} m_1=\frac{1}{2}(q-r), n_1=0, k_1=\frac{1}{3}(p_1-p_2), m_2=-\frac{1}{2}(q-r)_{x}, \\n_2=-\frac{1}{3}(q+r+p_1+p_2)(p_1-p_2)-\frac{1}{2}(p_1+p_2)(q-r),\\ k_2=-\frac{2}{9}(p_1-p_2)_{x}, \\m_3=\frac{1}{2}(q-r)_{xx}-\frac{1}{8}(q+r)(q-r)^{2}, \\ n_3=-\frac{2}{3}(q+r+p_1+p_2)_x(p_1-p_2)-\frac{4}{9}(q+r+p_1+p_2)(p_1-p_2)_x-\\ (p_1+p_2)_x(q-r)\\ \quad\quad \quad-\frac{1}{2}(p_1+p_2)_x(q-r)_x, \\ k_3=\frac{4}{27}(p_1-p_2)_{xx}-\frac{1}{9}(q+r+p_1+p_2)(q-r+p_1-p_2)(p_1-p_2)_x-\\ \frac{1}{6}(q-r+p_1-p_2)\\ \quad\quad \quad\times(p_1+p_2)(q-r)-\frac{1}{18}(p_1-p_2)(q^{2}-r^{2}), \cdots. \end{array}% \end{eqnarray} $

and a recursion formula for $m_j, n_j, k_j, $

$ \begin{eqnarray*} &&m_{j+1}=-m_{j, x}-\frac{1}{2}(q+r)\sum^{j-1}_{l=1}m_lm_{j-l}, \\ &&n_{j+1}=2n_{j, x}-\frac{1}{2}(q+r+p_1+p_2)k_j-\frac{1}{2}(p_1+p_2)m_j-\frac{1}{2}(q+r)\sum^{j-1}_{l=1}m_ln_{j-l}, \\ &&k_{j+1}=-\frac{2}{3}k_{j, x}+\frac{1}{3}(q-r+p_1-p_2)n_j+\frac{1}{3}(q+r)\sum^{j-1}_{l=1}k_lm_{j-l}. \end{eqnarray*} $

(4.5)

Because of

$ \begin{eqnarray} \left\{% \begin{array}{ll} \frac{\partial}{\partial t}[\lambda^{-1}+(q+r)M]=\frac{\partial}{\partial x}[a+(b+c)M], \\ \frac{\partial}{\partial t}[2N\lambda^{-1}+(p_1+p_2)M+(q+r+p_1+p_2)K]\\ \quad=\frac{\partial}{\partial x}[e+(f+g)M+(a+c)N+(b+c+f+g)K], \\ \end{array}% \right. \end{eqnarray} $

(4.6)

where

$ \begin{eqnarray} \begin{array}{ll} a=\lambda^{-2}+\frac{1}{2}(r^{2}-q^{2}), b=q\lambda^{-1}+r_x, c=r\lambda^{-1}+q_x, \\ e=\lambda^{-2}+\frac{1}{2}(r^{2}-q^{2})-2qp_1+2rp_2-p_1^{2}+p_2^{2}, \\ f=(q+2p_1)\lambda^{-1}+r_x+p_{2, x}, g=(r+2p_2)\lambda^{-1}+q_x+p_{1, x}. \end{array}% \end{eqnarray} $

(4.7)

Assume that

$ \begin{eqnarray} \begin{array}{ll} \sigma=\lambda^{-1}+(q+r)M, \\ \theta=a+(b+c)M, \\ \rho=2N\lambda^{-1}+(p_1+p_2)M+(q+r+p_1+p_2)K, \\ \delta=e+(f+g)M+(a+c)N+(b+c+f+g)K.\\ \end{array}% \end{eqnarray} $

(4.8)

Then the (4.6) can be written as $\sigma_{t}=\theta_{x}, \rho_{t}=\delta_{x}$, which are the right form of conservation laws. We expand $\sigma, \theta, \rho$ and $\delta$ as series in powers of $\lambda$ with the coefficients, which are called conserved densities and currents respectively

$ \begin{eqnarray*} &&\sigma=\lambda^{-1}+\sum^{\infty}_{j=1}\sigma_{j}\lambda^{j}, \theta=\lambda^{-2}+\sum^{\infty}_{j=1}\theta_{j}\lambda^{j}, \rho=\sum^{\infty}_{j=1}\rho_{j}\lambda^{j}, \\ &&\delta=\lambda^{-2}-\frac{1}{3}(q^{2}-r^{2})-\frac{2}{3}(p_1^{2}-p_2^{2})+2(p_2r-p_1q) +\sum^{\infty}_{j=1}\delta_{j}\lambda^{j}. \end{eqnarray*} $

(4.9)

The first few conserved densities and currents are read

$ \begin{eqnarray}\nonumber \begin{array}{ll} \sigma_{1}=\frac{1}{2}(q^{2}-r^{2}), \theta_{1}=qr_{x}-q_{x}r, \rho_{1}=\\ -\frac{1}{2}(p_1+p_2)(q-r)-\frac{1}{3}(p_1-p_2)(q+r+p_1+p_2), \\\delta_{1}=-\frac{4}{3}(p_1-p_2)_{x}(q+r+p_1+p_2)+\frac{1}{2}(q^{2}-r^{2})_x\\ -\frac{2}{3}(q+r)_x(p_1-p_2)-(p_1+p_2)_x(p_1-p_2)\\ \quad\quad\quad-\frac{3}{2}(p_1+p_2)_x(q-r), \\ \sigma_{2}=-\frac{1}{2}(q+r)(q-r)_x, \theta_{2}=\frac{1}{2}[(q+r)(q-r)_{xx}-(q_x^{2}-r_x^{2})-\frac{1}{4}(q^{2}-r^{2})], \\\rho_{2}=-\frac{3}{2}(p_1+p_2)(q-r)_x-\frac{10}{9}(p_1-p_2)_x(q+r+p_1+p_2)\\ -\frac{4}{3}(p_1-p_2)(q+r+p_1+p_2)_x -2(p_1+p_2)_x(q-r). \end{array}% \end{eqnarray} $

The recursion relations for $\sigma_j, \theta_j, \rho_j$ and $\delta_j$ are

$ \begin{eqnarray} \begin{array}{ll} \sigma_{j}=(q+r)m_j, \\ \theta_{j}=(q+r)m_{j+1}+(q+r)_xm_j, \\ \rho_{j}=(p_1+p_2)m_j+2n_{j+1}+(q+r+p_1+p_2)k_{j}, \\ \delta_{j}=(q+r+2p_1+2p_2)m_{j+1}+(q+r+p_1+p_2)_xm_{j}+2n_{j+2}+(r^{2}-q^{2}-2qp_1\\ \quad\quad\quad+2rp_2 -p_1^{2}+p_2^{2})n_j+2(q+r+p_1+p_2)k_{j+1}+(2r+2q+p_1+p_2)_xk_j.\\ \end{array}% \end{eqnarray} $

(4.10)

where $m_j, n_j$ and $k_j$ can be calculated from (4.5). The infinite conservation laws of nonlinear integrable couplings (2.19) can be easily obtained in (4.2)-(4.10), respectively.

By changing the nonlinear coupling spectral matrix of equations, more nonlinear integrable couplings and infinite conservation laws with physical meaning can be obtained. So, with the help of this method, more meaningful results of other integrable hierarchies can be generated.