# Low-Rank Tensor Decomposition over Finite Fields

Jason Yang

###### Abstract

We show that finding rank-1, rank-2, and rank-3 decompositions of a 3D tensor
over a fixed finite field can be done in polynomial time. However, if some
cells in the tensor are allowed to have arbitrary values, then rank-2 is NP-
hard over the integers modulo 2. We also explore rank-1 decomposition of a 3D
tensor and of a matrix where some cells are allowed to have arbitrary values.

##  1 Introduction

Given a $p\times q\times s$ tensor (multidimensional array) $T$, a rank-$R$
decomposition is a collection of variables $A_{a,b},B_{c,d},C_{e,f}$ such that

$T_{i,j,k}=\sum_{r=0}^{R-1}A_{r,i}B_{r,j}C_{r,k}\ \forall i,j,k.$

This equation is also written as $T=\sum_{r=0}^{R-1}A_{r,:}\times
B_{r,:}\times C_{r,:}$, where $M_{r,:}$ denotes the slice $[M_{r,j}]_{j}$ and
$\times$ denotes the outer product: for arbitrary tensors $X$ and $Y$,
$(X\times Y)_{i_{0},\dots i_{k-1},j_{0},\dots j_{\ell-1}}:=X_{i_{0},\dots
i_{k-1}}Y_{j_{0},\dots j_{\ell-1}}$.

Tensor decomposition is at the heart of fast matrix multiplication, a problem
that is the primary bottleneck of numerous linear algebra and graph
combinatorics algorithms, such as matrix inversion and triangle detection. All
asymptotically fast algorithms for matrix multiplication use a divide-and-
conquer scheme, and finding an efficient scheme is equivalent to finding a
decomposition of a certain tensor with low rank [2].

While significant progress has been made on this problem over the decades,
most of the asymptotically fastest algorithms discovered use indirect
constructions and convergence arguments that result in enormously large
divide-and-conquer schemes, making the constant factors for these algorithms
far too big for them to be useful in practice. For this reason, there has also
been interest in finding explicit decompositions for the tensors of small
divide-and-conquer schemes. With the introduction of AlphaTensor [3] and flip
graphs [5] just over a year ago, several new low-rank decompositions have been
discovered, especially over finite fields of characteristic 2 [1].

Decompositions using only integers are especially preferred, since they avoid
having to use floating point arithmetic, and most decompositions of small
divide-and-conquer schemes that have been found have happened to use only
integers. This makes it useful to search for decompositions over finite fields
(e.g., integers modulo 2), as doing so makes the search space finite and
allows one to rule out large classes of integer solutions.

It is known that deciding whether a given tensor has a decomposition of rank
$\leq R$ is NP-hard over finite fields via a reduction from 3SAT [4]. However,
since we only want low-rank decompositions, we are more interested in hardness
results for small fixed $R$.

###  1.1 Results

Our main result in this report is a polynomial-time algorithm for 3D tensor
decomposition with rank $\leq 3$. For simplicity of asymptotic analysis, we
focus on $n\times n\times n$ tensors $T$, although our algorithm easily
generalizes to non-cube tensors.

######  Theorem 1.

For $R\leq 3$, finding a rank-$R$ decomposition over a finite field
$\mathbb{F}$ of an $n\times n\times n$ tensor, or determining that it does not
exist, can be done in $O(f(|\mathbb{F}|)n^{3})$ time, for superexponential
$f$.

However, when allowing certain cells of the target tensor to have arbitrary
value (“wildcards”), we have a hardness result over the integers mod 2 via a
reduction from Not-All-Equal 3SAT (NAE-3SAT):

######  Theorem 2.

If some cells of a $p\times q\times s$ tensor $T$ are allowed to have
arbitrary value, then determining if $T$ has a rank-2 decomposition over
$\mathbb{Z}/2\mathbb{Z}$ is NP-hard.

We also sketch out polynomial-time algorithms for rank-1 decomposition with
wildcards over a 3D tensor and a matrix, over an arbitrary fixed finite field.

##  2 Without wildcards

In this section we describe our algorithm for Theorem 1. The inspiration for
this algorithm comes from analyzing rank-$\leq 2$ decomposition over the
integers mod 2.

Rank 1 is easy since $abc=1$ forces $a=1$, $b=1$, and $c=1$, so $T=a\times
b\times c$ if and only if $(a_{i}=\lor_{j,k}T_{i,j,k},\
b_{j}=\lor_{i,k}T_{i,j,k},\ c_{k}=\lor_{i,j}T_{i,j,k})$ is a solution.

Rank 2 at first glance is not clear, but we notice that if $T=a_{0}\times
b_{0}\times c_{0}+a_{1}\times b_{1}\times c_{1}$, then all slices $T_{i,:,:}$
are of the form $u(b_{0}\times c_{0})+v(b_{1}\times c_{1})$ for
$u,v\in\mathbb{Z}/2\mathbb{Z}$, of which there are at most four possibilities
(including the all-zero matrix). These slices also have rank at most 2 (in
fact their ranks are even more restrictive than that), suggesting that one
could do exhaustive casework over the ranks of these slices.

This general idea motivates our main algorithm, which generalizes to rank-3.
Suppose we want to find a rank-$R$ decomposition of $T$ as
$T_{i,j,k}=\sum_{r=0}^{R-1}A_{r,i}B_{r,j}C_{r,k}\ \forall i,j,k$ over the
finite field $\mathbb{F}$. Our idea is to extract distinct nonzero slices
$T_{i,:,:}$, then fix $A$ and solve for $B,C$.

To show why we only need to consider distinct $T_{i,:,:}$, let $I$ be some
list of distinct indices s.t. $\left\\{T_{i,:,:}\right\\}_{i\in
I}=\left\\{T_{i,:,:}\right\\}_{0\leq i<n}\setminus O$ (where $O$ denotes the
$n\times n$ all-zero matrix), and let
$T^{\prime}:=\begin{bmatrix}T_{i,:,:}\end{bmatrix}_{i\in I}$. If $|I|=0$ then
$T$ and $T^{\prime}$ are both all-zero, which is a trivial case. Otherwise, we
can construct an indexing function $f$ s.t. for all $i$ where $T_{i,:,:}\neq
O$, $0\leq f(i)<|I|$ and $T_{i,:,:}=T_{I_{f(i)},:,:}$; then if $T^{\prime}$
has a rank-$R$ decomposition $\sum_{r=0}^{R-1}U^{\prime}_{r,:}\times
V_{r,:}\times W_{r,:}$ we can construct
$U_{r,i}=\left(\begin{cases}0&T_{i,:,:}=O\\\
U^{\prime}_{r,f(i)}&\textrm{else}\end{cases}\right)$ s.t.
$T=\sum_{r}U_{r,:}\times V_{r,:}\times W_{r,:}$. Conversely, if
$T=\sum_{r}A_{r,:}\times B_{r,:}\times C_{r,:}$, then
$T^{\prime}=\sum_{r}\left[A_{r,i}\right]_{i\in I}\times B_{r,:}\times
C_{r,:}$.

Thus, WLOG we only need to solve $T_{i,:,:}=\sum_{r}A_{r,i}(B_{r,:}\times
C_{r,:})\ \forall i\in I$, which we can rewrite into block matrix form as
$\begin{bmatrix}\vdots\\\ T_{i,:,:}\\\ \vdots\end{bmatrix}_{i\in
I}=A^{\intercal}\begin{bmatrix}\vdots\\\ B_{r,:}\times C_{r,:}\\\
\vdots\end{bmatrix}_{r}$, where we implicitly redefine $A$ to only include
column indices $i\in I$.

If $T$ has a rank-$R$ decomposition, there can be at most $|\mathbb{F}|^{R}$
distinct slices, and thus $|\mathbb{F}|^{R}-1$ distinct nonzero slices, so
having $|I|\leq|\mathbb{F}|^{R}-1$ is a necessary condition for the existence
of a solution. This means there are at most
$|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}$ possible $A$, so we avoid exponential
dependence on $n$.

If $|I|=0$, then $T$ must be all-zero, so it trivially has a rank-$R$
decomposition (e.g., with all elements set to 0). Otherwise, for each possible
$A$, we can do Gaussian elimination on
$\left[\begin{array}[]{c|c}A^{\intercal}&\begin{bmatrix}\vdots\\\ T_{i,:,:}\\\
\vdots\end{bmatrix}_{i\in I}\end{array}\right]$ to obtain

$M:=\left[\begin{array}[]{c|c}\Xi&\begin{bmatrix}\vdots\\\ D_{i}\\\
\vdots\end{bmatrix}_{0\leq i<|I|}\end{array}\right]$ for
$\Xi:=\mathrm{rref}\left(A^{\intercal}\right)$ and $n\times n$ matrices
$D_{i}$. If there is some row $r$ s.t. $\Xi_{r,:}$ is all-zero but $D_{r}$ is
not, then there is no solution. Otherwise, since all all-zero rows of a
reduced row-echelon form (“rref”) matrix must be at the bottom, we can remove
everything except the first $N$ rows of $M$, where $N$ is the number of
nonzero rows of $\Xi$, to get
$\left[\begin{array}[]{c|c}E&\begin{bmatrix}\vdots\\\ D_{i}\\\
\vdots\end{bmatrix}_{0\leq i<N}\end{array}\right]$, where $E:=\Xi_{:N,:}$.

At this point, we analyze all possible forms of the truncated rref matrix $E$
with $R$ columns and between $1$ and $R$ rows, for $R\leq 3$:

For $R=1$, we must have $E=\begin{bmatrix}x\end{bmatrix}$, so we simply check
whether $\mathrm{rank}\left(D_{0}\right)\leq\begin{cases}1&x\neq 0\\\
0&\textrm{else}\end{cases}$.

For $R=2$, we either have $E=\begin{bmatrix}1\\\ &1\end{bmatrix}$, in which
case we check whether $\mathrm{rank}\left(D_{0}\right)\leq 1\ \land\
\mathrm{rank}\left(D_{1}\right)\leq 1$, or
$E=\begin{bmatrix}1&x\end{bmatrix},\begin{bmatrix}0&1\end{bmatrix}$, in which
case we check whether $\mathrm{rank}\left(D_{0}\right)\leq(\\#\textrm{ of
nonzeros in }E_{0,:})$.

For $R=3$, there are a few more cases, of which one is nontrivial:

  * •

$E=\begin{bmatrix}1&x&y\end{bmatrix},\begin{bmatrix}0&1&x\end{bmatrix},\begin{bmatrix}0&0&1\end{bmatrix}$:
check whether $\mathrm{rank}\left(D_{0}\right)\leq(\\#\textrm{ of nonzeros in
}E_{0,:})$;

  * •

$E=\begin{bmatrix}1\\\ &1\\\ &&1\end{bmatrix}$: check whether
$\mathrm{rank}\left(D_{i}\right)\leq 1\ \forall i=0,1,2$;

  * •

$E=\begin{bmatrix}1&x\\\ &&1\end{bmatrix},\begin{bmatrix}0&1\\\
&&1\end{bmatrix}$: check whether
$\mathrm{rank}\left(D_{i}\right)\leq(\\#\textrm{ of nonzeros in }E_{i,:})\
\forall i=0,1$;

  * •

$E=\begin{bmatrix}1&&x\\\ &1&y\end{bmatrix}$: this is the most difficult case.

If $x=0$ or $y=0$, we can do the same thing as in $\begin{bmatrix}1&x\\\
&&1\end{bmatrix}$ or $\begin{bmatrix}0&1\\\ &&1\end{bmatrix}$, as the sets of
nonzero terms of different rows in $E$ are pairwise disjoint.

Otherwise, check whether $\exists i\textrm{ s.t.
}\mathrm{rank}\left(D_{i}\right)=2$ (notice the $=$ instead of $\leq$ this
time):

    * –

if such $i$ exists, enumerate all rank-2 decompositions of $D_{i}=u_{0}\times
v_{0}+u_{1}\times v_{1}$, and for each decomposition, check whether
$\mathrm{rank}\left(D_{j}-u_{1}\times v_{1}\right)\leq 1$ for the other index
$j\neq i$;

    * –

otherwise, $\mathrm{rank}\left(D_{0}\right)\leq 1$ and
$\mathrm{rank}\left(D_{1}\right)\leq 1$, so we can force
$B_{2,:},C_{2,:}=\vec{0}$ and solve for $D_{0}=B_{0,:}\times C_{0,:},\
D_{1}=B_{1,:}\times C_{1,:}$;

In contrast to general tensors, rank decomposition of a $n\times n$ matrix
(with the decomposition having smallest possible rank) can be computed in
$O(n^{3})$ time via Gaussian elimination, so if all relevant checks in our
algorithm pass, it is easy to construct a solution for $B,C$.

It remains to count how many rank-2 decompositions there are for a rank-2
matrix. Fortunately we have such a bound for arbitrary matrix rank which does
not have exponential dependence on $n$:

######  Theorem 3.

Consider an $n\times m$ matrix $M$ over an arbitrary field with rank
$r:=\mathrm{rank}\left(M\right)>0$. Let
$F_{0}=\mathrm{rref}\left(M\right)_{:r,:}$ and $C_{0}$ satisfy $C_{0}F_{0}=M$.
Then for any $n\times r$ matrix $C$ and $r\times m$ matrix $F$ s.t. $CF=M$,
there is some invertible $r\times r$ $X$ s.t. $C=C_{0}X^{-1}$ and $F=XF_{0}$.

###### Proof.

First we note that $C_{0}$ exists and can be constructed via Gaussian
elimination by maintaining the net linear transformation of row operations in
a separate matrix (and then truncating to the first $r$ columns).
Additionally, $C_{0}$ is unique, because $\mathrm{rank}\left(F_{0}\right)=r$
and $r\leq\min\left(n,m\right)$ force the rows of $F_{0}$ to be linearly
independent.

We have
$\mathrm{rowspan}\left(M\right)=\mathrm{rowspan}\left(CF\right)\subseteq\mathrm{rowspan}\left(F\right)$
and $\mathrm{rank}\left(F\right)\leq r$, so
$\mathrm{rowspan}\left(M\right)=\mathrm{rowspan}\left(F\right)$. Since row
reduction and removing all-zero rows preserves row-span,
$\mathrm{rowspan}\left(F_{0}\right)=\mathrm{rowspan}\left(\mathrm{rref}\left(F\right)\right)$.

We recall a property of rref matrices. Let $W$ be an arbitrary rref matrix
with no all-zero rows. By definition of rref, all terms directly above and
below any leading 1 must be 0, so in any linear combination of rows of $W$,
leading 1s cannot be canceled out with any other terms. Thus, denoting $l_{i}$
as the column position of the leading 1 at row $i$, the position of the
leftmost nonzero element of the row vector $\vec{v}:=\sum_{i}c_{i}W_{i,:}$ for
scalars $c_{i}$ is $l_{k}$ for $k=\min{\left\\{i:c_{i}\neq 0\right\\}}$.
Furthermore, $\vec{v}_{l_{i}}=c_{i}$, so given an arbitrary linear combination
$\vec{v}$ of the rows of $W$ it is possible to determine the sequence of
coefficients $c_{i}$ used to obtain the linear combination. Since the rows of
$W$ are linearly independent, there is exactly one such sequence of
coefficients $c_{i}$.

Since $\mathrm{rowspan}\left(\mathrm{rref}\left(F\right)\right)$ must contain
every row of $F_{0}$, it is possible for a row vector in
$\mathrm{rowspan}\left(\mathrm{rref}\left(F\right)\right)$ to have its
leftmost nonzero element be at any leading 1’s column in $F_{0}$, forcing
$\mathrm{rref}\left(F\right)$ to have its set of column indices of leading 1s
be the same as that of $F_{0}$, and thus forcing the leading 1s of
$\mathrm{rref}\left(F\right)$ to be in the same set of locations as
$\mathrm{rref}\left(F_{0}\right)$. Then, since
$\mathrm{rowspan}\left(\mathrm{rref}\left(F\right)\right)$ must contain every
row of $F_{0}$, and each row of $F_{0}$ has a 1 at the column index its own
leading 1 and a 0 at all column indices of all other leading 1s,
$\mathrm{rref}\left(F\right)$ itself must contain each row of $F_{0}$. Because
leading 1s in a rref matrix must appear left-to-right as one goes from the top
to the bottom row, this is enough to force
$\mathrm{rref}\left(F\right)=F_{0}$.

By Gaussian elimination, there is some invertible $r\times r$ matrix $X$ s.t.
$F=X\mathrm{rref}\left(F\right)=XF_{0}$, so $M=CF=CXF_{0}$; since the equation
$M=AF_{0}$ has exactly one solution for $A$, $CX$ is forced to equal $C_{0}$,
so $C=C_{0}X^{-1}$. ∎

The upshot of this result is that over a finite field $\mathbb{F}$, there are
exactly $|\mathrm{GL}_{r}\left(\mathbb{F}\right)|=O(|\mathbb{F}|^{r^{2}})$
many rank-$r$ decompositions of a rank-$r$ matrix, and this quantity is
independent of the number of rows or columns the matrix has. One can generate
all of these decompositions by first constructing $(C_{0},F_{0})$ via Gaussian
elimination in $O(n^{3})$ time, then iterate over every $r\times r$ matrix
$X$, check whether $X$ is invertible in $O(r^{3})$ time (e.g., by checking if
$\det\left(X\right)\neq 0$), then compute $(C,F)=(C_{0}X^{-1},F_{0}X^{-1})$ in
$O((n+m)r^{2})$ time. Total running time is then
$O(n^{3}+|\mathbb{F}|^{r^{2}}(n+m)r^{5})$. If a list of invertible $r\times r$
matrices and their inverses are precomputed, total running time is
$O(n^{3}+|\mathbb{F}|^{r^{2}}(n+m)r^{2})$ instead.

In contrast, the number of rank-$(r+1)$ decompositions of $M$ might already be
exponential in the dimensions of $M$111Consider the case where every element
in $M$ is 1..

###  2.1 Runtime

In our algorithm, we first create the set of all distinct nonzero slices
$T_{i,:,:}$, which takes $O(n^{3})$ time.

Then we enumerate at most $|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}$ possible $A$;
for each $A$, we perform Gaussian elimination on the block matrix $M$, which
has dimensions $(\leq|\mathbb{F}|^{R}-1)\times R$, and each element of $M$ is
itself a $n\times n$ matrix, so we incur
$O\left(\left(|\mathbb{F}|^{R}-1\right)R^{2}n^{2}\right)=O(|\mathbb{F}|^{R}R^{2}n^{2})$
runtime222We can compute the row-reduced echelon form of a $n\times m$ matrix
in $O(nm^{2})$ many operations, using the following algorithm:  • initialize
$r$ to 0; • for $j$ from $0$ to $m-1$ inclusive: – find some row index $i\geq
r$ s.t. the element at row $r$ and column $j$ is nonzero; – swap rows $i$ and
$r$; – eliminate the term in the $j$-th column for each row $r+1,\dots n-1$; –
increment $r$; .

We then check whether there is some row $r$ s.t.
$\mathrm{rref}\left(A^{T}\right)_{r,:}=\vec{0}$ and $D_{r}\neq O$, which takes
$O((|\mathbb{F}|^{R}-1)(n^{2}+R))$ time. If no such row exists, then we
proceed to the rref case analysis, where we check the ranks of one or more
matrices $D_{i}$, each of which takes $O(n^{3})$ time via Gaussian
elimination.

If $R=3$ we also have to enumerate all rank-2 decompositions of a rank-2
$n\times n$ matrix $D_{i}$ in the worst case, which takes
$O(n^{3}+|\mathbb{F}|^{4}nR^{2})$ time.

If all relevant checks pass, then constructing a solution $A,B,C$ consists of
computing matrix rank factorizations of $\leq R$ many $D_{i}$, each of which
takes $O(n^{3})$ time.

Since there are at most $R$ many $D_{i}$ we analyze once we get to the rref
casework phase, our final runtime is
$O\left(n^{3}+\underbrace{|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}}_{\textrm{enumerating
}A}\left(\underbrace{|\mathbb{F}|^{R}R^{2}n^{2}}_{\textrm{row-reducing
}M}+\underbrace{(|\mathbb{F}|^{R}-1)(n^{2}+R)}_{\textrm{checking rows of
}M}+Rn^{3}+\left(\begin{cases}n^{3}+|\mathbb{F}|^{4}nR^{2}&R=3\\\ 0&R\leq
2\end{cases}\right)\right)+Rn^{3}\right)\\\
=O\left(|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}\left(|\mathbb{F}|^{R}R^{2}n^{2}+Rn^{3}+\left(\begin{cases}n^{3}+|\mathbb{F}|^{4}nR^{2}&R=3\\\
0&R\leq 2\end{cases}\right)\right)\right)\\\
=\begin{cases}O\left(|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}\left(|\mathbb{F}|^{R}R^{2}n^{2}+(R+1)n^{3}+|\mathbb{F}|^{4}nR^{2}\right)\right)&R=3\\\
O\left(|\mathbb{F}|^{(|\mathbb{F}|^{R}-1)R}\left(|\mathbb{F}|^{R}R^{2}n^{2}+Rn^{3}\right)\right)&R\leq
2\end{cases}$.

Thus, for any fixed $\mathbb{F}$ and fixed $R$ our algorithm runs in
$O(n^{3})$ time, with hidden constant superexponential in $|\mathbb{F}|$ and
doubly exponential in $R$; for $\mathbb{F}=\mathbb{Z}/2\mathbb{Z}$ and $R=3$,
the coefficient on the $n^{3}$ term is $2^{3(2^{3}-1)}\cdot(3+1)=4\cdot
2^{21}\approx 8.4\times 10^{6}$.

###  2.2 Rank 1

It is worth mentioning a much simpler and faster algorithm for rank-1
decomposition: if $T=a\times b\times c$, then all nonzero slices $T_{i,:,:}$
must have rank exactly 1 and all be scalar multiples of the same matrix (which
is $b\times c$). Thus we can collect all nonzero slices $T_{i,:,:}$, scale
each of them s.t. their first element in row-major order is 1, then check
whether all resulting slices are equal and have rank 1. One can see that this
algorithm runs in $O(n^{3})$ time with constant factor independent of
$|\mathbb{F}|$.

###  2.3 Beyond rank 3

It is not clear how to make our algorithm run in $\mathrm{poly}(n)$ time for
$R=4$. Consider the case where $E=\begin{bmatrix}1&&w&x\\\
&1&y&z\end{bmatrix}$ and
$\mathrm{rank}\left(D_{0}\right)=2,\mathrm{rank}\left(D_{1}\right)=2$: unlike
when $E=\begin{bmatrix}1&&x\\\ &1&y\end{bmatrix}$, we cannot simply force
$B_{3,:},C_{3,:}=\vec{0}$, as the system of equations might become too
restrictive. One might also notice that $D_{0}+D_{1}=B_{0,:}\times
C_{0,:}+B_{1,:}\times C_{1,:}$, so $D_{0}+D_{1}$ has rank $\leq 2$; however,
if it happens to have rank 1, then there still might be exponentially many
solutions for $(B_{0,:},C_{0,:},B_{1,:},C_{1,:})$, so it is too expensive to
simply iterate over all such solutions and check for feasibility of
$(B_{2,:},C_{2,:},B_{3,:},C_{3,:})$.

##  3 With wildcards

When certain elements in $T$ are wildcards (i.e., allowed to have arbitrary
value), tensor decomposition becomes much more difficult, and we prove rank-2
over $\mathbb{Z}/2\mathbb{Z}$ is NP-hard via a reduction from NAE-3SAT.

Let $\mathcal{I}:=\begin{bmatrix}\begin{bmatrix}1&0\\\
0&0\end{bmatrix}&\begin{bmatrix}0&0\\\ 0&1\end{bmatrix}\end{bmatrix}$ be the
$2\times 2\times 2$ tensor with 1s along the main diagonal and 0s everywhere
else. It can be verified by brute force that the only rank-2 decompositions of
$\mathcal{I}$ are $\begin{bmatrix}1\\\ 0\end{bmatrix}^{\times
3}+\begin{bmatrix}0\\\ 1\end{bmatrix}^{\times 3}$ and $\begin{bmatrix}0\\\
1\end{bmatrix}^{\times 3}+\begin{bmatrix}1\\\ 0\end{bmatrix}^{\times 3}$ (even
when allowing all-zero vectors). We can parameterize both of these solutions
as $\begin{bmatrix}a\\\ \neg a\end{bmatrix}^{\times 3}+\begin{bmatrix}\neg
a\\\ a\end{bmatrix}^{\times 3}$, yielding a variable gadget. By concatenating
$n$ copies of $\mathcal{I}$ diagonally333For an example, the result of
concatenating $\begin{bmatrix}\begin{bmatrix}a_{0}&a_{1}\\\
a_{2}&a_{3}\end{bmatrix}&\begin{bmatrix}a_{4}&a_{5}\\\
a_{6}&a_{7}\end{bmatrix}\end{bmatrix}$ and
$\begin{bmatrix}\begin{bmatrix}b_{0}&b_{1}\\\
b_{2}&b_{3}\end{bmatrix}&\begin{bmatrix}b_{4}&b_{5}\\\
b_{6}&b_{7}\end{bmatrix}\end{bmatrix}$ diagonally is
$\begin{bmatrix}\begin{bmatrix}a_{0}&a_{1}\\\ a_{2}&a_{3}\\\ &&&\\\
&&&\end{bmatrix}&\begin{bmatrix}a_{4}&a_{5}\\\ a_{6}&a_{7}\\\ &&&\\\
&&&\end{bmatrix}&\begin{bmatrix}&&&\\\ &&&\\\ &&&\\\ &&b_{0}&b_{1}\\\
&&b_{2}&b_{3}\end{bmatrix}&\begin{bmatrix}&&&\\\ &&&\\\ &&&\\\
&&b_{4}&b_{5}\\\ &&b_{6}&b_{7}\end{bmatrix}\end{bmatrix},$ where the unmarked
cells can be filled with arbitrary things. to form the target tensor $T$ and
setting all other elements to be wildcards, we force all rank-2 decompositions
of $T$ to be of the form $\begin{bmatrix}\vdots\\\ v_{i}\\\ \neg v_{i}\\\
\vdots\end{bmatrix}^{\times 3}+\begin{bmatrix}\vdots\\\ \neg v_{i}\\\ v_{i}\\\
\vdots\end{bmatrix}^{\times 3}$, for arbitrary variables $v_{i}$.

Since $\mathrm{NAE}\left(a,b,c\right)=\neg(abc\lor(\neg a)(\neg b)(\neg c))$,
and the AND-clauses $abc$ and $(\neg a)(\neg b)(\neg c)$ cannot both be equal
at the same time, $\mathrm{NAE}\left(a,b,c\right)=\neg(abc\oplus(\neg a)(\neg
b)(\neg c))=(abc+(\neg a)(\neg b)(\neg c)=0)$, where $\oplus$ denotes XOR and
is isomorphic to addition mod 2. Because $T$ already forces all elements of
vectors in the second outer product to be negations of vectors in the first
outer product, we can encode $\mathrm{NAE}\left(v_{i},v_{j},v_{k}\right)$ by
setting $T_{2i,2j,2k}=0$; for negative literals, we add 1 to the corresponding
tensor index, e.g., $\mathrm{NAE}\left(v_{i},v_{j},\neg v_{k}\right)$ can be
encoded by setting $T_{2i,2j,2k+1}=0$.

Thus, we convert a NAE-3SAT instance with $n$ variables and $m$ clauses into a
rank-2 decomposition instance on a $2n\times 2n\times 2n$ tensor with $8n+m$
fixed elements and all other elements being wildcards.

Over arbitrary finite fields, there is an issue of scale invariance, i.e.,
each outer product $a\times b\times c$ can be replaced with $(\alpha
a)\times(\beta b)\times(\gamma c)$ for any $\alpha,\beta,\gamma$ s.t.
$\alpha\beta\gamma=1$, allowing any one tensor decomposition to be transformed
into many other tensor decompositions that evaluate to the same tensor. This
scaling freedom makes constructing gadgets less obvious for NP-hardness.

###  3.1 Rank 1 with wildcards

Here we outline a polynomial-time algorithm for rank-1 decomposition with
wildcards over a fixed finite field $\mathbb{F}$.

First, we rewrite the problem as follows: given an arbitrary subset
$S\subseteq\left\\{0,\dots n-1\right\\}^{3}$ and set of elements
$T_{i,j,k}\in\mathbb{F}\ \forall(i,j,k)\in S$, find variables $a_{0\dots
n-1},b_{0\dots n-1},c_{0\dots n-1}$ s.t. $a_{i}b_{j}c_{k}=T_{i,j,k}\
\forall(i,j,k)\in S$.

When $\mathbb{F}=\mathbb{Z}/2\mathbb{Z}$, $abc=1$ implies $(a,b,c)=(1,1,1)$,
so a solution exists if and only if setting

$a_{i}=\lor_{j,k:(i,j,k)\in S}T_{i,j,k},\ b_{j}=\lor_{i,k:(i,j,k)\in
S}T_{i,j,k},\ c_{k}=\lor_{i,j:(i,j,k)\in S}T_{i,j,k}$

yields a solution; this is because this construction produces the minimal set
of elements in $a\times b\times c$ that are forced to be 1s by the values
$T_{i,j,k}$, and this set covers all $(i,j,k)\in S$ that have $T_{i,j,k}=1$.
Thus for $\mathbb{F}=\mathbb{Z}/2\mathbb{Z}$ there is an algorithm for rank-1
wildcard decomposition that runs in $O(n^{3})$ time with low constant factor.

For an arbitrary fixed finite field $\mathbb{F}$, we outline a polynomial-time
algorithm.

Let $S^{\prime}:=\left\\{(i,j,k)\in S:T_{i,j,k}\neq 0\right\\}$ denote all
points with nonzero values and $I:=\left\\{i:\exists(i,j,k)\in
S^{\prime}\right\\}$, $J:=\left\\{j:\exists(i,j,k)\in S^{\prime}\right\\}$,
$K:=\left\\{k:\exists(i,j,k)\in S^{\prime}\right\\}$ denote the sets of
coordinate values of the nonzero points along each dimension. Since
$T_{i,j,k}\neq 0\Rightarrow a_{i},b_{j},c_{k}\neq 0$, we have that $a_{i}\neq
0\ \forall i\in I$, $b_{j}\neq 0\ \forall j\in J$, and $c_{k}\neq 0\ \forall
k\in K$, so we check whether $T_{i,j,k}\neq 0\ \forall(i,j,k)\in(I\times
J\times K)\cap S$. If this check fails, there is no solution.

Otherwise, we proceed with the observation that all the equations
$a_{i}b_{j}c_{k}=T_{i,j,k}$ involve pure products, i.e., expressions of the
form
$\left(\prod_{i}a_{i}^{p_{i}}\right)\left(\prod_{j}b_{j}^{q_{j}}\right)\left(\prod_{k}c_{k}^{r_{k}}\right)$,
for integers $p_{i},q_{j},r_{k}$. If we allow arbitrary multiplication and
division between such expressions, all resulting expressions will still be of
this form, although some $p_{i},q_{j},r_{k}$ might be negative.

This observation suggests that we can solve this problem somewhat like a
linear system of equations: construct a $|S^{\prime}|\times 3n$ integer matrix
$E$ and a $|S^{\prime}|\times 1$ $\mathbb{F}$-vector $W$ s.t. each equation
$a_{i}b_{j}c_{k}=T_{i,j,k}$ for $(i,j,k)\in S^{\prime}$ is encoded as
$\prod_{m}v_{m}^{E_{l,m}}=W_{l}$ for some $l$, where the variables have been
re-indexed as $v_{0,\dots n-1}=a_{0,\dots n-1},\ v_{n,\dots 2n-1}=b_{0,\dots
n-1},\ v_{2n,\dots 3n-1}=c_{0,\dots n-1}$. Since we restrict ourselves to
$S^{\prime}$, each $W_{l}$ is nonzero. During row reductions, elements in $W$
have to be treated differently: when adding between two rows, multiply their
values in $W$; when scaling a row by $s$, raise its value in $W$ to the $s$-th
power.

However, once we reduce $\left[\begin{array}[]{c|c}E&W\end{array}\right]$ to
some normal form elements in $E$ might not be $\pm 1$, so extracting a
solution for the variables $v_{i}$ may involve computing $r$-th roots in
$\mathbb{F}$. An explicit example of this phenomenon is the tensor

$T=\begin{bmatrix}\begin{bmatrix}1&?\\\ ?&1\end{bmatrix},\begin{bmatrix}?&2\\\
1&?\end{bmatrix}\end{bmatrix},$

where ? represents wildcards (equivalently,
$S=\left\\{(0,0,0),(1,1,0),(0,1,1),(1,0,1)\right\\}$, $T_{0,0,0}=1,\
T_{1,1,0}=1,\ T_{0,1,1}=2,\ T_{1,0,1}=1$). Solving this instance is equivalent
to solving the equations

$\left\\{a_{0}b_{0}c_{0}=1,\ a_{1}b_{1}c_{0}=1,\ a_{0}b_{1}c_{1}=2,\
a_{1}b_{0}c_{1}=1\right\\};$

by combining these equations, we get
$\frac{(a_{0}b_{1}c_{1})(a_{1}b_{0}c_{1})}{(a_{0}b_{0}c_{0})(a_{1}b_{1}c_{0})}=\frac{2\cdot
1}{1\cdot 1}\Rightarrow\left(\frac{c_{1}}{c_{0}}\right)^{2}=2$, so solving
this equation requires us to compute a square root of 2.

To resolve this issue, we use the fact that the multiplicative group of any
finite field is cyclic, so if $\mathbb{F}$ is fixed we can find a primitive
element (i.e., multiplicative generator) $\mu$ via brute force and define
$u_{i},w_{i}$ s.t. $v_{i}=\mu^{u_{i}},\ W_{i}=\mu^{w_{i}}$; then our problem
reduces to solving the linear system of equations $Eu=w$ for the vector $u$
over $\mathbb{Z}/(|\mathbb{F}|-1)\mathbb{Z}$.

From [6], solving the matrix-vector equation $xA=b$ over
$\mathbb{Z}/N\mathbb{Z}$ for $n\times m$ $A$ can be done in $O(nm^{\omega-1})$
many operations, where $\omega\leq 2.371552$ [8] is the exponent for matrix
multiplication, with all integers throughout the algorithm having magnitude
$\leq N$. By transposing our system to
$u^{\intercal}E^{\intercal}=w^{\intercal}$ and noting that $E$ has dimensions
$|S^{\prime}|\times 3n$, we see that solving our system takes
$O(n|S^{\prime}|^{\omega-1})$ operations, which is within $O(n^{3\omega-2})$
since $|S^{\prime}|\leq n^{3}$. For fixed $\mathbb{F}$, since all integers
used in this algorithm have magnitude $\leq|\mathbb{F}|-1$, each operation
takes $O(1)$ time, so our algorithm runs in polynomial time. If we account for
the number of wildcards in $T$ in the time complexity, then by letting $k$ be
the number of wildcards, $|S^{\prime}|\leq|S|=n^{3}-k$, so we get
$O(n(n^{3}-k)^{\omega-1})$.

The main bottleneck in this algorithm is finding a primitive element of
$\mathbb{F}$, for which an efficient algorithm is not known [7]. On the other
hand, for $\mathbb{F}=\mathbb{F}_{p^{k}}$, the total input size of a problem
instance would be linear w.r.t. $\log(p)$ and $k$444Every element in
$\mathbb{F}_{p^{k}}$ is a linear combination of $1,\alpha,\dots,\alpha^{k-1}$
with coefficients in $\mathbb{Z}/p\mathbb{Z}$, for some element $\alpha$, so
every element can be specified with a $k$-element vector over
$\mathbb{Z}/p\mathbb{Z}$.. Because of this, we are not sure what the
computational hardness of rank-1 wildcard decomposition is w.r.t. input size
when $p,k$ are not fixed parameters.

###  3.2 Rank-1 with wildcards for matrices

One may wonder whether matrix rank factorization with wildcards is NP-hard.
The issue is that matrix factorizations have an immense degree of freedom, as
$CF=M$ implies $(CX^{-1})(XF)=M$, so constructing gadgets seems very difficult
even for rank $R=2$.

Additionally, for rank $R=1$, we sketch an $O(n^{3})$ algorithm that does not
contain $|\mathbb{F}|$ in its asymptotic running time, in contrast to the
analogous problem for 3D tensors. Again we formulate the problem as follows:
given a subset $S\subseteq\left\\{0,\dots n-1\right\\}^{2}$ and set of
elements $T_{i,j}\in\mathbb{F}\ \forall(i,j)\in S$, find variables $a_{0,\dots
n-1},b_{0,\dots n-1}$ s.t. $a_{i}b_{j}=T_{i,j}$.

Like before, define $S^{\prime}:=\left\\{(i,j)\in S:T_{i,j}\neq 0\right\\}$,
$I:=\left\\{i:(i,j)\in S^{\prime}\right\\}$, $J:=\left\\{j:(i,j)\in
S^{\prime}\right\\}$. Then $a_{i}\neq 0\ \forall i\in I$ and $b_{j}\neq 0\
\forall j\in J$, so we check whether $T_{i,j}\neq 0\ \forall(i,j)\in(I\times
J)\cap S$; if the check fails then there is no solution.

Otherwise, we notice that since every equation $a_{i}b_{j}=T_{i,j}$ constrains
the product of two variables, knowing one variable immediately forces the
value of the other variable. This motivates us to construct a graph $G$ with
vertex set $S^{\prime}$ where there is an edge between every pair of vertices
$(i,j),(i^{\prime},j^{\prime})$ s.t. $i=i^{\prime}$ or $j=j^{\prime}$.

Choose an arbitrary $(i_{0},j_{0})\in S^{\prime}$: WLOG we can set
$a_{i_{0}}=1,b_{j_{0}}=T_{i_{0},j_{0}}$ because any solution $(a,b)$ to the
original problem can be scaled to $(\alpha a,\frac{1}{\alpha}b)$ for any
$\alpha\neq 0$. Then do a breadth-first search over $S^{\prime}$ starting from
$(i_{0},j_{0})$: when we visit each new vertex $(i,j)\in S$, we will have
already determined $a_{i}$ and $b_{j}$, so for each neighbor of the form
$(i^{\prime},j)$ we set $a_{i^{\prime}}=T_{i^{\prime},j}/b_{j}$, and for each
neighbor of the form $(i,j^{\prime})$ we set
$b_{j^{\prime}}=T_{i,j^{\prime}}/a_{i}$. Each neighbor that was not already
visited is added to a queue of vertices to visit. Note that because every
$T_{i,j}$ we deal with is nonzero, we will never encounter division by 0 in
this procedure.

If at any point we update a variable $a_{i}$ or $b_{j}$ that already has a
assigned value and set it to a different value, then we have encountered an
inconsistency and there is no solution to the original problem. Otherwise, at
the end of this procedure we will have searched through the entire connected
component of $G$ that contains $(i_{0},j_{0})$, so we do another breadth-first
search starting at an arbitrary unvisited vertex in $G$, and repeat until all
vertices in $G$ have been visited. If no inconsistencies were ever
encountered, we can set all variables $a_{i}$ and $b_{j}$ that were not
already assigned a value to 0, and we have a solution.

Since breadth-first search over a graph with $V$ vertices and $E$ edges takes
$O(V+E)$ time, and $G$ has $O(n^{2})$ vertices and $O(n^{3})$ edges, our
algorithm runs in $O(n^{3})$ time.

The key for correctness is that if two vertices
$(i,j),(i^{\prime},j^{\prime})$ are in different connect components, then
$(i^{\prime},j^{\prime})$ is not adjacent to any vertex $(x,y)$ in the
connected component containing $(i,j)$, so $i^{\prime}\neq x$ and
$j^{\prime}\neq y$. This means that $a_{i^{\prime}}$ and $b_{j^{\prime}}$ are
not directly or indirectly involved in any equations with elements $T_{x,y}$
for each $(x,y)$ in the connected component of $(i,j)$. Thus, after we finish
searching over each connected component, all remaining unassigned variables
$a_{i}$ and $b_{j}$ are completely independent of all previously assigned
variables, so as we search through the next connected component we will never
encounter inconsistencies with variables assigned during a previous connected
component. This independence between connected components also allows us to
assume WLOG that $a_{i_{0}}=1$ and $b_{j_{0}}=T_{i_{0},j_{0}}$ at each
individual connected component.

If we tried generalizing this algorithm to 3D tensors, we would encounter a
problem: two vertices would have to match 2 out of 3 coordinates in order to
be adjacent to each other, because each equation $a_{i}b_{j}c_{k}=T_{i,j,k}$
constrains the product of three variables instead of two, so two variables
have to be known beforehand to force the value of the third one. But then if
$(i,j,k)$ and $(i^{\prime},j^{\prime},k^{\prime})$ occupy different
components, it is still possible for $(i^{\prime},j^{\prime},k^{\prime})$ to
share 1 out of 3 coordinates with a vertex contained in the connected
component containing $(i,j,k)$, so different connected components can still
share common variables and thus interact with each other in a nontrivial
manner.

##  4 Future directions

We conclude with some open questions:

  * •

For what $R\geq 4$ can rank-$R$ tensor decomposition (without wildcards) be
computed in polynomial time over a fixed finite field?

  * •

For $R\leq 3$, how low can the constant factor of rank-$R$ decomposition be
w.r.t. to the size of the finite field and w.r.t. $R$?

  * •

What is the parameterized complexity of rank-2 wildcard decomposition w.r.t.
the number of wildcard elements in the target tensor? On a related note, how
asymptotically low can the number of wildcards in a tensor be while keeping
rank-2 wildcard decomposition NP-hard?

  * •

What is the parameterized complexity of rank-1 wildcard decomposition w.r.t.
the size or the characteristic of the underlying finite field?

  * •

Which of our algorithms be generalized beyond finite fields, such as the
integers?

## Acknowledgments

We thank Prof. Erik Demaine and TAs Josh Brunner, Lily Chung, and Jenny
Diomidova for running MIT 6.5440 Fall 2023, which made this project possible.
We also thank them for giving feedback on an earlier version of this report,
some of which were requests for certain details and open questions to be
addressed, including rank-1 wildcard decomposition.

## References

  * [1] Arai, Y., Ichikawa, Y., & Hukushima, K. (2023). Adaptive Flip Graph Algorithm for Matrix Multiplication. https://arxiv.org/abs/2312.16960
  * [2] Bläser, M. (2013). Fast Matrix Multiplication. Theory of Computing. https://theoryofcomputing.org/articles/gs005/
  * [3] Fawzi, A., Balog, M., Huang, A., Hubert, T., Romera-Paredes, B., Barekatain, M., Novikov, A., Ruiz, F. J. R., Schrittwieser, J., Swirszcz, G., Silver, D., Hassabis, D., & Kohli, P. (2022). Discovering faster matrix multiplication algorithms with reinforcement learning. Nature, 610, 47-53. https://doi.org/10.1038/s41586-022-05172-4
  * [4] Håstad, J. (1990). Tensor rank is NP-complete. Journal of Algorithms, 11(4), 644-654. https://doi.org/10.1016/0196-6774(90)90014-6
  * [5] Kauers, M., & Moosbauer, J. (2022). Flip graphs for matrix multiplication. https://arxiv.org/abs/2212.01175
  * [6] Storjohann, A., & Mulders, T. (1998). Fast Algorithms for for Linear Algebra Modulo $N$. European Symposium on Algorithms, 139-150. https://doi.org/10.1007/3-540-68530-8_12
  * [7] Von Zur Gathen, J., & Shparlinski, I. (1996). Orders of Gauss Periods in Finite Fields. Applicable Algebra in Engineering, Communication and Computing, 9, 15-24. https://doi.org/10.1007%2Fs002000050093
  * [8] Williams, V. V., Xu, Y., Xu, Z., & Zhou, R. (2023). New Bounds for Matrix Multiplication: from Alpha to Omega. https://arxiv.org/abs/2307.07970


††thanks: ∗ Corresponding author<EMAIL_ADDRESS>
# Deformed two-dimensional rogue waves in the $(2+1)$-dimensional Korteweg-de
Vries equation

Yulei Cao1 1 Institute for Advanced Study, Shenzhen University, Shenzhen,
Guangdong 518060, P. R. China
<EMAIL_ADDRESS>Peng-Yan Hu2∗ 2∗ College of Mathematics and
Statistics, Shenzhen University, Shenzhen 518060, China
<EMAIL_ADDRESS>Yi Cheng3 3 School of Mathematical Sciences, USTC, Hefei,
Anhui 230026, P. R. China
<EMAIL_ADDRESS>Jingsong He1 1 Institute for Advanced Study, Shenzhen
University, Shenzhen, Guangdong 518060, P. R. China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Within the $(2+1)$-dimensional Korteweg-de Vries equation framework, new
bilinear B$\rm\ddot{a}$cklund transformation and Lax pair are presented based
on the binary Bell polynomials and gauge transformation. By introducing an
arbitrary function $\phi(y)$, a family of deformed soliton and deformed
breather solutions are presented with the improved Hirota’s bilinear method.
Choosing the appropriate parameters, their interesting dynamic behaviors are
shown in three-dimensional plots. Furthermore, novel rational solutions are
generated by taking the limit of obtained solitons. Additionally, two
dimensional [2D] rogue waves (localized in both space and time) on the soliton
plane are presented, we refer to it as deformed 2D rogue waves. The obtained
deformed 2D rogue waves can be viewed as a 2D analog of the Peregrine soliton
on soliton plane, and its evolution process is analyzed in detail. The
deformed 2D rogue wave solutions are constructed successfully, which are
closely related to the arbitrary function $\phi(y)$. This new idea is also
applicable to other nonlinear systems.

Keywords: 2D KdV equation $\cdot$ Bilinear method $\cdot$ B$\rm\ddot{a}$cklund
transformation $\cdot$ Lax pair $\cdot$ Deformed 2D rogue wave

PACS numbers: 05.45.Yv $\cdot$ 02.30.Jr $\cdot$ 02.30.Ik

##  1 Introduction

Rogue waves (RWs) in the ocean, also known as killer waves, monster waves,
freak waves and extreme waves, are special waves with extremely destructive,
which should be responsible for some maritime disasters. The main feature of
RW is that ”appear from nowhere and disappear without a trace”[1]. Peregrine
was the first one who obtained the RW solution from one-dimensional
systems[2], so such RW solution is also called ”Peregrine soliton”. These RWs
are localized in both space and time in one-dimensional systems[3]-[9], and
its dynamic behaviors are similar to that of lumps in high-dimensional
systems[10]-[14]. RWs in high-dimensional systems are often called line RWs,
which are merely localized in time[15]-[18]. The study of RW has been widely
used in diverse areas of theoretical and applied physics, including plasma
physics[19, 20], Bose-Einstein condensates[21, 22], atmosphere physics[23],
optics and photonics[24, 25, 26] and superfluids[27].

Now, a fundamental problem occurs: Is it possible that other types of RW
solutions exist in nonlinear systems? As well as seeking new exact solution
for nonlinear systems is an open problem and a challenging work. Furthermore,
in integrable system, the Peregrine soliton only appears in the one-
dimensional systems. But, the RWs, which are localized in time and space, have
never been found in the high-dimensional system. Therefore, searching for the
RW solution that is local in time and space of high-dimensional systems has
always been an open problem. Until recently, Guo et al.[28] obtained RW
solutions, which are localized in both space and time, from a two-dimensional
nonlinear Schr$\rm\ddot{o}$dinger model by the even-fold Darboux
transformation. Now the long-standing problem of the construction of two-
dimensional [2D] rogue wave on zero background is addressed[28]. This is a
pioneering work, which inspires us to explore the structure of 2D RW solutions
of other high-dimensional soliton equations. Additionally, RW solutions of
complex nonlinear equations have been widely researched, but the RW solutions
of real nonlinear equations are rarely mentioned, which is another motivation.

Inspired by the above considerations, in this paper, we focus on the
$(2+1)$-dimensional [2D] Korteweg-de Vries [KdV] equation

$\displaystyle U_{t}+U_{xxx}-3(UV)_{x}=0,\qquad U_{x}=V_{y},$ (1)

which was first derived by Boiti-Leon-Manna-Pempinelli with a weak Lax
pair[29], and can be reduced to the celebrated KdV equation if $x=y$.
Equation(1) is an asymmetric part of the Nizhnik-Novikov-Veselov [NNV]
equation[30], thus equation(1) is also called ANNV equation; this equation(1)
can also be derived by using the inner parameter-dependent symmetry constraint
of the KP equation[30], and can be considered as a model of an incompressible
fluid where $U$ and $V$ are the component of velocity[31]. Equation(1) was
also regarded to be a generalization[32] of the results of Hirota and
Satsuma[33] in the $(2+1)$ dimension.

Letting $U=W_{y},V=W_{x}$, then equation(1) reduce to the following single
equation,

$\displaystyle W_{yt}+W_{xxxy}-3W_{xx}W_{y}-3W_{x}W_{xy}=0.$ (2)

This equation(2) known as BLMP equation was widely investigated. The Cauchy
problem, associated with initial data decaying sufficiently rapidly at
infinity, was solved by means of inverse scattering transformation[29]. The
spectral transformation, nonclassical symmetries and Painlev$\acute{e}$
property for equation(2) has been researched in Refs.[29, 34, 35]. And a
series of B$\rm\ddot{a}$cklund transformations[36, 37, 38], Lax pairs[38],
supersymmetric[39], soliton-like[40]-[43], breather[44, 45] and lump-type[46,
47, 48] solutions are derived.

In this letter, we are committed to exploring new B$\rm\ddot{a}$cklund
transformations, Lax pairs and new RW solutions of 2D KdV equation(1) by the
Hirota method. The new RW solution is called deformed 2D RW because its
formula involves an arbitrary function $\phi(y)$ which provides a deformed
background for rogue wave. This arbitrary function $\phi(y)$ comes from the
crucial bilinear transformation which maps $U$ to a bilinear form. In
section2, new B$\rm\ddot{a}$cklund transformations and Lax pair are derived
based on the binary Bell polynomials and linearizing the bilinear equation. In
section3, we obtain deformed kink soliton and deformed breather solutions of
2D KdV equation by means of the Hirota method. In section4, a family of new
rational solutions and deformed 2D RW solutions of the equation(1) are
presented by introducing an arbitrary function $\phi(y)$. The main results of
the paper are summarized and discussed in section5.

##  2 Bilinear B$\rm\ddot{a}$cklund transformation and associated Lax pair

In this section, we mainly focus on the new B$\rm\ddot{a}$cklund
transformation and Lax pair of 2D KdV equation(1). B$\rm\ddot{a}$cklund
transformation and lax pair are recognized as the main characteristics of
integrability, which can be used to obtain exact solutions of nonlinear
systems. We first introduce the following variables transformation

$\displaystyle U=-q_{xy},\quad V=-q_{xx},$ (3)

then 2D KdV equation(1) becomes

$\displaystyle E(q)=q_{yt}+q_{xxxy}+3q_{xx}q_{xy}=0,$ (4)

two different solutions $q=2\ln G$ and $\widetilde{q}=2\ln F$ of equation(1)
are introduced, the corresponding two-field condition is as follows[49]-[53]

$\displaystyle
E(\widetilde{q})-E(q)=(\widetilde{q}-q)_{yt}+(\widetilde{q}-q)_{xxxy}+3(\widetilde{q}_{xx}\widetilde{q}_{xy}-q_{xx}q_{xy})=0.$
(5)

We introduce two new variables

$\displaystyle v=\frac{\widetilde{q}-q}{2}=\ln\frac{F}{G},\quad
w=\frac{\widetilde{q}+q}{2}=\ln FG.$ (6)

Based on the binary Bell polynomials, equation(5) can be rewritten into
$\mathcal{Y}$-polynomials from Bell polynomial theory as follows and rewrite
the equation(5) into the form (see Refs. [51, 52, 53] for details)

$\displaystyle E(\widetilde{q})-E(q)$
$\displaystyle=2[v_{yt}+v_{xxxy}+3w_{xx}v_{xy}+3w_{xy}v_{xx}]$ (7)
$\displaystyle=2[\partial_{y}(\mathcal{Y}_{xxx}+\mathcal{Y}_{t})+3w_{xy}v_{xx}-3w_{xxy}v_{x}-3v_{x}^{2}v_{xy}].$

Through the following constraints

$\displaystyle w_{xy}+v_{x}v_{y}-\lambda v_{x}=0,$ (8)

where $\lambda$ is the constant parameter. Then, the following coupled system
of $\mathcal{Y}$-polynomials of equation(1) is deduced

$\displaystyle\mathcal{Y}_{xy}-\lambda\mathcal{Y}_{x}=0,$ (9)
$\displaystyle\mathcal{Y}_{xxx}+\mathcal{Y}_{t}=0.$

According to the relationship between $\mathcal{Y}$-polynomials[51, 52, 53]
and bilinear operators $D$[54]

$\displaystyle\mathcal{Y}_{n_{1}x_{1},\cdots,n_{l}x_{l}}(v=\ln\frac{F}{G},w=\ln
FG)$ $\displaystyle=$ $\displaystyle(FG)^{-1}D_{x_{1}}^{n_{1}}\cdots
D_{x_{l}}^{n_{l}}F\cdot G,$ (10)

the following bilinear B$\rm\ddot{a}$cklund transformation is derived for the
2D KdV equation(1)

$\displaystyle(D_{x}D_{y}-\lambda D_{x})F\cdot G=0,$ (11)
$\displaystyle(D^{3}_{x}+D_{t})F\cdot G=0.$

Then logarithmic linearization of $\mathcal{Y}$-polynomials under the Hopf-
Cole transformation $v=\ln\psi$, and Bell polynomial formula

$\displaystyle(FG)^{-1}D_{x_{1}}^{n_{1}}\cdots D_{x_{l}}^{n_{l}}F\cdot
G\bigg{|}_{G=e^{\frac{q}{2}},\frac{F}{G}=\psi}$ (12) $\displaystyle=$
$\displaystyle\sum\limits_{r_{1}+\cdots+r_{l}=even}\sum\limits_{r_{1}=0}^{n_{1}}\cdots\sum\limits_{r_{l}=0}^{n_{l}}\prod\limits_{i=1}^{l}\begin{pmatrix}n_{i}\\\
r_{i}\end{pmatrix}P_{r_{1}x_{1}\cdots r_{l}x_{l}}(q)$
$\displaystyle\times\mathcal{Y}_{(n_{1}-r_{1})x_{1},\cdots,(n_{l}-r_{l})x_{l}}(v),$

the $\mathcal{Y}$-polynomials can be written as

$\displaystyle\mathcal{Y}_{n_{1}x_{1},\cdots,n_{l}x_{l}}(v)=\frac{\psi_{n_{1}x_{1},\cdots,n_{l}x_{l}}}{\psi},$
(13)

thus the system(9) is then linearized into a Lax pair

$\displaystyle\psi_{xy}+q_{xy}\psi-\lambda\psi_{x}=0,$ (14)
$\displaystyle\psi_{xxx}+3q_{xx}\psi_{x}+\psi_{t}=0,$

which is equivalent to the Lax pair of equation(1)[36]

$\displaystyle\psi_{xy}-U\psi-\lambda\psi_{x}=0,$ (15)
$\displaystyle\psi_{xxx}-3V\psi_{x}+\psi_{t}=0.$

In order to get a new B$\rm\ddot{a}$cklund transformation, we introduce the
following gauge transformation to the bilinear B$\rm\ddot{a}$cklund
transformation(11)

$\displaystyle F\rightarrow e^{\xi}F,\quad G\rightarrow e^{\eta}G,$ (16)

where $\xi=m_{1}x+n_{1}y+l_{1}t$ and $\eta=m_{2}x+n_{2}y+l_{2}t$. System(11)
is transformed into

$\displaystyle\left[(D_{x}+m_{1}-m_{2})(D_{y}+n_{1}-n_{2})-\lambda(D_{x}+m_{1}-m_{2})\right]F\cdot
G=0,$ (17)
$\displaystyle\left[(D_{x}+m_{1}-m_{2})^{3}+(D_{t}+l_{1}-l_{2})\right]F\cdot
G=0,$

under the following constraints

$\displaystyle m_{1}-m_{2}=\lambda,\quad n_{1}-n_{2}=\lambda,\quad
l_{2}-l_{1}=\lambda^{3}.$ (18)

Then, a new B$\rm\ddot{a}$cklund transformation of equation(1) is derived as
follows

$\displaystyle(D_{x}D_{y}+\lambda D_{y})F\cdot G=0,$ (19)
$\displaystyle(D^{3}_{x}+3\lambda D^{2}_{x}+3\lambda^{2}D_{x}+D_{t})F\cdot
G=0.$

Similarly, the Lax pair of system(25) is as follows

$\displaystyle\psi_{xy}+\lambda\psi_{y}+q_{xy}\psi=0,$ (20)
$\displaystyle\psi_{xxx}+2\lambda\psi_{xx}+3(\lambda^{2}+q_{xx})\psi_{x}+3\lambda
q_{xx}\psi+\psi_{t}=0,$

or equivalently

$\displaystyle(L_{1}-U)\psi=$
$\displaystyle\psi_{xy}+\lambda\psi_{y}-U\psi=0,$ (21)
$\displaystyle(\partial_{t}+L_{2})\psi=$
$\displaystyle\psi_{xxx}+2\lambda\psi_{xx}+3(\lambda^{2}-V)\psi_{x}-3\lambda
V\psi+\psi_{t}=0.$

It is easy to check that the integrability condition
$[L_{1}-U,\partial_{t}+L_{2}]\psi=0$ is satisfied.

Additionally, another form of B$\rm\ddot{a}$cklund transformation and lax pair
of 2D KdV equation(1) are also derived, when taking

$\displaystyle U=-2\partial_{xy}\ln f(x,y,t)+\phi(y),\quad
V=-2\partial_{xx}\ln f(x,y,t),$ (22)

2D KdV equation(1) is then translated into the following bilinear form,

$\displaystyle[D^{3}_{x}D_{y}+D_{y}D_{t}-3\phi(y)D^{2}_{x}]f(x,y,t)\cdot
f(x,y,t)=0,$ (23)

$\phi(y)$ is an arbitrary function of $y$, $f$ is a real function and $D$ is
the Hirota’s bilinear differential operator[54], based on the binary Bell
polynomials[49, 50, 51], equation(1) admits the following B$\rm\ddot{a}$cklund
transformation

$\displaystyle\left[D_{x}D_{y}-\phi(y)\right]F\cdot G=0,$ (24)
$\displaystyle\left[D^{3}_{x}+D_{t}+k\right]F\cdot G=0.$

If using gauge transformation(16) to the above bilinear B$\rm\ddot{a}$cklund
transformation(24), a new B$\rm\ddot{a}$cklund transformation of equation(1)
is derived

$\displaystyle\left[D_{x}D_{y}+\lambda(D_{x}+D_{y})+\lambda^{2}-\phi(y)\right]F\cdot
G=0,$ (25) $\displaystyle\left[D^{3}_{x}+3\lambda
D^{2}_{x}+3\lambda^{2}D_{x}+D_{t}+k\right]F\cdot G=0.$

Then, we derive the corresponding linear system under the Hopf-Cole
transformation(13)

$\displaystyle\psi_{xy}+\lambda(\psi_{x}+\psi_{y})+(\lambda^{2}+q_{xy}-\phi(y))\psi=0,$
(26)
$\displaystyle\psi_{xxx}+3\lambda\psi_{xx}+3(\lambda^{2}+q_{xx})\psi_{x}+\psi_{t}+(3\lambda
q_{xx}+k)\psi=0,$

which is equivalent to the Lax system of 2D KdV equation(1)

$\displaystyle\psi_{xy}+\lambda(\psi_{x}+\psi_{y})+(\lambda^{2}-U)\psi=0,$
(27)
$\displaystyle\psi_{xxx}+3\lambda\psi_{xx}+3(\lambda^{2}-V)\psi_{x}+\psi_{t}+(k-3\lambda
V)\psi=0.$

##  3 Deformed multi-kink soliton and deformed breather solutions of the 2D
KdV equation

Based on the above new bilinear equation(23), some exact solutions of 2D KdV
equation(1), including deformed multi-solitons and deformed breathers, are
generated. First, deformed one-soliton solutions take the forms

$\displaystyle U_{[1]}=-2\partial_{xy}\ln(1+e^{\eta_{1}})+\phi(y),\quad
V_{[1]}=-2\partial_{xx}\ln(1+e^{\eta_{1}}),$ (28)

where
$\eta_{1}=p_{1}x+\Phi(y)q_{1}y+\Omega_{1}t+\eta^{0}_{1},\;\eta^{0}_{1}=0$ and
$\phi(y)=\frac{d\Phi(y)}{dy}$. Substituting (28) into (1), the following
dispersion relation of 2D KdV equation is obtained

$\Omega_{1}=-\frac{p_{1}^{2}(p_{1}q_{1}-3)}{q_{1}}.$

The exact expression of $U_{[1]}$ and $V_{[1]}$ are as follows

$\displaystyle U_{[1]}$
$\displaystyle=\phi(y)\frac{1+\cosh\left(p_{1}x+q_{1}\Phi(y)+(\frac{3p^{2}_{1}}{q_{1}}-p^{3}_{1})t\right)-p_{1}q_{1}}{1+\cosh\left(p_{1}x+q_{1}\Phi(y)+(\frac{3p^{2}_{1}}{q_{1}}-p^{3}_{1})t\right)},$
(29) $\displaystyle V_{[1]}$
$\displaystyle=-\frac{p^{2}_{1}}{1+\cosh\left(p_{1}x+q_{1}\Phi(y)+(\frac{3p^{2}_{1}}{q_{1}}-p^{3}_{1})t\right)}.$

As can be seen from the above expression, the extremum lines
$p_{1}x+q_{1}\Phi(y)+(\frac{3p^{2}_{1}}{q_{1}}-p^{3}_{1})t$ and velocity
$\frac{3p_{1}-q_{1}p_{1}^{2}}{q_{1}}$ of $U_{[1]}$ and $V_{[1]}$ are the same,
and the corresponding amplitudes of $U_{[1]}$ and $V_{[1]}$ in $(x,y)$-plane
are $\phi(y)(1-\frac{p_{1}q_{1}}{2})$ and $\frac{p_{1}^{2}}{2}$ respectively.
Thus, $U_{[1]}$ is a deformed bright soliton when $\phi(y)(2-p_{1}q_{1})>0$
and a deformed dark soliton when $\phi(y)(2-p_{1}q_{1})<0$. However, $V_{[1]}$
is always a deformed dark soliton see Fig.1.

(a) $\Phi(y)=\operatorname{sech}(y)$

(b) $\Phi(y)=\operatorname{sech}(y)$

(c) $\Phi(y)=\sin(y)$

(d) $\Phi(y)=\sin(y)$

Figure 1: Deformed one-soliton $U_{[1]}$ and $V_{[1]}$ of equation(1) in the
$(x,y)$-plane with parameters $p_{1}=1,q_{1}=-4$ and displayed at $t=0$.

In order to obtain deformed two-solitons of 2D KdV equation(1), taking

$\displaystyle
U_{[2]}=-2\partial_{xy}\ln(1+e^{\eta_{1}}+e^{\eta_{2}}+e^{\eta_{1}+\eta_{2}+A_{12}})+\phi(y),$
(30) $\displaystyle
V_{[2]}=-2\partial_{xx}\ln(1+e^{\eta_{1}}+e^{\eta_{2}}+e^{\eta_{1}+\eta_{2}+A_{12}}),$

where $\eta_{1}=p_{1}x+\Phi(y)q_{1}y+\Omega_{1}t+\eta^{0}_{1}$ and
$\eta_{2}=p_{2}x+\Phi(y)q_{2}y+\Omega_{2}t+\eta^{0}_{2}$. Substituting (30)
into (1), the following relations are obtained

$\displaystyle\Omega_{1}$
$\displaystyle=-\frac{p_{1}^{2}(p_{1}q_{1}-3)}{q_{1}},\quad\Omega_{2}=-\frac{p_{2}^{2}(p_{2}q_{2}-3)}{q_{2}},$
(31) $\displaystyle\exp(A_{12})$
$\displaystyle=\frac{(p_{1}q_{2}-p_{2}q_{1})^{2}+p_{1}p_{2}q_{1}q_{2}(p_{1}-p_{2})(q_{1}-q_{2})}{(p_{1}q_{2}-p_{2}q_{1})^{2}+p_{1}p_{2}q_{1}q_{2}(p_{1}+p_{2})(q_{1}+q_{2})}.$

Further, taking
$p_{1}=3,q_{1}=3,p_{2}=\frac{5}{2},q_{2}=3,\eta^{0}_{1}=-\eta^{0}_{2}=\frac{\pi}{2}$,
the analytical expressions $U_{[2]}$ and $V_{[2]}$ of deformed two-soliton
solutions are as follows

$\displaystyle U_{[2]}$
$\displaystyle=\frac{\begin{Bmatrix}&\cosh(2\kappa_{11})-\sinh(2\kappa_{11})+982081\left(\sum\limits_{j=2}^{4}[\cosh(2\kappa_{1j})+\sinh(2\kappa_{1j})]\right)\\\
&+\left(\sum\limits_{m=1}^{5}\Gamma_{m}[\cosh(\kappa_{m})-(-1)^{m}\sinh(\kappa_{m})]\right)\end{Bmatrix}}{\frac{1}{\phi(y)}\left(\cosh(\kappa_{11})-\sinh(\kappa_{11})+991\sum\limits_{j=2}^{4}[\cosh(\kappa_{1j})+\sinh(\kappa_{1j})]\right)^{2}},$
(32) $\displaystyle V_{[2]}$
$\displaystyle=\frac{\begin{Bmatrix}&\left[\cosh(3\Phi(y))+\sinh(3\Phi(y))\right](\sum\limits_{m=1}^{5}\gamma_{m}[\cosh(\kappa_{m}-3(-1)^{m+1}\Phi(y))]\\\
&-\sum\limits_{m=1}^{5}(-1)^{m}[\sinh(\kappa_{m}-3(-1)^{m+1}\Phi(y))])\end{Bmatrix}}{-1982\left(\frac{1}{991}[\cosh(\kappa_{11})-\sinh(\kappa_{11})]+\sum\limits_{j=2}^{4}[\cosh(\kappa_{1j})+\sinh(\kappa_{1j})]\right)^{2}},$

where

$\displaystyle\kappa_{11}$
$\displaystyle=-\frac{11}{2}x-6\Phi(y)+\frac{75}{8}t-\frac{\pi}{2},\quad\kappa_{12}=\frac{5}{2}x+3\Phi(y)+\frac{69}{8}t,\quad\kappa_{13}=18t+\frac{\pi}{2},\quad\kappa_{14}=3x+3\Phi(y),$
(33) $\displaystyle\kappa_{1}$
$\displaystyle=\kappa_{12}+\kappa_{14},\quad\kappa_{2}=\kappa_{11}-\kappa_{12},\quad\kappa_{3}=\kappa_{12}+\kappa_{13},\quad\kappa_{4}=\kappa_{11}-\kappa_{14},\quad\kappa_{5}=\kappa_{13}+\kappa_{14},\quad\Gamma_{1}=1900738,$
$\displaystyle\Gamma_{2}$
$\displaystyle=-15856,\quad\Gamma_{3}=-12767053,\quad\Gamma_{4}=-12883,\quad\Gamma_{5}=15713296,\quad\gamma_{1}=1112,\quad\gamma_{2}=36,$
$\displaystyle\gamma_{3}$
$\displaystyle=24775,\quad\gamma_{4}=25,\quad\gamma_{5}=35676.$

The dynamic behaviors of deformed two-soliton solutions are more complex and
interesting by choosing the appropriate parameter $\phi(y)$ see Fig.2.
Furthermore, as can be seen the two-dimensional plots of Fig.2, the
interaction of deformed two-soliton solutions is an elastic collision.

(a) $\Phi(y)=\operatorname{sech}(y)$

(b) $\Phi(y)=\operatorname{sech}(y)$

(c) $\Phi(y)=\sin(y)$

(d) $\Phi(y)=\sin(y)$

(e) $\Phi(y)=\operatorname{sech}(y)$

(f) $\Phi(y)=\operatorname{sech}(y)$

(g) $\Phi(y)=\sin(y)$

(h) $\Phi(y)=\sin(y)$

Figure 2: Deformed two-soliton solutions $U_{[2]}$ and $V_{[2]}$ of
equation(1) in the $(x,y)$-plane displayed at $t=0$. Panels (a),(b),(c),(d)
are the two-dimensional plots of (e),(f),(g),(h) respectively

Similarly, $N$-soliton solutions $U_{[N]}$ and $V_{[N]}$ are given in
equation(22) of the 2D KdV equation(1), in which $f$ can be written as
follows:

$\displaystyle
f=f_{[N]}=\sum_{\mu=0,1}\exp(\sum_{j<k}^{(N)}\mu_{j}\mu_{k}A_{jk}+\sum_{j=1}^{N}\mu_{j}\eta_{j}),$
(34)

here

$\displaystyle\eta_{j}=p_{j}x+\Phi(y)q_{j}y+\Omega_{j}t+\eta^{0}_{j},\quad\Omega_{j}=-\frac{p_{j}^{2}(p_{j}q_{j}-3)}{q_{j}},$
(35)
$\displaystyle\exp(A_{jk})=\frac{(p_{j}q_{k}-p_{k}q_{j})^{2}+p_{j}p_{k}q_{j}q_{k}(p_{j}-p_{k})(q_{j}-q_{k})}{(p_{j}q_{k}-p_{k}q_{j})^{2}+p_{j}p_{k}q_{j}q_{k}(p_{j}+p_{k})(q_{j}+q_{k})},$

where $p_{j},q_{j}$ are arbitrary real parameters, $\eta^{0}_{j}$ is a complex
constant, $\phi(y)$ is an arbitrary function of $y$, and the subscript $j$
denotes an integer. For higher order soliton solutions, more deformed solitons
will be generated from equation(34) with appropriate parameter $\phi(y)$, and
its dynamic behavior will be more complex, we did not discuss it here.

##  4 Deformed 2D RW of the 2D KdV equation

In order to obtain the deformed 2D RW solutions of the 2D KdV equation(1),
taking

$\displaystyle q_{1}=\lambda_{1}p_{1},\quad
q_{2}=\lambda_{2}p_{2},\quad\eta_{1}^{0}=\eta_{2}^{0}=i\pi,$ (36)

in equation(30) and take a suitable limit as $p_{1},p_{2}\rightarrow 0$, then
we further take $\lambda_{1}=a+bi,\lambda_{2}=a-bi$ in (36). The solutions
$U^{[1]}$ and $V^{[1]}$ are given as follows

$\displaystyle U^{[1]}(x,y,t)$
$\displaystyle=\phi(y)+4b^{2}\phi(y)\frac{ab^{2}(a^{2}+b^{2})A_{21}^{2}-\frac{b^{4}}{a}A_{22}^{2}-A_{23}}{(A_{11}^{2}+A_{12}^{2}+A_{13})^{2}},$
(37) $\displaystyle V^{[1]}(x,y,t)$
$\displaystyle=4b^{2}\frac{b^{2}(a^{2}-b^{2})(a^{2}+b^{2})^{2}B_{11}^{2}-\frac{b^{4}}{a^{2}-b^{2}}B_{12}^{2}-B_{13}}{(A_{11}^{2}+A_{12}^{2}+A_{13})^{2}},$

with

$\displaystyle A_{11}$
$\displaystyle=\frac{3b(a^{2}-b^{2})}{a^{2}+b^{2}}t+abx+b(a^{2}+b^{2})\Phi(y),\quad
A_{12}=b^{2}x+\frac{6ab^{2}}{a^{2}+b^{2}}t,$ (38) $\displaystyle A_{13}$
$\displaystyle=a(a^{2}+b^{2})^{2},\quad
A_{21}=\frac{a^{2}+b^{2}}{a}x+(a^{2}+b^{2})\Phi(y)+3t,$ $\displaystyle A_{22}$
$\displaystyle=6at+(a^{2}+b^{2})x,A_{23}=a^{2}(a^{2}+b^{2})^{3},\quad
B_{11}=\frac{a}{a^{2}-b^{2}}x+\Phi(y)+\frac{3}{a^{2}-b^{2}}t,$ $\displaystyle
B_{12}$ $\displaystyle=6at+(a^{2}+b^{2})x,\quad B_{13}=a(a^{2}+b^{2})^{3}.$

As can be seen from the above expression, for ensure that the above solutions
$U^{[1]}$ and $V^{[1]}$ are smooth, parameter $a>0$ must be held.

###  4.1 Fundamental rational solutions

The fundamental rational solutions $U^{[1]}$ and $V^{[1]}$ of 2D KdV
equation(1) are derived if $\phi(y)$ is a polynomial function. Without loss of
generality, take $\phi(y)=2y$. The trajectories of $U^{[1]}$ and $V^{[1]}$ are
as follows

$\displaystyle l_{1}$
$\displaystyle=ax+(a^{2}+b^{2})y^{2}+\frac{3(a^{2}-b^{2})}{a^{2}+b^{2}}t,$
(39) $\displaystyle l_{2}$ $\displaystyle=x+\frac{6a}{a^{2}+b^{2}}t.$

$\lim\limits_{(x,y)\rightarrow\infty}U^{[1]}=2y$ and
$\lim\limits_{(x,y)\rightarrow\infty}V^{[1]}=0$, which show that the
background planes of $U^{[1]}$ and $V^{[1]}$ are $2y$ and 0, respectively.
Further, taking $a=2$ and $b=2$, the exact expression of fundamental rational
solutions $U^{[1]}$ and $V^{[1]}$ can be obtained. In order to describe the
evolution process of fundamental rational solutions $U^{[1]}$ more clearly, we
removed the background plane $2y$ when plotting. As shown in Fig.3. the
rational solution $U^{[1]}$ removing the background plane $2y$ appears from a
constant plane when $t\ll 0$. With time evolution an arc line $x=-4y^{2}$
appears on the constant background around at $t=0$ [see the Fig.3(a)], and
finally the arc line fission into a bright lump and a dark lump [see the
Fig.3(c)]. The trajectories of fission are $x+4y^{2}=0$ and $2x+3t=0$ [see the
Fig.3(d)(e)].

(a) t=0

(b) t=5

(c) t=10

(d) t=10

Figure 3: The temporal evolution of fundamental rational solution $U_{[1]}$
removing the background plane $\phi(y)$ of the equation(1) in the
$(x,y)$-plane with parameters $\phi(y)=2y,a=2$ and $b=2$; Panels (d) and (e)
are two-dimensional plots of $U_{[1]}$.

However, the dynamic behavior of the rational solution $V^{[1]}$ is different
from $U^{[1]}$. For given $\phi(y)=2y$, the solution $V^{[1]}$ has the
following nine critical points in $(x,y)$-plane

$\displaystyle\Lambda_{1}=(x_{1},y_{1})=(\frac{3at}{a^{2}+b^{2}},y=0);$ (40)
$\displaystyle\Lambda_{2,3}=(x_{2,3},y_{2,3})=(\frac{-3abt\pm\sqrt{3a(a^{2}+b^{2})^{3}+27b^{4}t^{2}}}{b(a^{2}+b^{2})},y=0);$
$\displaystyle\Lambda_{4,5}=(x_{4,5},y_{4,5})=(-\frac{6at}{a^{2}+b^{2}},y=\pm\sqrt{\frac{3t}{a^{2}+b^{2}}});$
$\displaystyle\Lambda_{6,7,8,9}=(x_{6,7,8,9},y_{6,7,8,9})=(\frac{-6abt\pm\sqrt{3a(a^{2}+b^{2})^{3}}}{b(a^{2}+b^{2})},y=\pm\sqrt{\frac{3t}{a^{2}+b^{2}}}),$

letting

$\displaystyle\triangle(x,y)=\frac{\partial^{2}V^{[1]}}{\partial
x^{2}};\quad\mathbf{H}(x,y)=\frac{\partial^{2}V^{[1]}}{\partial
x^{2}}\frac{\partial^{2}V^{[1]}}{\partial
y^{2}}-(\frac{\partial^{2}V^{[1]}}{\partial x\partial y})^{2},$ (41)

yield

$\displaystyle\triangle(x,y)\Bigg{|}_{\Lambda_{1}}=\frac{24b^{4}(a^{2}+b^{2})^{4}}{\left[a(a^{2}+b^{2})^{3}+9b^{4}t^{2}\right]^{2}};\qquad\mathbf{H}(x,y)\Bigg{|}_{\Lambda_{1}}=-\frac{1152b^{10}(a^{2}+b^{2})^{7}t}{\left[a(a^{2}+b^{2})^{3}+9b^{4}t^{2}\right]^{4}},$
(42)
$\displaystyle\triangle(x,y)\Bigg{|}_{\Lambda_{2,3}}=-\frac{3b^{4}(a^{2}+b^{2})^{4}}{4[a(a^{2}+b^{2})^{3}+9b^{4}t^{2}]^{2}};\quad\mathbf{H}(x,y)\Bigg{|}_{\Lambda_{2,3}}=-\frac{9b^{10}(a^{2}+b^{2})^{7}t}{2[a(a^{2}+b^{2})^{3}+9b^{4}t^{2}]^{4}},$
$\displaystyle\triangle(x,y)\Bigg{|}_{\Lambda_{4,5}}=\frac{24b^{4}}{a^{2}(a^{2}+b^{2})^{2}};\qquad\qquad\qquad\mathbf{H}(x,y)\Bigg{|}_{\Lambda_{4,5}}=\frac{2304b^{1}0t}{a^{4}(a^{2}+b^{2})^{5}},$
$\displaystyle\triangle(x,y)\Bigg{|}_{\Lambda_{6,7,8,9}}=-\frac{3b^{4}}{4a^{2}(a^{2}+b^{2})^{2}};\qquad\qquad\mathbf{H}(x,y)\Bigg{|}_{\Lambda_{6,7,8,9}}=\frac{9b^{10}t}{a^{4}(a^{2}+b^{2})^{5}}.$

Based on the above analysis, the evolution of rational solution $V^{[1]}$ of
2D KdV equation(1) can be divided into the following three stages.

  * •

(i) When $t<0$, $\lim\limits_{(x,y)\rightarrow\infty}V^{[1]}=0$, which show
that rational solution $V^{[1]}$ appears from a constant background. Now the
maximum value $V^{[1]}_{max}(x,y)$ and minimum value $V^{[1]}_{min}(x,y)$ of
$V^{[1]}$ are obtained at $\Lambda_{2,3}$ and $\Lambda_{1}$

$\displaystyle
V^{[1]}_{max}(x,y)=V^{[1]}(x,y)\Bigg{|}_{\Lambda_{2,3}}=\frac{b^{2}(a^{2}+b^{2})^{2}}{2a(a^{2}+b^{2})^{3}+18b^{4}t^{2}},$
(43) $\displaystyle
V^{[1]}_{min}(x,y)=V^{[1]}(x,y)\Bigg{|}_{\Lambda_{1}}=-\frac{4b^{2}(a^{2}+b^{2})^{2}}{a(a^{2}+b^{2})^{3}+9b^{4}t^{2}}.$

Obviously, the global extreme values of rational solution $V^{[1]}$ changes
with time, and

$\lim\limits_{t\rightarrow-\infty}V^{[1]}_{max}(x,y)=\lim\limits_{t\rightarrow-\infty}V^{[1]}_{min}(x,y).$

  * •

(ii) When $t=0$, rational solution $V^{[1]}$ has three extreme points

$\displaystyle A_{1}(x,y)=(0,0),\quad
A_{2}(x,y)=(\frac{\sqrt{3a^{2}+3ab^{2}}}{b},0),\quad
A_{3}(x,y)=(-\frac{\sqrt{3a^{2}+3ab^{2}}}{b},0).$

Maximum and minimum values are as follows

$\displaystyle
V^{[1]}_{max}(x,y)=V^{[1]}(x,y)\Bigg{|}_{A_{2}}=V^{[1]}(x,y)\Bigg{|}_{A_{3}}=\frac{b^{2}}{2a(a^{2}+b^{2})},$
(44) $\displaystyle
V^{[1]}_{min}(x,y)=V^{[1]}(x,y)\Bigg{|}_{A_{1}}=-\frac{4b^{2}}{a(a^{2}+b^{2})}.$

From the above analysis, we can seen that its dynamic behavior is similar to
the RW in one-dimensional systems. The amplitude of the RW in one-dimensional
systems is three times that of the background plane. However, the amplitude of
RW of 2D KdV equation(1) is controlled by parameters $a$ and $b$.

  * •

(iii) When $t>0$, rational solution $V^{[1]}$ has six extreme points
$\Lambda_{4,5}$ and $\Lambda_{6,7,8,9}$. And it has four maximums and two
minimums

$\displaystyle
V^{[1]}_{max}(x,y)=V^{[1]}(x,y)\Bigg{|}_{\Lambda_{6,7,8,9}}=\frac{b^{2}}{2a(a^{2}+b^{2})},$
(45) $\displaystyle
V^{[1]}_{min}(x,y)=V^{[1]}(x,y)\Bigg{|}_{\Lambda_{4,5}}=-\frac{4b^{2}}{a(a^{2}+b^{2})}.$

Obviously, the extreme values of $V^{[1]}$ at $t=0$ is equal to the extreme
values at $t>0$.

Through the above analysis, $|V^{[1]}_{max}(x,y)|<|V^{[1]}_{min}(x,y)|$ always
satisfied, which indicated that $V^{[1]}$ is dark. When $t<0$, we know that
the rational solution $V^{[1]}$ appears from a constant plane, and the
amplitude above the background plane is $8$ times of that below the background
plane. $V^{[1]}$ reaches the maximum amplitude
$\frac{9b^{2}}{2a(a^{2}+b^{2})}$ at $t=0$. With the time involution, dark lump
fission into two identical dark lumps along the curve
$ax+(a^{2}+b^{2})y^{2}+\frac{3(a^{2}-b^{2})}{a^{2}+b^{2}}t$. Interestingly,
after the dark lump fission into two identical lumps, its amplitude does not
change, which means the energy does not change when $t>0$ [see the Fig.4].

(a) t=-50

(b) t=-9

(c) t=0

(d) t=12

Figure 4: The temporal evolution of fundamental rational solution $V^{[1]}$ of
the equation(1) in the $(x,y)$-plane with parameters $\phi(y)=2y,a=2$ and
$b=2$.

###  4.2 Fundamental deformed 2D RW solutions

The fundamental deformed 2D RW $U^{[1]}$ and $V^{[1]}$ of 2D KdV equation(1)
are obtained when $\phi(y)=\operatorname{sech}(y)$.
$\lim\limits_{(x,y)\rightarrow\infty}U^{[1]}=\phi(y)$, which indicated that
the background plane of deformed 2D RW $U^{[1]}$ is $\phi(y)$. The
trajectories of $U^{[1]}$ are

$\displaystyle l^{{}^{\prime}}_{1}$
$\displaystyle=ax+(a^{2}-b^{2})\operatorname{sech}(y)+3t,$ (46) $\displaystyle
l^{{}^{\prime}}_{2}$ $\displaystyle=x+2a\operatorname{sech}(y).$

As can be seen in Fig.5, in order to better observe the evolution of the
deformed 2D RW $U^{[1]}$, we remove the background plane of $U^{[1]}$. Four
panels describe the appearance and annihilation of 2D RW in $(x,y)$-plane
along the curve $x+2\operatorname{sech}(y)$. This is the first time to obtain
such deformed 2D RW in high-dimensional systems. The dynamic behavior of
fundamental deformed 2D RW $V^{[1]}$ is more complicated and interesting.

(a) t=-30

(b) t=-2

(c) t=3

(d) t=30

(e) t=-2

(f) t=3

Figure 5: The temporal evolution of fundamental deformed 2D RW solution
$U^{[1]}$ removing the background plane $\phi(y)$ of the equation(1) in the
$(x,y)$-plane with parameters $\Phi(y)=\operatorname{sech}(y),a=1$ and $b=4$.

Through simple calculation and analysis, the evolution of deformed 2D RW
$V^{[1]}$ can be divided into the following four stages

  * •

(i) When $t<0$, a line RW appears from constant background plane, with the
time involution, the amplitude of the line RW gradually increases.

  * •

(ii) When $t=0$, the amplitude of the line RW reaches its maximum value. The
maximum amplitude $V^{[1]}_{Amp}$ is

$\displaystyle
V^{[1]}_{Amp}=V^{[1]}_{Max}-V^{[1]}_{Min}=\frac{b^{2}}{2a(a^{2}+b^{2})}-\frac{-4b^{2}}{a(a^{2}+b^{2})}=\frac{9b^{2}}{2a(a^{2}+b^{2})}.$
(47)

  * •

(iii) When $0<t\leq\frac{a^{2}+b^{2}}{3}$, the line RW disappears gradually.
Simultaneously, a Peregrine-type soliton is generated, but the amplitude of
the RW $V^{[1]}$ is fixed in this process of evolution

  * •

(iv) When $t>\frac{a^{2}+b^{2}}{3}$, two identical maximum values
$V^{[1]}_{Max}$ and one minimum value $V^{[1]}_{Min}$ of the Peregrine-type
soliton are as follows

$\displaystyle
V^{[1]}_{Max}=\frac{b^{2}(a^{2}+b^{2})^{2}}{2b^{4}(3t-a^{2}-b^{2})^{2}+2a(a^{2}+b^{2})^{3}},$
(48) $\displaystyle
V^{[1]}_{Min}=-\frac{4b^{2}(a^{2}+b^{2})^{2}}{b^{4}(3t-a^{2}-b^{2})^{2}+a(a^{2}+b^{2})^{3}}.$

It can be seen from the above expression that the amplitude of the Peregrine-
type soliton decays rapidly in a very short time. When
$t\gg\frac{a^{2}+b^{2}}{3}$, RW $V^{[1]}$ uniformly approaches a constant
background see Fig.6.

(a) t=-50

(b) t=-6

(c) t=-2

(d) t=0

(e) t=2

(f) t=7

(g) t=50

Figure 6: The temporal evolution of fundamental deformed 2D RW solution
$V^{[1]}$ of the equation(1) in the $(x,y)$-plane with parameters
$\Phi(y)=\operatorname{sech}(y),a=1$ and $b=4$.

###  4.3 High-order deformed 2D RW solutions

When taking

$\displaystyle N=4,\quad
q_{j}=\lambda_{j}p_{j},\quad\eta_{j}^{0}=i\pi(j=1,2,3,4),$ (49)

in equation(34). Further taking

$\displaystyle\lambda_{1}=\lambda_{2}^{*}=1+3i,\quad\lambda_{3}=\lambda_{4}^{*}=\frac{1}{2}-4i,\quad\Phi(y)=\operatorname{sech}(y),$
(50)

the second-order deformed 2D RW solutions $U^{[2]}=-2\ln f^{[2]}_{xy}+\phi(y)$
and $V^{[2]}=-2\ln f^{[2]}_{xx}$ are generated, in which $f^{[2]}$ can be
written as,

$\displaystyle f^{[2]}(x,y,t)$
$\displaystyle=x^{4}+[6\Phi(y)+\frac{27}{10}t]x^{3}+[21\Phi(y)^{2}+\frac{21}{5}\Phi(y)t+\frac{189}{40}t^{2}+\frac{149877}{1156}]x^{2}+[36\Phi(y)^{3}+\frac{27}{10}t\Phi(y)^{2}$
(51)
$\displaystyle+\frac{81}{20}t^{2}\Phi(y)+\frac{81}{20}t^{3}+\frac{2924139}{11560}t]x+\frac{81}{40}t^{4}-\frac{81}{20}\Phi(y)t^{3}+[\frac{801}{40}\Phi(y)^{2}+\frac{688257}{46240}]t^{2}+[-\frac{144}{5}\Phi(y)^{3}$
$\displaystyle+\frac{1008648}{1445}\Phi(y)]t+40\Phi(y)^{4}-\frac{122850}{289}\Phi(y)^{2}+\frac{1158885}{89}.$

As can be seen in Fig.7, remove the background plane $\phi(y)$, the second-
order deformed 2D RW solution $U^{[2]}$ describes that four Peregrine-type
solitons appear and annihilate from the constant background plane. However,
the dynamic behaviors of the second-order RW solution $V^{[2]}$ are similar to
that of the fundamental RW solution $V^{[1]}$. The RW solution $V^{[2]}$ is
uniformly approach to a constant background plane when
$t\rightarrow\pm\infty$. Two line RWs appear from the constant plane under the
time evolution, and their amplitudes increased rapidly. Then the amplitude of
the two line RW attenuated rapidly, at the same time, two Peregrine-type
solitons are produced. Finally, these two Peregrine-type solitons disappeared
in a very short time without a trace see Fig.8.

(a) t=-20

(b) t=-2

(c) t=4

(d) t=20

Figure 7: The temporal evolution of second-order deformed 2D RW solution
$U^{[2]}$ removing the background plane $\phi(y)$ of equation(1) in the
$(x,y)$-plane.

(a) t=-20

(b) t=-2

(c) t=0

(d) t=4

(e) t=20

Figure 8: The temporal evolution of second-order deformed 2D RW solution
$V^{[2]}$ of equation(1) in the $(x,y)$-plane.

For larger N, higher-order solutions are generated, taking the parameters

$\displaystyle N=2n,\quad
q_{j}=\lambda_{j}p_{j},\quad\eta_{j}^{0}=i\pi\quad(1\leq j\leq N),$ (52)

in equation(34) and take a suitable long wave limit as $p_{j}\rightarrow 0$,
the functions $f$ defined in (34) become a polynomial-type function containing
an arbitrary function $\phi(y)$. Therefore, the general $n$-th rational-type
functions $U^{[n]}=-2\ln f_{xy}+\phi(y)$ and $V^{[n]}=-2\ln f_{xx}$ of 2D KdV
equation(1) can be derived[55], in which $f$ can be written as follows,

$\displaystyle
f=f^{[n]}=\prod_{k=1}^{N}\theta_{k}+\frac{1}{2}\sum_{k,j}^{(N)}\alpha_{kj}\prod_{l\neq
k,j}^{N}\theta_{l}+\cdots+\frac{1}{M!2^{M}}\sum_{i,j,...,m,n}^{(N)}\overbrace{\alpha_{kj}\alpha_{kl}\cdots\alpha_{mn}}^{M}\prod_{p\neq
k,j,...m,n}^{N}\theta_{p}+\cdots,$ (53)

with

$\displaystyle\theta_{j}=\frac{\lambda_{j}^{2}\phi(y)+\lambda_{j}x+3t}{\lambda_{j}},\quad\alpha_{jk}=-\frac{2\lambda_{j}\lambda_{k}(\lambda_{j}+\lambda_{k})}{(\lambda_{j}-\lambda_{k})^{2}},$
(54)

the two positive integers $k$ and $j$. We must emphasize that $\lambda_{j}$ is
a complex constant and $\lambda^{*}_{j}=\lambda_{n+j}=a_{n}+ib_{n}$. When
$a_{n}>0$, the above general $n$-th rational-type functions $U^{[n]}$ and
$V^{[n]}$ are smooth.

Remark 1 The above solutions $U^{[n]}$ and $V^{[n]}$ are rational solutions,
if $\phi(y)$ is a nonzero polynomial function. For example, when $\phi(y)=2y$,
the rational solution $U^{[n]}$ describes the fission of $n$-bright lumps and
$n$-dark lumps from the background plane $2y$. The rational solution $V^{[n]}$
describes the fission of $2n$-dark lumps from a constant background plane.

Remark 2 The above solutions $U^{[n]}$ and $V^{[n]}$ are deformed 2D RW
solutions, if $\Phi(y)=\operatorname{sech}(y)$. The solution $U^{[n]}$
describes $2n$ Peregrine-type solitons appear and annihilate from a kink-
soliton plan. The solution $V^{[n]}$ shows that $n$-line RWs appear and decay
rapidly from a constant background plane, and fission into $n$ Peregrine-type
solitons, and finally disappear in the constant background plane without a
trace.

##  5 Summary

In this letter, new bilinear B$\rm\ddot{a}$cklund transformation and lax pair
of the 2D KdV equation(1) are derived, which are different from the
B$\rm\ddot{a}$cklund transform and lax pair in Ref.[36, 37]. N-soliton
solutions are presented by means of the improved Hirota’s bilinear method.
Deformed soliton and deformed breather solutions of elastic collision are
generated by selecting the appropriate free parameter $\phi(y)$ see Fig.1 and
Fig.2. When $\Phi(y)$ is a non-constant polynomial function, a family of new
rational solutions of the 2D KdV equation are generated using the long wave
limit. When $\Phi(y)$ ia a non-zero constant, the rational solutions $U^{[n]}$
and $V^{[n]}$ reduce to the rational solution of 2D KdV equation in Ref.[47].
When $\Phi(y)=\operatorname{sech}(y)$, the deformed 2D RW solution $U$
describes a family of Peregrine-type solitons appear and annihilate from a
kink-soliton plan. In order to better observe the evolution of the deformed 2D
RW solutions, we remove the background plane of $U$ when plotting see Figs.5
and 7. The deformed 2D RW solution $V$ shows that a series of line RWs appear
and decay rapidly from a constant background plane, and fission into
Peregrine-type solitons, and finally annihilate in the constant background
plane without a trace see Figs.6 and 8. This paper successfully constructed
the deformed 2D RW solutions, which are closely related to the introduced
arbitrary function $\Phi(y)$. These novel phenomena have never been reported
before in nonlinear systems. Our presented work not only provides a new
reference method for seeking new exact solutions of nonlinear partial
differential equations, but also may be helpful to promote a deeper
understanding of nonlinear phenomena.

## Acknowledgments

This work is supported by the NSF of China under Grants No.11671219,
No.11871446, No.12071304 and No.12071451.

## References

  * [1] Akhmediev N, Ankiewicz A and Taki M 2009 Phys. Lett. A 373 675-678
  * [2] Peregrine D H 1983 J. Aust. Math. Soc. B 25 16-43
  * [3] Guo L J, Wang L H, Cheng Y and He J S 2017 Commun. Nonlinear Sci. Numer. Simulat. 52 11-23
  * [4] Zhang Y S, Guo L J, Chabchoub A and He J S 2017 Rom. J. Phys. 62 102
  * [5] Zhang G Q and Yan Z Y 2018 Commun. Nonlinear Sci. Numer. Simulat. 62 117-133
  * [6] He J S, Xu S W and Porsezian K 2012 Phys. Rev. E 86 066603
  * [7] Li R M, Geng X G and Xue B 2020 J. Nonlinear Math. Phys. 27 279-294
  * [8] Yang J and Zhu Z N 2018 chaos 28 093103
  * [9] Ling L M, Feng B F and Zhu Z N 2016 Physica D 327 13-29
  * [10] Cao Y L, He J S, Cheng Y and Mihalache D 2020 Nonlinear Dyn. 99 3013-3028
  * [11] Gai L T, Ma W X and Li M C 2020 Phys. Lett. A 384 126178
  * [12] Rao J G, Zhang Y S, Fokas A S and He J S 2018 Nonlinerity 31 4090-4107
  * [13] Liu W H, Zhang Y F and Shi D D 2019 Phys. Lett. A 383 97-102
  * [14] Rao J G, Wang L H, Liu W and He J S 2017 Theor. Math. Phys. 193 1783-1800
  * [15] Cao Y L, Malomed B A and He J S 2018 Chaos Soliton and Fract. 114 99-107
  * [16] Ohta Y and Yang J K 2013 J. Phys. A: Math. Theor. 46 105202
  * [17] Rao J G, Mihalache D, Cheng Y and He J S 2019 Phys. Lett. A 383 1138-1142
  * [18] Zhang Y S, Rao J G, Porsezian K and He J S 2019 Nonlinear Dyn. 95 1133-1146
  * [19] Elawady E and Moslem W M 2011 Phys. Plasmas 18 082306
  * [20] Bailung H, Sharma S K and Nakamura Y 2011 Phys. Rev. Lett. 107 255005
  * [21] Bludov Y V, Konotop V V and Akhmediev N 2009 Phys. Rev. A 80 2962-2964
  * [22] Bludov Y V, Konotop V V and Akhmediev N 2010 Eur. Phys. J. Spec. Top. 185 169-180
  * [23] Stenflo L and Marklund M 2010 J. Plasma Phys. 76 293-295
  * [24] Montina A, Bortolozzo U, Residori S and Arecchi F T 2009 Phys. Rev. Lett. 103 173901
  * [25] Solli D R, Ropers C, Koonath P and Jalali B 2007 Nature(London) 450 1054-1057
  * [26] Mihalache D 2017 Rom. Rep. Phys. 69 403
  * [27] Ganshin A N, Efimov V B, Kolmakov G V, Mezhov-Deglin L P and Mcclintock P V E 2008 Phys. Rev. Lett. 101 065303
  * [28] Guo L J, He J S, Wang L H, Cheng Y, Frantzeskakis D J, Bremer T S and Kevrekidis P G 2020 Physical Review Research 2 033376
  * [29] Boiti M, Leon J J P, Manna M and Pempinelli F 1986 Inverse Probl. 2 271
  * [30] Lou S Y and Hu X B 1997 J. Math. Phys. 38 6401-6427
  * [31] Estevez P G and Leble S 1995 Inverse Probl. 11 925
  * [32] Leble S B and Ustinov N V 1991 Inverse Probl. 10 617
  * [33] Hirota R and Satsma J 1994 J. Phys. Soc. Jpn. 40 611
  * [34] Ablowitz M J and Clarkson P A 1991 Solitons, Nonlinear Evolution Equations and Inverse Scattering (Cambridge: Cambridge University Press)
  * [35] Delisle L and Mosaddeghi M 2013 J. Phys. A: Math. Theor. 46 115203
  * [36] Lue X, Tian B, Sun K and Wang P 2010 J. Math. Phys. 51 113506
  * [37] Luo L 2011 Phys. Lett. A 375 1059-1063
  * [38] Wang C J 2017 Nonlinear Dyn. 87 2635-2642
  * [39] Clarkson P A and Mansfield E L 1994 Nonlinearity 7 795
  * [40] Tian B and Gao Y T 1996 Chaos Soliton Fract. 7 1497-1499
  * [41] Lou S Y 1995 J. Phys. A: Math. Theor. 28 7227-7232
  * [42] Hu H C, Tang X Y, Lou S Y and Liu Q P 2004 Chaos Soliton Fract. 22 327
  * [43] Chen Y, Wang Q and Li B 2004 Commun. Theor. Phys. 42 655-660.
  * [44] Fan E G 2009 J. Phys. A 42 095206
  * [45] Luo L 2010 Commun. Theor. Phys. 54 208-214
  * [46] Chen Y R, Song M and Liu Z R 2015 Nonlinear Dyn. 82 333-347
  * [47] Liu Y Q and Wen X Y 2019 Adv. Diff. Equa. 2019 332
  * [48] Wang C J 2016 Nonlinear Dyn. 84 697-702
  * [49] Gilson C, Lambert F, Nimmo J and Willox R 1996 Proc. Roy. Soc. Lond. A 452 223
  * [50] Bell E T 1935 Ann. Math. 35 258-277
  * [51] Fan E G 2011 Phys. Lett. A 375 493-497
  * [52] Luo L 2019 Appl. Math. Lett. 94 94-98
  * [53] Lambert F and Springael J 2001 Chaos Soliton Fract. 12 2821
  * [54] Hirota R 2004 The direct method in soliton theory (Cambridge: Cambridge University Press)
  * [55] Ablowitz M J and Satsuma J 1979 J. Math. Phys. 20 1496-1503