# Pattern transformation in higher-order lumps of the Kadomtsev-Petviashvili I
equation

Bo Yang and Jianke Yang  Department of Mathematics and Statistics, University
of Vermont, Burlington, VT 05405, U.S.A.

###### Abstract

Pattern formation in higher-order lumps of the Kadomtsev-Petviashvili I
equation at large time is analytically studied. For a broad class of these
higher-order lumps, we show that two types of solution patterns appear at
large time. The first type of patterns comprise fundamental lumps arranged in
triangular shapes, which are described analytically by root structures of the
Yablonskii–Vorob’ev polynomials. As time evolves from large negative to large
positive, this triangular pattern reverses itself along the $x$-direction. The
second type of patterns comprise fundamental lumps arranged in non-triangular
shapes in the outer region, which are described analytically by nonzero-root
structures of the Wronskian–Hermit polynomials, together with possible
fundamental lumps arranged in triangular shapes in the inner region, which are
described analytically by root structures of the Yablonskii–Vorob’ev
polynomials. When time evolves from large negative to large positive, the non-
triangular pattern in the outer region switches its $x$ and $y$ directions,
while the triangular pattern in the inner region, if it arises, reverses its
direction along the $x$-axis. Our predicted patterns at large time are
compared to true solutions, and excellent agreement is observed.

##  I Introduction

The Kadomtsev-Petviashvili (KP) equation was derived as a two-dimensional
generalization of the Korteweg-de Vries equation for the evolution of weakly
nonlinear plasma waves and shallow water waves KP1970 ; Ablowitz1979 . In the
water wave context, this equation reads Ablowitz1979

$\left[2f_{t}+3ff_{x}+(\frac{1}{3}-T)f_{xxx}\right]_{x}+f_{yy}=0,$ (1)

where the spatial coordinate $x$ is relative to a certain moving frame,
$f(x,y,t)$ represents the water surface elevation, and $T$ is a dimensionless
surface tension parameter. If the surface tension is large, i.e., $T>1/3$,
which corresponds to very thin sheets of water, this equation is called KP-I.
In this case, rescaling variables by

$y=\frac{\hat{y}}{\sqrt{3(T-\frac{1}{3}})},\quad
t=-\frac{2\hat{t}}{T-\frac{1}{3}},\quad f=-2(T-\frac{1}{3})u$ (2)

and dropping the hats, this equation becomes

$\left(u_{t}+6uu_{x}+u_{xxx}\right)_{x}-3u_{yy}=0.$ (3)

Note that the KP-I equation also arises in other branches of physics, such as
nonlinear optics Pelinovsky1995  and Bose-Einstein condensates
Tsuchiya_BEC2008 .

The KP-I equation (3) is solvable by the inverse scattering transform
Zakharov_book ; Ablowitz_book . It admits stable fundamental lump solutions
that are bounded rational functions decaying in all spatial directions
Petviashvili1976 ; Manakov1977 ; Ablowitz_Satsuma1979 . These lumps are the
counterparts of solitons in the Korteweg-de Vries equation. In the water wave
context, these lumps physically correspond to dips on the water surface due to
the negative sign in the $f$ scaling above. The KP-I equation also admits a
broad class of rational solutions that describe the interactions of these
lumps. If individual lumps have distinct asymptotic velocities, then they
would pass through each other without change in velocities or phases
Manakov1977 ; Ablowitz_Satsuma1979 . But if they have the same asymptotic
velocities, they would undergo novel anomalous scattering, where the lumps
would separate from each other in new spatial directions that are very
different from their original incoming directions Peli93b ; Ablowitz97 ;
Ablowitz2000 . In this article, we are concerned with this latter type of
solutions, which we will call higher-order lumps (they are also called multi-
pole lumps in the literature Ablowitz97 ; Ablowitz2000 ).

Analytical expressions of higher-order lumps have been derived by a wide
variety of methods before Peli93a ; Peli93b ; Ablowitz97 ; Peli98 ;
Ablowitz2000 ; Dubard2010 ; Dubard2013 ; Ma2015 ; Chen2016 ; ClarksonDowie2017
; Gaillard2018 ; Chang2018 . Gorshkov, et al. Peli93b  reported a second-order
lump solution that describes the interaction and anomalous scattering of two
lumps. Ablowitz et al. Ablowitz97 ; Ablowitz2000  derived higher-order lumps
by the inverse scattering transform and Darboux transformation, and reproduced
the solution in Peli93b  as a special case. They also showed that when
$|t|\to\infty$, these higher-order lumps generically split into a certain
number of fundamental lumps, whose relative spatial separations grow in
proportion to $|t|^{q}$, where $\frac{1}{3}\leq q\leq\frac{1}{2}$. In
addition, some new lump patterns such as squares at large time were reported.
Pelinovsky and Stepanyants Peli93a  reported a class of higher-order lump
solutions that are stationary in a moving frame. Pelinovsky Peli98  studied
rational solutions of the KP hierarchy and linked them to the dynamics of the
Calogero–Moser hierarchy (but his Wronskian-form solutions for KP-I were not
made real-valued and thus were not physical solutions). Dubard et al.
Dubard2010 ; Dubard2013  constructed a class of higher-order KP-I lump
solutions from higher-order rogue waves of the nonlinear Schrödinger equation,
and graphically showed that such second- and third-order lump solutions split
into triangles of fundamental lumps when $|t|\to\infty$. Chen et al. Chen2016
considered a certain class of higher-order lump solutions, and graphically
observed that these solutions evolve from a vertical line of fundamental lumps
to a horizontal line of fundamental lumps in the $(x,y)$ plane when time goes
from negative infinity to infinity. They also predicted the locations of
fundamental lumps inside the solution complex at $t=0$ by roots of certain
polynomial equations; but such polynomial equations were not justified.
Clarkson and Dowie ClarksonDowie2017  derived a second-order lump solution
which incorporates the ones in Dubard2010 ; Dubard2013 ; Peli93b  as special
cases. Gaillard Gaillard2018  studied a special class of higher-order lump
solutions and reported lump patterns such as triangles and pentagons at $t=0$
when some internal parameters in such solutions get large. Chang Chang2018
studied the large-time asymptotics of higher-order lumps and showed that, for
some special solutions, all lumps are located on a vertical line in the
$(x,y)$ plane at large negative time but rotate to a horizontal line at large
positive time. Ma Ma2015  derived a fundamental lump solution which contains
more free parameters; but that solution can be made equivalent to the original
lump solution as reported in Manakov1977 ; Ablowitz_Satsuma1979 . We note by
passing that non-rational KP-I solutions in the form of a linear periodic
chain of lumps, and those that describe the resonant collision between lumps
and line solitons, have also been reported recently Zakharov2021 ;
He_KPI_rogue_lump .

In this article, we study pattern formation in higher-order lumps of the KP-I
equation (3). This work is motivated by our earlier work on pattern formation
of rogue waves in various integrable systems YangYang21a ; YangYang21b , where
we showed that universal rogue patterns appear when one of the internal
parameters in rogue waves gets large, and those rogue patterns are
analytically described by root structures of the Yablonskii–Vorob’ev
polynomial hierarchy. For higher-order lumps of the KP-I equation, however, we
will focus on their pattern formation at large time rather than at large
parameters. In particular, we are interested to know how their patterns at
large positive time relate to their patterns at large negative time. For a
broad class of higher-order lump solutions, we will show that two types of
lump patterns appear at large time. The first type of patterns comprise
fundamental lumps arranged in triangular shapes, which are described
analytically by root structures of the Yablonskii–Vorob’ev polynomials. As
time evolves from large negative to large positive, this triangular pattern
reverses itself along the $x$-direction. The second type of patterns comprise
fundamental lumps arranged in non-triangular shapes in the outer region, which
are described analytically by nonzero-root structures of the Wronskian–Hermit
polynomials, together with possible fundamental lumps arranged in triangular
shapes in the inner region, which are described analytically by root
structures of the Yablonskii–Vorob’ev polynomials. When time evolves from
large negative to large positive, the non-triangular pattern in the outer
region switches its $x$ and $y$ directions, plus some rescaling along each
direction, while the triangular pattern in the inner region, if it arises,
reverses its direction along the $x$-axis. These dramatic pattern
transformations with the elapse of time are fascinating. We have also compared
these predicted patterns with true solutions, and excellent agreement is
observed.

This paper is organized as follows. In Sec. 2, we present general higher-order
lump solutions in the KP-I equation through Schur polynomials, and introduce
Yablonskii–Vorob’ev and Wronskian–Hermit polynomials. In Sec. 3, we present
our main analytical results on solution patterns at large time, and explain
how these patterns transform from large negative time to large positive time.
In Sec. 4, we illustrate our pattern predictions and compare them with true
solutions. In Sec. 5, we provide proofs for our analytical results in Sec. 3.
The last section summarizes our results, together with some discussions. In
the Appendix, a brief derivation of our general higher-order lump solutions in
Sec. 2 is given.

##  II Preliminaries

The KP equation (3) admits three important invariances. The first one is that
it is invariant when $(x,t)\to(-x,-t)$. This invariance is important because
it shows that KP-solution patterns are reversible in time (albeit with a sign
switch in $x$). In earlier works Chen2016 ; Chang2018 , the authors showed
that certain higher-order KP lumps evolve from a vertical line of fundamental
lumps to a horizontal line of fundamental lumps in the $(x,y)$ plane when time
goes from negative infinity to infinity. The above invariance indicates that a
reverse pattern transformation could also occur, i.e., those higher-order KP
lumps can also evolve from a horizontal line to a vertical line when time goes
from negative infinity to infinity.

The second invariance of the KP equation (3) is the Galilean invariance
Weiss1985 ; Chen2016 , i.e., when

$(x,y,t)\to(x+2ky+12k^{2}t,\,y+12kt,\,t),$ (4)

the KP solution $u(x,y,t)$ remains a solution. Here, $k$ is an arbitrary real
constant. This invariance indicates that, if the overall solution complex has
a $y$-direction velocity $12k$, then we can apply this invariance to remove
that $y$-direction velocity. In doing so, the solution pattern in the $(x,y)$
plane would change as well through a linear transformation of shear type. This
Galilean invariance is important, because it allows us to remove the overall
$y$-direction velocity in a higher-order lump solution. More will be said on
it later in this section.

The third invariance of the KP equation is scaling invariance, i.e., when

$(x,y,t,u)\to(\alpha x,\alpha^{2}y,\alpha^{3}t,\alpha^{-2}u),$ (5)

the KP equation remains invariant. Here, $\alpha$ is any nonzero real
constant. This invariance is useful since, when combined with the Galilean
invariance above, it allows us to normalize the spectral parameter in the KP-
lump solutions to be unity without any loss of generality. This we will do in
Sec. 3.

###  II.1 Explicit expressions of higher-order lumps

In this paper, we consider pattern formation of higher-order lumps in the KP-I
equation (3). General higher-order lump solutions have been derived by
Ablowitz et al. Ablowitz2000  through Darboux transformation. Their solutions
were given through determinants whose matrix elements involve differential
operators with respect to the spectral parameter. For our analysis, those
solution expressions are not explicit enough. Thus, we have derived these
higher-order lumps again by the bilinear method. To present our solutions, we
first introduce elementary Schur polynomials $S_{k}(\mbox{\boldmath$x$})$ with
$\emph{{x}}=\left(x_{1},x_{2},\ldots\right)$, which are defined by the
generating function

$\sum_{n=0}^{\infty}S_{n}(\mbox{\boldmath$x$})\epsilon^{n}=\exp\left(\sum_{n=1}^{\infty}x_{n}\epsilon^{n}\right).$
(6)

More explicitly,

$S_{0}(\emph{{x}})=1,\quad S_{1}(\mbox{\boldmath$x$})=x_{1},\quad
S_{2}(\mbox{\boldmath$x$})=\frac{1}{2}x_{1}^{2}+x_{2},\quad\cdots,\quad
S_{n}(\mbox{\boldmath$x$})=\sum_{l_{1}+2l_{2}+\cdots+ml_{m}=n}\left(\
\prod_{j=1}^{m}\frac{x_{j}^{l_{j}}}{l_{j}!}\right).$

Under these notations, our general higher-order KP-I lumps are given by the
following theorem.

> Theorem 1 _General higher-order lumps of the KP-I equation ( 3) are _
>
> $\displaystyle u_{\Lambda}(x,y,t)=2\partial_{x}^{2}\ln\sigma,$ (7)
>
> _where_
>
> $\sigma(x,y,t)=\det_{1\leq i,j\leq N}\left(m_{ij}\right),$ (8)
> $m_{i,j}=\sum_{\nu=0}^{\min(n_{i},n_{j})}\left[\frac{|p|^{2}}{(p+p^{*})^{2}}\right]^{\nu}\hskip
> 1.70709ptS_{n_{i}-\nu}(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}}_{i})\hskip
> 1.70709ptS_{n_{j}-\nu}[(\textbf{\emph{x}}^{+})^{*}+\nu\textbf{\emph{s}}^{*}+\textbf{\emph{a}}_{j}^{*}],$
> (9)
>
> _$N$ is an arbitrary positive integer, $\Lambda\equiv(n_{1},n_{2},\cdots
> n_{N})$ is a vector of arbitrary positive integers,_ _$p$ is an arbitrary
> non-imaginary complex number, the asterisk ‘*’ represents complex
> conjugation, the vector
> $\textbf{\emph{x}}^{+}=\left(x_{1}^{+},x_{2}^{+},\cdots\right)$ is defined
> by _
>
> $\displaystyle
> x_{k}^{+}=p\frac{1}{k!}x+p^{2}\frac{2^{k}}{k!}\textrm{i}y+p^{3}\frac{3^{k}}{k!}(-4)t,$
> (10)
>
> _the vector $\textbf{\emph{s}}=(s_{1},s_{2},\cdots)$ is defined through the
> expansion_
>
>
> $\displaystyle\ln\left[\frac{1}{\kappa}\left(p+p^{*}\right)\left(\frac{e^{\kappa}-1}{p\hskip
> 1.42271pte^{\kappa}+p^{*}}\right)\right]=\sum_{j=1}^{\infty}s_{j}\hskip
> 1.42271pt\kappa^{j},$ (11)
>
> _vectors $\textbf{\emph{a}}_{i}$ are_
>
> $\textbf{\emph{a}}_{i}=\left(a_{i,1},a_{i,2},\cdots,a_{i,n_{i}}\right),$
> (12)
>
> _and $a_{i,j}\hskip 1.42271pt(1\leq i\leq N,1\leq j\leq n_{i})$ are free
> complex constants. _

The proof of this theorem will be given in the Appendix.

Remark 1. In this theorem, positive integers $(n_{1},n_{2},\cdots n_{N})$ do
not have to be distinct if their corresponding vectors $\textbf{\emph{a}}_{i}$
are different. In such cases, by first rewriting the $\sigma$ determinant (8)
as a larger determinant as was done in Ref. OhtaJY2012 , then linking Schur
polynomials with different $\textbf{\emph{a}}_{i}$ vectors in that larger
determinant by relations similar to Eq. (167) in Ref. YangYang3wave , and
finally applying row operations and parameter redefinitions to the resulting
determinant, we can show that this $\sigma$ determinant (8) with non-distinct
integers $(n_{1},n_{2},\cdots n_{N})$ can be reduced to one where the new
integers $(\hat{n}_{1},\hat{n}_{2},\cdots\hat{n}_{N})$ become distinct. Thus,
in this paper, we will require positive integers $(n_{1},n_{2},\cdots n_{N})$
to be distinct without loss of generality. In this case, we will also arrange
them in the ascending order, i.e., $n_{1}<n_{2}<\cdots<n_{N}$.

Remark 2. The higher-order lumps in Theorem 1 contain free complex parameters
$p$ and $\textbf{\emph{a}}_{i}$ ($1\leq i\leq N$), totaling
$1+n_{1}+n_{2}+\cdots n_{N}$. However, using techniques similar to that
outlined in Remark 1, we can show that $N(N-1)/2$ of those parameters in
$\\{\textbf{\emph{a}}_{i}\\}$ can be eliminated. Thus, the number of free
complex parameters in these higher-order lumps can be reduced to $1+\rho$,
where

$\rho=\sum_{i=1}^{N}n_{i}-\frac{N(N-1)}{2}.$ (13)

This number of free parameters matches that given in Ref. Ablowitz2000  for
solutions produced by Darboux transformation. In fact, from the derivation of
Theorem 1 in the Appendix, we can see that our higher-order lumps given in
this theorem by the bilinear method are equivalent to those derived in
Ablowitz2000  by Darboux transformation, except that our expressions are more
explicit.

Remark 3. The fundamental lump can be derived by taking $N=1$ and $n_{1}=1$ in
Eq. (8). Through a shift of the $(x,y)$ axes, we can normalize $a_{1,1}=0$.
Then, the resulting $\sigma_{1}(x,y,t)$ function can be reduced to

$\displaystyle\sigma_{1}=\left|x+2\textrm{i}py-12p^{2}t\right|^{2}+\frac{1}{(p+p^{*})^{2}}$
$\displaystyle\hskip 8.5359pt=\left[x-2p_{i}\hskip
0.85355pty-12(p_{r}^{2}-p_{i}^{2})\hskip
0.85355ptt\right]^{2}+\left[2p_{r}(y-12p_{i}\hskip
0.85355ptt)\right]^{2}+\frac{1}{4p_{r}^{2}},$ (14)

where $p_{r}$ and $p_{i}$ are the real and imaginary parts of the spectral
parameter $p$. The corresponding solution $u_{1}(x,y,t)$ through Eq. (7) moves
at $x$-direction velocity of $12|p|^{2}$ and $y$-direction velocity of
$12p_{i}$. By applying the Galilean invariance (4) with $k=p_{i}$, we can
remove the $y$-direction velocity $12p_{i}$ and reduce $\sigma_{1}(x,y,t)$ to

$\displaystyle\sigma_{1}=\left(x-12p_{r}^{2}t\right)^{2}+\left(2p_{r}y\right)^{2}+\frac{1}{4p_{r}^{2}}.$
(15)

This means that, under Galilean invariance, we can take $p$ in the original
fundamental lump to be purely real without loss of generality. Then, by
utilizing the scaling invariance (5) with $\alpha=p_{r}$, we can further
normalize $p_{r}$ in the above $\sigma_{1}$ to be unity. The final simplified
fundamental-lump expression is

$u_{1}(x,y,t)=2\partial_{x}^{2}\ln\left[\left(x-12t\right)^{2}+4y^{2}+\frac{1}{4}\right].$
(16)

This is a moving single lump with peak amplitude $16$, which is attained at
the spatial location of $(x,y)=(12t,\hskip 1.42271pt0)$.

Remark 4. In the general higher-order lump of Theorem 1, the whole solution
complex moves at $x$-direction velocity $12|p|^{2}$ and $y$-direction velocity
$12p_{i}$, plus some possible slower motion relative to that moving frame. In
this general case, we can also use the Galilean invariance (4) to remove the
$y$-direction velocity $12p_{i}$ of the complex, i.e., $p_{i}$ can be made to
be zero. In addition, we can use the scaling invariance (5) to normalize
$p_{r}$ to unity. Thus, without any loss of generality, we can choose $p$ in
the higher-order lump solution of Theorem 1 to be equal to one. For this
reason, we will set $p=1$ in the remainder of this paper.

Remark 5. In Peli93a , a class of higher-order lump solutions that are
stationary in a moving frame was reported. Those special solutions satisfy the
Boussinesq equation. Thus, they are special cases of Boussinesq rogue waves
YangYangBoussi . Those stationary higher-order lumps are also special cases of
our solutions in Theorem 1 when the index vector $(n_{1},n_{2},\cdots n_{N})$
and internal parameters $\\{\textbf{\emph{a}}_{i}\\}$ are properly chosen.
Indeed, rational solutions in Theorem 1 would be stationary if the $\sigma$
function in (8) satisfies the dimension-reduction condition
$\sigma_{t}-V\sigma_{x}=0$, where $V$ is the velocity of the moving frame
along the $x$-direction. In the bilinear derivation of Boussinesq rogue waves
YangYangBoussi , one needs to solve the bilinear $\tau$ equation of KP-I,
together with this $\tau$’s dimension reduction condition
$\tau_{x_{3}}-3\tau_{x_{1}}=C\tau$, where $x_{1}$ is proportional to $x$,
$x_{3}$ proportional to $t$, and $C$ is a constant. Since this $\tau$ function
turns out to be equal to $\sigma$ multiplying an exponential of a linear
function of $x$ and $t$ YangYangBoussi , we see that $\tau$’s dimension
reduction condition is equivalent to $\sigma$’s dimension reduction condition
after proper variable scalings. This means that constraints from $\tau$’s
dimension-reduction condition can be borrowed over and imposed on solutions of
Theorem 1 in order to obtain stationary higher-order KP lumps. One of such
constraints is on the index vector $(n_{1},n_{2},\cdots n_{N})$, where
$n_{i}=2i-1$ must be chosen YangYangBoussi . In addition, internal parameters
$\\{\textbf{\emph{a}}_{i}\\}$ also need to be constrained. For a different
choice of differential operators than those in Eq. (105) of the Appendix, this
parameter constraint was derived in YangYangBoussi . For the present choice of
differential operators in Eq. (105), this parameter constraint would be more
complex. In this case, such a parameter constraint was worked out in
ChenJunChao2018  for another integrable system under a different
parameterization.

###  II.2 Yablonskii–Vorob’ev polynomials and Wronskian-Hermit polynomials

We will show in later text that patterns of certain higher-order lump
solutions at large time are described by root structures of the
Yablonskii–Vorob’ev polynomials and Wronskian-Hermit polynomials. Thus, these
polynomials and their root structures will be introduced first.

####  II.2.1 Yablonskii–Vorob’ev polynomials and their root structures

Yablonskii-Vorob’ev polynomials arose in rational solutions of the second
Painlevé equation ($\mbox{P}_{\mbox{\scriptsize II}}$) Yablonskii1959 ;
Vorobev1965 . Later, a determinant expression for these polynomials was found
in Kajiwara-Ohta1996 . Let $p_{k}(z)$ be polynomials defined by

$\sum_{k=0}^{\infty}p_{k}(z)\epsilon^{k}=\exp\left(z\epsilon-\frac{4}{3}\epsilon^{3}\right).$
(17)

Then, Yablonskii-Vorob’ev polynomials $Q_{N}(z)$ are given by the $N\times N$
determinant Kajiwara-Ohta1996

$\displaystyle
Q_{N}(z)=c_{N}\left|\begin{array}[]{cccc}p_{1}(z)&p_{0}(z)&\cdots&p_{2-N}(z)\\\
p_{3}(z)&p_{2}(z)&\cdots&p_{4-N}(z)\\\ \vdots&\vdots&\vdots&\vdots\\\
p_{2N-1}(z)&p_{2N-2}(z)&\cdots&p_{N}(z)\end{array}\right|,$ (22)

where $c_{N}=\prod_{j=1}^{N}(2j-1)!!$, and $p_{k}(z)\equiv 0$ if $k<0$. This
determinant is a Wronskian since one can see from Eq. (17) that
$p^{\prime}_{k+1}(z)=p_{k}(z)$, where the prime represents differentiation.
Yablonskii-Vorob’ev polynomials are monic polynomials with integer
coefficients Clarkson2003-II , and the first four of them are

$\displaystyle Q_{1}(z)=z,$ $\displaystyle Q_{2}(z)=z^{3}+4,$ $\displaystyle
Q_{3}(z)=z^{6}+20z^{3}-80,$ $\displaystyle Q_{4}(z)=z(z^{9}+60z^{6}+11200).$

Root structures of these polynomials have been studied in Fukutani ; Taneda ;
Clarkson2003-II ; Miller2014 ; Bertola2016 , and the following facts are
known.

  1. 1.

The degree of the $Q_{N}(z)$ polynomial is $N(N+1)/2$, which can be easily
seen from Eq. (22).

  2. 2.

All roots of $Q_{N}(z)$ are simple Fukutani . Thus, $Q_{N}(z)$ has $N(N+1)/2$
simple roots.

  3. 3.

Zero is a root of $Q_{N}(z)$ if and only if $N\equiv 1\hskip
2.84544pt\mbox{mod}\hskip 2.84544pt3$ Taneda .

  4. 4.

$Q_{N}(z)$ can be factorized as $Q_{N}(z)=z^{m}f(\zeta)$, where $m=1$ when
$N\equiv 1\hskip 2.84544pt\mbox{mod}\hskip 2.84544pt3$ and $m=0$ otherwise,
$\zeta\equiv z^{3}$, and $f(\zeta)$ is a polynomial of $\zeta$ with integer
coefficients and a nonzero constant term Clarkson2003-II . This factorization
shows that the root structure of $Q_{N}(z)$ is invariant under
$120^{\circ}$-angle rotation in the complex $z$ plane.

  5. 5.

Roots of $Q_{N}(z)$ exhibit a triangular pattern in the complex plane for all
$N\geq 2$ Clarkson2003-II ; Miller2014 ; Bertola2016 . This fact is not
surprising given the $120^{\circ}$ rotational symmetry of $Q_{N}(z)$’s root
structure mentioned above.

  6. 6.

Roots of $Q_{N}(z)$ are also symmetric with respect to the real-$z$ axis,
since the coefficients of $Q_{N}(z)$ are real and thus complex roots appear in
conjugate pairs. This conjugate symmetry, together with the $120^{\circ}$
rotational symmetry, implies that one vertex of the triangular root structure
of $Q_{N}(z)$ is on the real-$z$ axis.

Due to importance of these root structures to our work, we reproduce some of
them in Fig. 1 for $2\leq N\leq 5$.

Figure 1: Root structures of Yablonskii-Vorob’ev polynomials $Q_{N}(z)$ in the
complex $z$ plane for $2\leq N\leq 5$.

####  II.2.2 Wronskian-Hermit polynomials and their root structures

Next, we introduce Wronskian-Hermit polynomials. Let $q_{k}(z)$ be polynomials
defined by

$\sum_{k=0}^{\infty}q_{k}(z)\epsilon^{k}=\exp\left(z\epsilon+\epsilon^{2}\right).$
(23)

These $q_{k}(z)$ polynomials are related to Hermit polynomials through simple
variable scalings. Then, for any positive integer $N$ and index vector
$\Lambda=(n_{1},n_{2},\dots,n_{N})$, where $\\{n_{i}\\}$ are positive and
distinct integers in ascending order, i.e., $n_{1}<n_{2}<\cdots<n_{N}$, the
Wronskian-Hermite polynomial $W_{\Lambda}(z)$ is defined as the Wronskian of
$q_{k}(z)$ polynomials

$W_{\Lambda}(z)=\mbox{Wron}\left[q_{n_{1}}(z),q_{n_{2}}(z),\ldots,q_{n_{N}}(z)\right],$
(24)

or equivalently,

$W_{\Lambda}(z)=\left|\begin{array}[]{cccc}q_{n_{1}}(z)&q_{n_{1}-1}(z)&\cdots&q_{n_{1}-N+1}(z)\\\
q_{n_{2}}(z)&q_{n_{2}-1}(z)&\cdots&q_{n_{2}-N+1}(z)\\\
\vdots&\vdots&\vdots&\vdots\\\
q_{n_{N}}(z)&q_{n_{N}-1}(z)&\cdots&q_{n_{N}-N+1}(z)\end{array}\right|,$ (25)

since we can see $q_{k+1}^{\prime}(z)=q_{k}(z)$ from the definition (23). In
the above determinant, $q_{k}(z)\equiv 0$ when $k<0$.

Regarding root structures of Wronskian-Hermite polynomials $W_{\Lambda}(z)$,
we have the following facts.

  1. 1.

The degree of the polynomial $W_{\Lambda}(z)$ is equal to $\rho$, where $\rho$
is given in Eq. (13). This fact can be seen from the definition (25).

  2. 2.

The multiplicity of the zero root in $W_{\Lambda}(z)$ is equal to $d(d+1)/2$,
where

$d=k_{odd}-k_{even},$ (26)

and $k_{odd}$, $k_{even}$ are the numbers of odd and even elements in the
index vector $(n_{1},n_{2},\dots,n_{N})$ respectively. This fact was mentioned
in Felder2012 ; Garcia_Ferrero2015  and proved in Bonneux2020 . If
$d(d+1)/2=0$, i.e., $d=0$ or $-1$, then zero is not a root of
$W_{\Lambda}(z)$.

  3. 3.

The number of nonzero roots (counting multiplicity) in $W_{\Lambda}(z)$, which
we denote as $N_{W}$, is

$N_{W}=\rho-\frac{d(d+1)}{2}.$ (27)

  4. 4.

The polynomial $W_{\Lambda}(z)$ can be factored as
$W_{\Lambda}(z)=z^{d(d+1)/2}f(\zeta)$, where $d$ is given in Eq. (26),
$\zeta\equiv z^{2}$, and $f(\zeta)$ is a polynomial of $\zeta$ with real
coefficients and a nonzero constant term Bonneux2020 .

  5. 5.

If $z_{0}$ is a root of $W_{\Lambda}(z)$, so are $-z_{0},z_{0}^{*}$ and
$-z_{0}^{*}$. This quartet root symmetry can be seen from the above
factorization of $W_{\Lambda}(z)$ and the fact that the coefficients of the
polynomial $W_{\Lambda}(z)$ are real. As a consequence of this quartet
symmetry, the root structure of $W_{\Lambda}(z)$ is non-triangular. This is a
big difference from Yablonskii-Vorob’ev polynomials, which feature triangular
root structures.

In addition, we have the following lemma.

> Lemma 1. The Wronskian-Hermite polynomial $W_{\Lambda}(z)$ has only zero
> roots, i.e., $N_{W}=0$, if and only if
> $(n_{1},n_{2},\dots,n_{N})=(1,3,5,\cdots,2N-1)$.

Proof. Since $k_{odd}+k_{even}=N$, we have from Eq. (27) that

$N_{W}=\sum_{i=1}^{N}n_{i}-k_{odd}^{2}-k_{even}(k_{even}-1).$ (28)

Since $\\{n_{i}\\}$ are distinct positive integers, their smallest possible
values, after reordering, are $\\{1,3,\cdots
2k_{odd}-1,2,4,\cdots,2k_{even}\\}$. Thus,

$\sum_{i=1}^{N}n_{i}\geq\left[1+3+\cdots(2k_{odd}-1)\right]+\left(2+4+\cdots+2k_{even}\right)=k_{odd}^{2}+k_{even}(k_{even}+1).$
(29)

Then,

$N_{W}\geq 2k_{even}.$ (30)

For $W_{\Lambda}(z)$ to have only zero roots, $N_{W}$ must be zero; so
$k_{even}=0$, i.e., all numbers in $\\{n_{i}\\}$ must be odd. In addition, for
the equality in (30) to hold, all these odd and distinct numbers must be the
lowest, i.e., $(n_{1},n_{2},\dots,n_{N})=(1,3,5,\cdots,2N-1)$. This completes
the proof.

This lemma tells us that, as long as $\Lambda\neq(1,3,5,\cdots,2N-1)$, the
Wronskian-Hermit polynomial $W_{\Lambda}(z)$ would always have nonzero roots.
This result is important to us, as we will show in later text that the
presence or absence of nonzero roots in $W_{\Lambda}(z)$ will have direct
consequences on the solution patterns of higher-order lumps.

On roots of Wronskian-Hermite polynomials, beside the above facts, the
following conjecture has also been proposed.

> Conjecture 1 Felder2012 . All roots of every Wronskian-Hermite polynomial
> $W_{\Lambda}(z)$ are simple, except possibly the zero root.

This conjecture will be useful, as we will show in later text that the
multiplicity of a root in the Wronskian-Hermite polynomial has direct
implications on the wave structure of higher-order lumps. Based on this
conjecture, $W_{\Lambda}(z)$ would have $N_{W}$ nonzero simple roots, where
$N_{W}$ is given in Eq. (27). We have checked this conjecture on a number of
examples of $W_{\Lambda}(z)$, and found it to always hold.

To illustrate root structures of Wronskian-Hermite polynomials, we choose two
index vectors

$\Lambda_{1}=(2,3,4,5),\quad\Lambda_{2}=(3,4,5,7,9).$ (31)

The corresponding polynomials are found to be

$\displaystyle W_{\Lambda_{1}}(z)=\frac{z^{8}-16z^{6}+120z^{4}+720}{2880},$
(32) $\displaystyle
W_{\Lambda_{2}}(z)=-\frac{z^{6}\left(z^{12}-12z^{10}+180z^{8}+672z^{6}-7056z^{4}-181440z^{2}-1270080\right)}{2743372800}.$
(33)

Root structures of these two polynomials are plotted in Fig. 2. It is seen
that for the first polynomial, its root structure is rectangular and does not
contain zero. For the second polynomial, its root structure is quasi-
rectangular with a zero root (of multiplicity six) in the center. All nonzero
roots in these two polynomials are simple, which is consistent with the
earlier conjecture.

Figure 2: Root structures of Wronskian-Hermite polynomials
$W_{\Lambda_{1}}(z)$ (left) and $W_{\Lambda_{2}}(z)$ (right) in the complex
plane $z$, where index vectors $\Lambda_{1}$ and $\Lambda_{2}$ are given in
Eq. (31).

##  III Patterns of higher-order lumps at large times

In this section, we study patterns of higher-order lumps at large times. In
this study, we will set the spectral parameter $p=1$ without loss of
generality (see Remark 4 in the previous section). In this case, the constant
factor in Eq. (9) simplifies to $1/4^{\nu}$. In addition, the definition (11)
of the s vector reduces to

$\displaystyle\ln\left[\frac{2}{\kappa}\tanh\left(\frac{\kappa}{2}\right)\right]=\sum_{k=1}^{\infty}s_{k}\hskip
1.42271pt\kappa^{k},$ (34)

which is identical to the s vector in the earlier work OhtaJY2012  on rogue
waves of the NLS equation. In particular, all odd-indexed elements $s_{odd}$
of s are zero since the function on the left side of the above equation is
even.

It turns out that pattern analysis of lumps depends on whether vectors of
internal parameters $\\{\textbf{\emph{a}}_{i}\\}$ are the same vectors (i.e.,
whether vector elements $a_{i,j}$ depend on the $i$ index). In this paper, we
only consider the case where these $\\{\textbf{\emph{a}}_{i}\\}$ vectors are
the same, i.e., $\textbf{\emph{a}}_{i}=\textbf{\emph{a}}$. In this case,
through a shift of the $(x,y)$ axes, we can make the first element of a to be
zero. Thus, our parameter choices are

$\textbf{\emph{a}}_{i}=\textbf{\emph{a}}=(0,a_{2},a_{3},\cdots).$ (35)

Under these parameters, we have two theorems on patterns of higher-order lumps
at large times, depending on whether the index vector $\Lambda$ is equal to
$(1,3,5,\dots,2N-1)$.

###  III.1 Large-time lump patterns when $\Lambda=(1,3,5,\dots,2N-1)$

Our first theorem is for the case where the index vector $\Lambda$ is equal to
$(1,3,5,\dots,2N-1)$.

> Theorem 2. If the index vector $\Lambda=(1,3,5,\dots,2N-1)$, then, when
> $|t|\gg 1$, the higher-order lump solution $u_{\Lambda}(x,y,t)$
> asymptotically separates into $N(N+1)/2$ fundamental lumps
> $u_{1}(x-x_{0},\hskip 1.13791pty-y_{0},t)$, where $u_{1}(x,y,t)$ is given in
> Eq. (16),
>
> $x_{0}=\Re(z_{0})\hskip 1.42271pt(12t)^{1/3},\quad
> y_{0}=\frac{\Im(z_{0})}{2}(12t)^{1/3},$ (36)
>
> $z_{0}$ is each of the $N(N+1)/2$ simple roots of the Yablonskii–Vorob’ev
> polynomial $Q_{N}(z)$, and $\Re$, $\Im$ represent the real and imaginary
> parts of a complex number. The peak of each fundamental lump is spatially
> located at $(x,y)=(12t+x_{0},y_{0})$. The absolute error of this
> fundamental-lump approximation is $O(|t|^{-1/3})$ when $z_{0}\neq 0$ and
> $O(|t|^{-1})$ when zero is a root and $z_{0}=0$. Expressed mathematically,
> when $(x,y)$ is in the neighborhood of each of these fundamental lumps,
> i.e., $(x-12t-x_{0})^{2}+(y-y_{0})^{2}=O(1)$, we have the following solution
> asymptotics for $|t|\gg 1$,
>
> $u_{\Lambda}(x,y,t)=\left\\{\begin{array}[]{ll}u_{1}(x-x_{0},\hskip
> 1.13791pty-y_{0},t)+O\left(|t|^{-1/3}\right),&\mbox{if}\hskip
> 1.9919ptz_{0}\neq 0,\\\ u_{1}(x-x_{0},\hskip 1.13791pty-
> y_{0},t)+O\left(|t|^{-1}\right),&\mbox{if}\hskip
> 1.9919ptz_{0}=0.\end{array}\right.$ (37)
>
> When $(x,y)$ is not in the neighborhood of any of these $N(N+1)/2$
> fundamental lumps, $u_{\Lambda}(x,y,t)$ asymptotically approaches zero as
> $|t|\to\infty$.

The proof of this theorem will be provided in Sec. V.

This theorem indicates that, wave patterns at large times are formed by
$N(N+1)/2$ fundamental lumps. Relative to the moving frame of $x$-direction
velocity $12$, positions $(x_{0},y_{0})$ of these fundamental lumps are just a
simple linear transformation of the root structure of the Yablonskii–Vorob’ev
polynomial $Q_{N}(z)$, i.e.,

$\left[\begin{array}[]{c}x_{0}\\\
y_{0}\end{array}\right]=(12t)^{1/3}\left[\begin{array}[]{cc}1&0\\\
0&\frac{1}{2}\end{array}\right]\left[\begin{array}[]{c}\Re(z_{0})\\\
\Im(z_{0})\end{array}\right].$ (38)

Since the transformation matrix is diagonal, this transformation is simply a
stretching along both horizontal and vertical directions. Recall that the
Yablonskii–Vorob’ev root structure is triangular (see Fig. 1). The resulting
lump pattern is then triangular as well. When $t\gg 1$, this triangular lump
pattern preserves the same orientation of the original triangle of the
Yablonskii–Vorob’ev root structure. But when $t\ll-1$, the triangular lump
pattern would be oriented opposite of the Yablonskii–Vorob’ev root structure.
Indeed, it is easy to see from Eq. (38) that, when time changes from large
negative to large positive, i.e., from $-t$ to $+t$, their lump positions
would be related as

$\left[\begin{array}[]{c}x_{0}^{+}\\\
y_{0}^{+}\end{array}\right]=-\left[\begin{array}[]{c}x_{0}^{-}\\\
y_{0}^{-}\end{array}\right].$ (39)

Thus, these triangular lump patterns have reversed directions along the
$x$-axis (the $y$-direction reversal does not matter since the pattern is
symmetric in $y$). This $x$-direction reversal of triangular lump patterns
when time changes from large negative to large positive is a dramatic pattern
transformation in the KP-I equation. This phenomenon has been graphically
reported in Dubard2013  on several low-order solution examples. Here, we
established this fact for the general case.

Theorem 2 also indicates that, at large time, fundamental lumps in the
solution complex separate from each other in proportion to $|t|^{1/3}$. This
rate of separation is very slow, relative to the overall (linear) speed $12$
of the whole complex.

One more feature of Theorem 2 is that, positions $(x_{0},y_{0})$ in Eq. (36)
for individual fundamental lumps in the solution complex are independent of
the solution’s internal parameters a. This means that, when $|t|\to\infty$,
solutions $u_{\Lambda}(x,y,t)$ with different internal parameters a would
approach the same limit solution.

###  III.2 Large-time lump patterns when $\Lambda\neq(1,3,5,\dots,2N-1)$

When $\Lambda\neq(1,3,5,\dots,2N-1)$, the Wronskian-Hermite polynomial
$W_{\Lambda}(z)$ has a zero root of multiplicity $d(d+1)/2$, with $d$ given in
Eq. (26), as well as nonzero roots that are conjectured to be all simple (see
Sec. II.2.2). Note that the zero root would be absent if $d=0$ or $-1$; but
nonzero roots always exist. In this case, our results on solution patterns at
large time are summarized in the following theorem.

> Theorem 3.  Suppose the index vector $\Lambda\neq(1,3,5,\dots,2N-1)$, and
> all nonzero roots of the Wronskian-Hermite polynomial $W_{\Lambda}(z)$ are
> simple. Then, for $|t|\gg 1$, the following asymptotics for the solution
> $u_{\Lambda}(x,y,t)$ holds.
>
>   1. 1.
>
> In the outer region — the region that is $O\left(|t|^{1/2}\right)$ away from
> the wave center $(x,y)=(12t,0)$, or
> $\sqrt{(x-12t)^{2}+y^{2}}=O\left(|t|^{1/2}\right)$, the higher-order lump
> $u_{\Lambda}(x,y,t)$ asymptotically separates into $N_{W}$ fundamental lumps
> $u_{1}(x-x_{0},\hskip 1.13791pty-y_{0},t)$, where $N_{W}$ is given in Eq.
> (27), $u_{1}(x,y,t)$ is given in Eq. (16),
>
> $x_{0}=\Re\left[z_{0}(-12t)^{1/2}\right]+O(1),\quad
> y_{0}=\frac{\Im\left[z_{0}(-12t)^{1/2}\right]}{2}+O(1),$ (40)
>
> and $z_{0}$ is each of the $N_{W}$ nonzero simple roots of $W_{\Lambda}(z)$.
> The absolute error of this fundamental-lump approximation is
> $O(|t|^{-1/2})$. Expressed mathematically, when $(x,y)$ is in the
> neighborhood of each of these fundamental lumps, i.e.,
> $(x-12t-x_{0})^{2}+(y-y_{0})^{2}=O(1)$, we have the following solution
> asymptotics
>
> $u_{\Lambda}(x,y,t)=u_{1}(x-x_{0},\hskip 1.13791pty-
> y_{0},t)+O\left(|t|^{-1/2}\right),\quad|t|\gg 1.$ (41)
>
>   2. 2.
>
> If zero is a root of $W_{\Lambda}(z)$, i.e., $d\neq 0$ and $-1$, then in the
> inner region — the region that is within $O(|t|^{1/3})$ of the wave center
> $(x,y)=(12t,0)$, or $\sqrt{(x-12t)^{2}+y^{2}}\leq O\left(|t|^{1/3}\right)$,
> lies $d(d+1)/2$ fundamental lumps $u_{1}(x-x_{0},\hskip 1.13791pty-
> y_{0},t)$, where $u_{1}(x,y,t)$ is given in Eq. (16),
>
> $x_{0}=\Re(z_{0})\hskip 1.42271pt(12t)^{1/3}+O(1),\quad
> y_{0}=\frac{\Im(z_{0})}{2}(12t)^{1/3}+O(1),$ (42)
>
> and $z_{0}$ is each of the $d(d+1)/2$ simple roots of the
> Yablonskii–Vorob’ev polynomial $Q_{\hat{d}}(z)$, with $\hat{d}$ defined as
>
> $\hat{d}=\left\\{\begin{array}[]{ll}d,&\mbox{when}\hskip 2.84544ptd\geq
> 0,\\\ |d|-1,&\mbox{when}\hskip 2.84544ptd\leq-1.\end{array}\right.$ (43)
>
> Notice that $d(d+1)/2=\hat{d}(\hat{d}+1)/2$. The absolute error of this
> fundamental-lump approximation is $O(|t|^{-1/3})$ when $z_{0}\neq 0$ and
> $O(|t|^{-1})$ when zero is a root of $Q_{\hat{d}}(z)$ and $z_{0}=0$.
> Expressed mathematically, when $(x,y)$ is in the neighborhood of each of
> these fundamental lumps, i.e., $(x-12t-x_{0})^{2}+(y-y_{0})^{2}=O(1)$, with
> $(x_{0},y_{0})$ given in (42), we have the following solution asymptotics
> for $|t|\gg 1$,
>
> $u_{\Lambda}(x,y,t)=\left\\{\begin{array}[]{ll}u_{1}(x-x_{0},\hskip
> 1.13791pty-y_{0},t)+O\left(|t|^{-1/3}\right),&\mbox{if}\hskip
> 1.9919ptz_{0}\neq 0,\\\ u_{1}(x-x_{0},\hskip 1.13791pty-
> y_{0},t)+O\left(|t|^{-1}\right),&\mbox{if}\hskip
> 1.9919ptz_{0}=0.\end{array}\right.$ (44)
>
>   3. 3.
>
> When $(x,y)$ is not in the neighborhood of any of the above fundamental
> lumps, $u_{\Lambda}(x,y,t)$ asymptotically approaches zero as
> $|t|\to\infty$.
>
>

Remark 6. In this theorem, we assumed all nonzero roots of $W_{\Lambda}(z)$
simple, which is true for all examples we tested, such as the two in Eqs.
(32)-(33). In view of Conjecture 1 in the previous section, this assumption is
expected to hold in all cases. If this conjecture is false, i.e., some nonzero
roots of $W_{\Lambda}(z)$ are not simple, then this theorem for the outer
region, i.e., Eqs. (40)-(41), would still hold, but only for nonzero simple
roots $z_{0}$ of $W_{\Lambda}(z)$.

Now, we explain what Theorem 3 says regarding solution patterns at large times
when $\Lambda\neq(1,3,5,\dots,2N-1)$. In this case, Theorem 3 indicates that,
the whole wave field is generically split up into two regions featuring
different patterns.

  1. 1.

In the outer region — the region that is $O(|t|^{1/2})$ away from the wave
center $(x,y)=(12t,0)$, the wave field at large time comprises $N_{W}$
fundamental lumps, whose positions are given through the nonzero roots of the
Wronskian-Hermite polynomial $W_{\Lambda}(z)$. Specifically, relative to the
moving frame of $x$-direction velocity $12$, positions $(x_{0},y_{0})$ of
these fundamental lumps, to the leading order of large time, are just a linear
transformation of $W_{\Lambda}(z)$’s nonzero-root structure. The reader is
reminded from Sec. II.2.2 that when $\Lambda\neq(1,3,5,\dots,2N-1)$, nonzero
roots of $W_{\Lambda}(z)$ always exist, and their shape in the $z$-plane is
non-triangular. When $t$ is large negative, these fundamental-lump positions
to the leading order are

$\left[\begin{array}[]{c}x_{0}^{-}\\\
y_{0}^{-}\end{array}\right]=(12|t|)^{1/2}\left[\begin{array}[]{cc}1&0\\\
0&\frac{1}{2}\end{array}\right]\left[\begin{array}[]{c}\Re(z_{0})\\\
\Im(z_{0})\end{array}\right],$ (45)

where $z_{0}$ is any nonzero root of $W_{\Lambda}(z)$. However, when $t$ is
large positive, these lump positions become

$\left[\begin{array}[]{c}x_{0}^{+}\\\
y_{0}^{+}\end{array}\right]=(12|t|)^{1/2}\left[\begin{array}[]{cc}0&-1\\\
\frac{1}{2}&0\end{array}\right]\left[\begin{array}[]{c}\Re(z_{0})\\\
\Im(z_{0})\end{array}\right].$ (46)

In the former case, the wave pattern formed by these fundamental lumps is
simply a stretching of the Wronskian-Hermite nonzero-root structure along both
horizontal and vertical directions. But in the latter case, on top of this
stretching, the horizontal and vertical directions are also swapped. In both
cases, the resulting wave patterns from transformations (45)-(46) are non-
triangular since the root structure of $W_{\Lambda}(z)$ is non-triangular.

From the above two transformations, we see that fundamental lumps at large
negative time $-t$ and large positive time $+t$ in the outer region are
related as

$\left[\begin{array}[]{c}x_{0}^{+}\\\
y_{0}^{+}\end{array}\right]=\left[\begin{array}[]{cc}0&-2\\\
\frac{1}{2}&0\end{array}\right]\left[\begin{array}[]{c}x_{0}^{-}\\\
y_{0}^{-}\end{array}\right].$ (47)

Thus, when time goes from large negative to large positive, outer-region lump
patterns in the $(x,y)$ plane have swapped horizontal and vertical directions.
In addition, stretching of different amounts has also occurred along the two
directions. This swapping of horizontal and vertical directions is another
type of dramatic pattern transformation, and it is very different from the
triangular $x$-direction reversal that occurs when
$\Lambda=(1,3,5,\dots,2N-1)$. For certain single-line patterns of fundamental
lumps, a change from a vertical line to a horizontal line in the $(x,y)$ plane
has been graphically reported in Chen2016  and analytically explained in
Chang2018 . Here, we proved this fact for the general case, where patterns of
fundamental lumps based on Wronskian-Hermite root structures can be arbitrary,
not just lines (see the next section for examples).

In this outer region, fundamental lumps at large time separate from each other
in proportion to $|t|^{1/2}$. This is another big difference between the
present solutions and those with $\Lambda=(1,3,5,\dots,2N-1)$ in the previous
subsection, where lumps separate in proportion to $|t|^{1/3}$ instead.

  2. 2.

In the inner region — the region that is within $O(|t|^{1/3})$ of the wave
center $(x,y)=(12t,0)$, if $d\neq 0$ and $-1$, then the solution
$u_{\Lambda}(x,y,t)$ at large time would comprise $d(d+1)/2$ fundamental
lumps, whose positions are given through roots of the Yablonskii–Vorob’ev
polynomial $Q_{\hat{d}}(z)$, with $\hat{d}$ defined in Eq. (43). The reader is
reminded that $\hat{d}(\hat{d}+1)/2=d(d+1)/2$. Relative to the moving frame of
$x$-direction velocity $12$, positions $(x_{0},y_{0})$ of these fundamental
lumps, to the leading order of large time, are just a linear transformation of
$Q_{\hat{d}}(z)$’s root structure, i.e.,

$\left[\begin{array}[]{c}x_{0}\\\
y_{0}\end{array}\right]=(12t)^{1/3}\left[\begin{array}[]{cc}1&0\\\
0&\frac{1}{2}\end{array}\right]\left[\begin{array}[]{c}\Re(z_{0})\\\
\Im(z_{0})\end{array}\right],$ (48)

where $z_{0}$ is each of the $d(d+1)/2$ simple roots of $Q_{\hat{d}}(z)$. This
lump-position formula in the inner region is very similar to (38) of the
$\Lambda=(1,3,5,\dots,2N-1)$ case. Thus, the pattern of these $d(d+1)/2$
fundamental lumps in the inner region at large time is a simple stretching of
$Q_{\hat{d}}(z)$’s root structure, and the resulting pattern is triangular if
$\hat{d}>1$. In addition, as time evolves from large negative to large
positive, these triangular lump patterns would reverse direction along the
$x$-axis. Furthermore, fundamental lumps in this inner region separate from
each other in proportion to $|t|^{1/3}$ at large time, similar to the
$\Lambda=(1,3,5,\dots,2N-1)$ case in Theorem 2. If $\hat{d}=0$, i.e., $d=0$ or
$-1$, this inner region would be absent.

The above results reveal that, the pattern of the solution
$u_{\Lambda}(x,y,t)$ for $\Lambda\neq(1,3,5,\dots,2N-1)$ at large time is
richer, with the outer region exhibiting the non-triangular shape of the
stretched nonzero-root structure of the Wronskian-Hermite polynomial
$W_{\Lambda}(z)$, and with the inner region exhibiting the triangular shape of
the stretched root structure of the Yablonskii–Vorob’ev polynomial
$Q_{\hat{d}}(z)$. As time changes from large negative to large positive, the
outer pattern swaps horizontal and vertical directions, while the inner
pattern reverses the horizontal direction. These different types of pattern
transformations in the outer and inner regions of the same solution are
fascinating. When $\Lambda\neq(1,3,5,\dots,2N-1)$, the outer pattern is always
present since $W_{\Lambda}(z)$ always has nonzero roots, but the inner pattern
is present only when $d\neq 0$ and $-1$ and absent otherwise. When
$\Lambda=(1,3,5,\dots,2N-1)$, the outer pattern disappears, since
$W_{\Lambda}(z)$ has only zero roots (see Lemma 1). In this special case, our
results for the inner region in Theorem 3 are consistent with those in Theorem
2. However, Theorem 2 for this special case is stronger, since it shows that
positions $(x_{0},y_{0})$ of fundamental lumps now have no $O(1)$ shifts [see
Eq. (36)] — a more accurate prediction than Eq. (42) of Theorem 3 which shows
$O(1)$ shifts in general.

In the end, we note that in the earlier work Ablowitz2000 , it was reported
that at large time, fundamental lumps in the higher-order lump complex
separate from each other in proportion to $|t|^{q}$, where $\frac{1}{3}\leq
q\leq\frac{1}{2}$. Our results in Theorems 2 and 3 indicate that this $q$
value can only be $1/3$ or $1/2$, nothing in between.

##  IV Comparison between true lump patterns and analytical predictions

In this section, we compare our analytical predictions of lump patterns with
true solutions.

###  IV.1 Pattern transformation when $\Lambda=(1,3,5,\dots,2N-1)$

First, we do the comparison when $\Lambda=(1,3,5,\dots,2N-1)$, where a
triangular pattern of lumps at large time is predicted. To be specific, we
take $N=4$; so $\Lambda=(1,3,5,7)$. Root structure of the corresponding
Yablonskii-Vorob’ev polynomial $Q_{4}(z)$ has been displayed in Fig. 1. Using
those roots and formulae (36), predicted solutions from Theorem 2 at large
times $t=-10$ and $10$ are plotted in Fig. 3.

Figure 3: Predicted solutions $u_{\Lambda}(x,y,t)$ with $\Lambda=(1,3,5,7)$ at
time values of $t=-10$ (left) and $t=10$ (right).

Now, we compare these predicted solutions with true ones. In the true solution
$u_{\Lambda}(x,y,t)$, we select its internal parameters as
$\textbf{\emph{a}}=(0,0,0,0,0,0,0)$. Then, evolutions of this true solution,
at six time values of $t=-10,\ -1,\ 0,\ 0.2,\ 1$ and $10$, are plotted in Fig.
4. When comparing these true solutions at large times $t=\pm 10$ to those
predicted in Fig. 3, they clearly match each other. First, the true solutions
at $t=\pm 10$ indeed exhibit a triangular pattern, as the prediction says.
Second, the triangular pattern at $t=10$ is indeed a $x$-direction reversal of
the triangular pattern at $t=-10$, relative to a frame moving with
$x$-direction velocity $12$. Thirdly, we have quantitatively compared the
predicted and true solutions at $t=\pm 10$, and found them to match each other
as well. For example, we have quantitatively compared the difference between
the true peak location of individual lumps and its analytical prediction (36),
similar to what we did in Figs. 5-6 of Ref. YangYang21a  on the error analysis
of rogue-pattern predictions. This comparison shows that this difference is
indeed $O(|t|^{-1/3})$ when $z_{0}\neq 0$ and $O(|t|^{-1})$ when $z_{0}=0$, as
our analytical formula (37) says. Thus, our asymptotic theory on patterns of
higher-order lumps at large times for $\Lambda=(1,3,5,\dots,2N-1)$ is fully
confirmed.

In addition to large times, Fig. 4 also displays the true solution
$u_{\Lambda}(x,y,t)$ at intermediate times, where our asymptotic theory does
not apply. These intermediate panels shed light on how the dramatic
$x$-direction reversal of triangular patterns takes place as time changes from
large negative to large positive. We see that in this solution, as time
increases from $-10$ to 10, the triangle of fundamental lumps first approach
each other and shrink in size, then coalesce at $t=0$ and form a single lump
of extreme height that is ten times that of original fundamental lumps, and
then separate into a triangle of fundamental lumps again but with reversed
$x$-direction. This transformation process is fascinating.

Figure 4: The true solution $u_{\Lambda}(x,y,t)$ with $\Lambda=(1,3,5,7)$ and
$\textbf{\emph{a}}=(0,0,0,0,0,0,0)$ at various time values shown inside the
panels.

How will this $u_{\Lambda}(x,y,t)$ solution evolve if its internal parameters
a are different from $(0,0,0,0,0,0,0)$? Theorem 2 tells us that in this case,
the $u_{\Lambda}(x,y,t)$ solution would approach the same asymptotic state as
that shown in Fig. 4 at large times. At intermediate times, however, this
$u_{\Lambda}(x,y,t)$ solution could look very different from that in Fig. 4.
For instance, by suitably choosing the a values, we can get
$u_{\Lambda}(x,y,t)$ solutions whose graphs at $t=0$ exhibit very different
patterns such as a pentagon or a heptagon — a phenomenon that has been
reported in Gaillard2018 . Thus, although these $u_{\Lambda}(x,y,t)$ solutions
with different a values exhibit the same large-time triangular patterns, how
this triangular pattern at large negative time transforms to its $x$-reversed
pattern at large positive time is a process that strongly depends on the
choices of the internal a values.

###  IV.2 Pattern transformation when $\Lambda\neq(1,3,5,\dots,2N-1)$

Next, we perform the comparison when $\Lambda\neq(1,3,5,\dots,2N-1)$, where
the solution pattern at large time is determined by nonzero-root structure of
the Wronskian-Hermit polynomial $W_{\Lambda}(z)$ in the outer region, and by
root structure of the Yablonskii-Vorob’ev polynomial $Q_{\hat{d}}(z)$ in the
inner region (if $\hat{d}>0$). Since this inner region can be present or
absent depending on the $d$ value [see Eq. (26)], we will present two
examples, one for each case.

Our first example is $N=4$ and $\Lambda=(2,3,4,5)$. In this case, $d=0$, and
thus zero is not a root of $W_{\Lambda}(z)$ and the inner region is absent.
Root structure of the corresponding Wronskian-Hermit polynomial has been
displayed in Fig. 2 (the left panel). It was seen that this $W_{\Lambda}(z)$
admits eight simple nonzero roots which form a rectangle pattern. Using those
roots and leading-order terms in formulae (40), predicted solutions from
Theorem 3 at large times $t=-6$ and $6$ are plotted in Fig. 5. The predicted
patterns contain eight fundamental lumps which also form a rectangular shape
in the $(x,y)$ plane. At $t=-6$, this lump pattern is just a stretching of the
Wronskian-Hermit root structure. But at $t=6$, this lump pattern has swapped
its $x$ and $y$ directions and changed from its original $x$-direction
orientation to the new $y$-direction orientation.

Figure 5: Predicted solutions $u_{\Lambda}(x,y,t)$ with $\Lambda=(2,3,4,5)$ at
time values of $t=-6$ (left) and $t=6$ (right).

To confirm these asymptotic predictions, we plot in Fig. 6 the corresponding
true solution $u_{\Lambda}(x,y,t)$ at six time values of $t=-6,-2,-0.5,0,2$
and $6$. In this true solution, we have selected its internal parameters as
$\textbf{\emph{a}}=(0,0,0,0,800)$. It is seen that at large times of $t=\pm
6$, the true solutions indeed comprise eight fundamental lumps forming a
rectangular shape, and their orientations have changed from the $x$-direction
to the $y$-direction, exactly as our asymptotic theory has predicted. In
addition, quantitative comparisons between these true rectangular patterns and
the predicted ones in Fig. 5 show good agreement. Thus, our asymptotic theory
on patterns of higher-order lumps at large times is fully confirmed for
$\Lambda=(2,3,4,5)$.

Figure 6: True solutions $u_{\Lambda}(x,y,t)$ with $\Lambda=(2,3,4,5)$ and
$\textbf{\emph{a}}=(0,0,0,0,800)$, at various times whose values are shown
inside the panels.

By inspecting Fig. 6, we can also see how this dramatic rectangular-pattern
reorientation takes place as time increases. First, these eight fundamental
lumps of rectangular shape with $x$-direction orientation get closer to each
other and rearrange their shapes. At $t=0$, the solution has evolved into a
pentagon of five fundamental lumps surrounding a higher-peak lump near the
center. Afterwards, this pentagon structure further adjusts its shape in
significant ways, until eight new fundamental lumps emerge as a rectangular
with $y$-direction orientation in the end. Again, this transformation process
is amazing.

Our second example is $N=5$ and $\Lambda=(3,4,5,7,9)$. In this case, $d=3$,
and thus zero is a root of multiplicity six in $W_{\Lambda}(z)$, and the inner
region is present. Root structure of the corresponding Wronskian-Hermit
polynomial has been displayed in Fig. 2 (the right panel). It is seen that
this $W_{\Lambda}(z)$ admits 12 simple nonzero roots which form a quasi-
rectangular shape, plus the zero root of multiplicity six at the center of the
quasi-rectangle. Using those roots and leading-order terms in formulae (40)
and (42), predicted solutions from Theorem 3 at large times $t=-10$ and $10$
are plotted in Fig. 7. The predicted patterns contain 12 fundamental lumps
which also form a quasi-rectangular pattern in the outer region of the $(x,y)$
plane, plus six fundamental lumps which form a triangle in the inner region.
At $t=-10$, the outer lump pattern is a stretching of the Wronskian-Hermit
polynomial $W_{\Lambda}(z)$’s nonzero-root structure, while the inner lump
pattern is a stretching of the Yablonskii–Vorob’ev polynomial $Q_{3}(z)$’s
root structure. At $t=10$, however, the predicted outer lump pattern has
rotated by $90^{\circ}$ from its $t=-10$ state [plus additional
$(x,y)$-direction stretching], while the predicted inner triangular lump
pattern has reversed its direction along the $x$-axis.

Figure 7: Predicted solutions $u_{\Lambda}(x,y,t)$ with $\Lambda=(3,4,5,7,9)$
at time values of $t=-10$ (left) and $t=10$ (right).

To confirm these asymptotic predictions, we plot in Fig. 8 the corresponding
true solution $u_{\Lambda}(x,y,t)$ at six time values of $t=-10,-2,-0.2,0,2$
and $10$. In this true solution, we have selected all-zero internal parameters
of $\textbf{\emph{a}}=(0,0,0,0,0,0,0,0,0)$. It is seen that at large times of
$t=\pm 10$, the true solutions closely match our predictions in the previous
figure. Specifically, the true solutions at these large times also split into
outer and inner regions, with outer quasi-rectangular patterns and inner
triangular patterns closely resembling our predicted ones in Fig. 7.
Quantitative comparisons between these true patterns and predicted ones show
good agreement as well. Thus, our asymptotic theory on higher-order lump
patterns at large times is fully confirmed for $\Lambda=(3,4,5,7,9)$.

Figure 8: True solutions $u_{\Lambda}(x,y,t)$ with $\Lambda=(3,4,5,7,9)$ and
$\textbf{\emph{a}}=(0,0,0,0,0,0,0,0,0)$, at various times whose values are
shown inside the panels.

True solution graphs at intermediate time values in Fig. 8 reveal how these
striking pattern transformations in outer and inner regions take place. It is
seen that all fundamental lumps in the inner and outer regions at large
negative time first move toward each other. Then they merge and coalesce at
$t\approx 0$. Afterwards, all these fundamental lumps re-emerge and move away
from each other, but not returning to their pre-merging state. Instead, the
quasi-rectangular outer lumps have swapped their $x$ and $y$ directions, and
the triangular inner lumps have reversed the $x$-direction. These pattern
transformations are visually miraculous and mysterious. But due to our Theorem
3, they can now be completely understood from a mathematical point of view.

##  V Proofs of the two theorems

Now, we prove our two theorems stated in Sec. 3. The reader is reminded that
in these proofs, $p=1$ and $\textbf{\emph{a}}_{i}$ are chosen as (35) in the
higher-order lump solutions of Theorem 1, for reasons which have been
explained earlier in the paper. Thus, solution expressions in Theorem 1 can be
simplified. Notably, the constant factor in Eq. (9) simplifies to $1/4^{\nu}$,
and the s vector is real with $s_{odd}=0$ (see the beginning of Sec. 3).

###  V.1 Proof of Theorem 2

In this case, $\Lambda=(1,3,5,\dots,2N-1)$. First, we rewrite the determinant
(8) as a larger $3N\times 3N$ determinant OhtaJY2012 ; YangYang21a

$\sigma=\left|\begin{array}[]{cc}\textbf{O}_{N\times N}&\Phi_{N\times 2N}\\\
-\Psi_{2N\times N}&\textbf{I}_{2N\times 2N}\end{array}\right|,$ (49)

where

$\Phi_{i,j}=2^{-(j-1)}S_{2i-j}\left[\textbf{\emph{x}}^{+}+(j-1)\textbf{\emph{s}}+\textbf{\emph{a}}\right],\quad\Psi_{i,j}=2^{-(i-1)}S_{2j-i}\left[(\textbf{\emph{x}}^{+})^{*}+(i-1)\textbf{\emph{s}}+\textbf{\emph{a}}^{*}\right],$
(50)

$S_{j}\equiv 0$ for $j<0$, and vectors $\textbf{\emph{x}}^{+}$ and s are given
in Eqs. (10) and (34). This determinant can be further simplified. Indeed,
using the technique outlined in Appendix A of Ref. YangYang21a , we can
eliminate all $x_{even}^{+}$ and $a_{even}$ terms from the vectors
$\textbf{\emph{x}}^{+}$ and a, and reduce the above matrix element formulae to

$\Phi_{i,j}=2^{-(j-1)}S_{2i-j}\left[\hat{\textbf{\emph{x}}}^{+}+(j-1)\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}\right],\quad\Psi_{i,j}=2^{-(i-1)}S_{2j-i}\left[(\hat{\textbf{\emph{x}}}^{+})^{*}+(i-1)\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}^{*}\right],$
(51)

where

$\hat{\textbf{\emph{x}}}^{+}\equiv\left(x_{1}^{+},0,x_{3}^{+},0,x_{5}^{+},0,\cdots\right),\quad\hat{\textbf{\emph{a}}}\equiv\left(0,0,a_{3},0,a_{5},0,\cdots\right).$
(52)

The elimination of the solution’s dependence on $x_{2}^{+}$ is a key feature
of the index vector $\Lambda=(1,3,5,\dots,2N-1)$, and this feature is
responsible for the distinctive pattern behaviors described in Theorem 2.

Now, we analyze the large-time asymptotics of the above determinant $\sigma$.
For this purpose, we introduce a moving $x$-frame coordinate

$\hat{x}\equiv x-12t.$ (53)

Then, the elements $x_{k}^{+}$ in Eq. (10) become

$x_{1}^{+}=\hat{x}+2\textrm{i}y,\quad
x_{k}^{+}=\frac{1}{k!}\hat{x}+\frac{2^{k}}{k!}\textrm{i}y+T_{k},$ (54)

where $T_{k}\equiv 12(1-3^{k-1})t/k!$. In particular,

$T_{2}=-12t,\quad T_{3}=-16t.$ (55)

In this moving $x$-frame, when $|t|$ is large and
$\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$, we have the leading-order asymptotics
for
$S_{k}\left(\hat{\textbf{\emph{x}}}^{+}+\nu\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}\right)$
as

$S_{k}\left(\hat{\textbf{\emph{x}}}^{+}+\nu\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}\right)\sim
S_{k}(\textbf{v}),\quad\quad|t|\gg 1,$ (56)

where

$\textbf{v}=\left(x_{1}^{+},\hskip 1.13791pt0,\hskip 1.13791ptT_{3},\hskip
1.13791pt0,\hskip 1.13791pt0,\hskip 1.13791pt0,\cdots\right).$ (57)

By comparing the definition of Schur polynomials $S_{k}(\textbf{v})$ to the
definition of $p_{k}(z)$ polynomials in Eq. (17), we see that

$S_{k}(\textbf{v})=\left(-3T_{3}/4\right)^{k/3}p_{k}(z),$ (58)

where

$z=(-3T_{3}/4)^{-1/3}x_{1}^{+}=(-3T_{3}/4)^{-1/3}\left(\hat{x}+2\textrm{i}y\right).$
(59)

Using these formulae and the Laplace expansion of the $3N\times 3N$
determinant (49)

$\displaystyle\sigma=\sum_{0\leq\nu_{1}<\nu_{2}<\cdots<\nu_{N}\leq
2N-1}\det_{1\leq i,j\leq
N}\left[\frac{1}{2^{\nu_{j}}}S_{2i-1-\nu_{j}}(\hat{\textbf{\emph{x}}}^{+}+\nu_{j}\textbf{\emph{s}}+\hat{\textbf{\emph{a}}})\right]\times\det_{1\leq
i,j\leq
N}\left[\frac{1}{2^{\nu_{j}}}S_{2i-1-\nu_{j}}[(\hat{\textbf{\emph{x}}}^{+})^{*}+\nu_{j}\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}^{*}]\right],$
(60)

together with the fact that the highest order term of $|t|$ in this $\sigma$
comes from the index choices of $\nu_{j}=j-1$, we can readily show that the
highest $t$-power term of $\sigma$ is

$\sigma\sim|\alpha_{0}|^{2}\hskip
1.42271pt\left|3T_{3}/4\right|^{\frac{N(N+1)}{3}}\left|Q_{N}(z)\right|^{2},\quad\quad|t|\gg
1,$ (61)

where $\alpha_{0}=2^{-N(N-1)/2}c_{N}^{-1}$. Inserting this leading-order term
of $\sigma$ into Eq. (7), we see that the solution $u_{\Lambda}(x,y,t)$
approaches zero when $|t|\to\infty$, except at or near $(\hat{x},y)$ locations
$\left(x_{0},y_{0}\right)$, i.e., at or near $(x,y)$ locations
$\left(12t+x_{0},y_{0}\right)$, where

$z_{0}=(-3T_{3}/4)^{-1/3}\left(x_{0}+2\textrm{i}y_{0}\right)$ (62)

is a root of the polynomial $Q_{N}(z)$. Solving this equation, we get the
$(x_{0},y_{0})$ locations given by Eq. (36) in Theorem 2. Due to our
requirement of $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$, $z_{0}$ in the above
equation should be nonzero.

In order to derive the solution behavior near this
$(x,y)=\left(12t+x_{0},y_{0}\right)$ location, we need to perform a more
refined asymptotic analysis and calculate the next-order terms in $t$, since
the leading-order term in Eq. (61) vanishes at this point. Recalling
$s_{1}=0$, this refined analysis is very similar to that we did for rogue
waves in the NLS equation YangYang21a . For $z_{0}\neq 0$ in the
$(x_{0},y_{0})$ formula (36), i.e., if the $\left(12t+x_{0},y_{0}\right)$
location is $O(|t|^{1/3})$ away from the wave center $(12t,0)$, then in the
$O(1)$ neighborhood of $\left(12t+x_{0},y_{0}\right)$, i.e., when
$(x-12t-x_{0})^{2}+(y-y_{0})^{2}=O(1)$, we have an asymptotics more refined
than (56), which is

$S_{k}\left(\hat{\textbf{\emph{x}}}^{+}+\nu\textbf{\emph{s}}+\hat{\textbf{\emph{a}}}\right)=S_{k}(\textbf{v})\left[1+O(|t|^{-2/3})\right].$
(63)

This $O(|t|^{-2/3})$ relative error is due to our omission of
$\hat{x}/6+4\textrm{i}y/3$ relative to $T_{3}$ in $x_{3}^{+}$, and omission of
$x_{5}^{+}$ relative to $x_{3}^{+}$. Using this refined asymptotics and
repeating the same steps as in YangYang21a , we find that

$\displaystyle\sigma(x,y,t)=\left|\alpha_{0}\right|^{2}\hskip
1.70709pt\left|Q_{N}^{\prime}(z_{0})\right|^{2}|3T_{3}/4|^{\frac{N(N+1)-2}{3}}\left[\left(x-12t-x_{0}\right)^{2}+4(y-y_{0})^{2}+\frac{1}{4}\right]\left[1+O\left(|t|^{-1/3}\right)\right],$
(64)

where $\alpha_{0}$ is given below Eq. (61). For Yablonskii-Vorob’ev
polynomials $Q_{N}(z)$, all roots are simple. Thus, $Q_{N}^{\prime}(z_{0})\neq
0$.

In the $O(1)$ neighborhood of the wave center $\left(12t,0\right)$, where
$(x-12t)^{2}+y^{2}=O(1)$, we need to perform a separate asymptotic analysis,
because the earlier $S_{k}$ asymptotics (56) and (63) do not hold in this
region. In this case, due to Eq. (54), when we lump $T_{2k+1}$ and $a_{2k+1}$
together in Eq. (51) and recall $T_{2k+1}$ is proportional to $t$, the large-
time analysis of the present $\sigma$ determinant (49) is almost identical to
that in Appendix C of Ref. YangYang21a  for the analysis of NLS rogue patterns
when its internal parameters $(a_{3},a_{5},\cdots)$ are all large and of the
same order. Repeating that analysis, we find that if zero is a root of the
Yablonskii-Vorob’ev polynomial $Q_{N}(z)$, i.e., $N\equiv 1\hskip
2.84544pt\mbox{mod}\hskip 2.84544pt3$, then

$\displaystyle\sigma(x,y,t)=\beta_{0}\hskip
1.42271pt|t|^{\frac{N(N+1)-2}{3}}\left[\left(x-12t\right)^{2}+4y^{2}+\frac{1}{4}\right]\left[1+O\left(|t|^{-1}\right)\right],$
(65)

where $\beta_{0}$ is a $N$-dependent positive constant. If zero is not a root
of $Q_{N}(z)$, then $\sigma(x,y,t)\sim\beta_{0}\hskip
1.42271pt|t|^{\frac{(N+2)(N-1)}{3}}$.

Substituting the above two $\sigma$ asymptotics (64)-(65) into the solution
expression (7) and performing a little simplification, we then get the
asymptotics (37). Theorem 2 is then proved.

###  V.2 Proof of Theorem 3

In this case, $\Lambda\neq(1,3,5,\dots,2N-1)$. We first rewrite the
determinant $\sigma$ in (8) as a larger $(N+n_{N}+1)\times(N+n_{N}+1)$
determinant

$\sigma=\left|\begin{array}[]{cc}\textbf{O}_{N\times
N}&\Phi_{N\times(n_{N}+1)}\\\ -\Psi_{(n_{N}+1)\times
N}&\textbf{I}_{(n_{N}+1)\times(n_{N}+1)}\end{array}\right|,$ (66)

where

$\Phi_{i,j}=2^{-(j-1)}S_{n_{i}+1-j}\left[\textbf{\emph{x}}^{+}+(j-1)\textbf{\emph{s}}+\textbf{\emph{a}}\right],\quad\Psi_{i,j}=2^{-(i-1)}S_{n_{j}+1-i}\left[(\textbf{\emph{x}}^{+})^{*}+(i-1)\textbf{\emph{s}}+\textbf{\emph{a}}^{*}\right],$
(67)

and vectors $\textbf{\emph{x}}^{+}$ and s are given in Eqs. (10) and (34).
Unlike the previous case, we cannot eliminate $x_{2}^{+}$ from this solution
now. Our large-time asymptotics of this determinant proceeds as follows.

####  V.2.1 Proof for the outer region

First, we prove the asymptotics (40)-(41) for the outer region. In this
region, $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/2})$. Thus, we have the leading-
order asymptotics for
$S_{k}\left(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}}\right)$
as

$S_{k}\left(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}}\right)\sim
S_{k}(\textbf{w}),\quad\quad|t|\gg 1,$ (68)

where

$\textbf{w}=\left(x_{1}^{+},\hskip 1.13791ptT_{2},\hskip 1.13791pt0,\hskip
1.13791pt0,\hskip 1.13791pt0,\cdots\right),$ (69)

and $T_{2}$ is as given in Eq. (55). By comparing the definition of Schur
polynomials $S_{k}(\textbf{w})$ to the definition of $q_{k}(z)$ polynomials in
Eq. (23), we see that

$S_{k}(\textbf{w})=T_{2}^{k/2}q_{k}(z),$ (70)

where

$z=T_{2}^{-1/2}x_{1}^{+}=T_{2}^{-1/2}\left(\hat{x}+2\textrm{i}y\right).$ (71)

Using these formulae and the Laplace expansion of the determinant (66) for
$\sigma$, we can readily show that the highest $t$-power term of $\sigma$ is

$\sigma\sim|\mu_{0}|^{2}\hskip
1.42271pt|T_{2}|^{\rho}\left|W_{\Lambda}(z)\right|^{2},\quad\quad|t|\gg 1,$
(72)

where $\rho$ is given in Eq. (13), and $\mu_{0}=2^{-N(N-1)/2}$. Inserting this
leading-order term of $\sigma$ into Eq. (7), we see that the solution
$u_{\Lambda}(x,y,t)$ approaches zero when $|t|\to\infty$, except at or near
$(\hat{x},y)$ locations $\left(\hat{x}_{0},\hat{y}_{0}\right)$, i.e., at or
near $(x,y)$ locations $\left(12t+\hat{x}_{0},\hat{y}_{0}\right)$, where

$z_{0}=T_{2}^{-1/2}\left(\hat{x}_{0}+2\textrm{i}\hat{y}_{0}\right)$ (73)

is a root of the Wronskian-Hermit polynomial $W_{\Lambda}(z)$. Solving this
equation, we get

$\hat{x}_{0}=\Re\left[z_{0}T_{2}^{1/2}\right],\quad\hat{y}_{0}=\frac{\Im\left[z_{0}T_{2}^{1/2}\right]}{2},$
(74)

which are the leading-order terms of $(x_{0},y_{0})$ in Eq. (40) of Theorem 3.
Due to our requirement of $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/2})$, $z_{0}$ in
the above equation should be nonzero.

To derive the solution behavior near this
$(x,y)=\left(12t+\hat{x}_{0},\hat{y}_{0}\right)$ location, we perform a more
refined asymptotic analysis. Our starting point is a more accurate asymptotics
for
$S_{k}\left(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}}\right)$,

$S_{k}(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}})=S_{k}(\hat{\textbf{w}})\left[1+O\left(|t|^{-1}\right)\right],\quad\quad|t|\gg
1,$ (75)

where

$\hat{\textbf{w}}=\left(x_{1}^{+},\hskip 1.13791ptx_{2}^{+},\hskip
1.13791ptT_{3},\hskip 1.13791pt0,\hskip 1.13791pt0,\hskip
1.13791pt0,\cdots\right)=\textbf{w}+\left(0,\hat{x}_{2}^{+},T_{3},\hskip
1.13791pt0,\hskip 1.13791pt0,\hskip 1.13791pt0,\cdots\right),$ (76)

w is given in (69), and

$\hat{x}_{2}^{+}\equiv\frac{1}{2}\hat{x}+2\textrm{i}y.$ (77)

The asymptotics (75) holds since $a_{1}=s_{1}=0$. From the definition (6) of
Schur polynomials and the above equation, we can relate
$S_{k}(\hat{\textbf{w}})$ and $S_{k}(\textbf{w})$ as

$S_{k}(\hat{\textbf{w}})=\sum_{j=0}^{k}b_{j}S_{k-j}(\textbf{w}),$ (78)

where $b_{j}$ are the coefficients in the expansion

$e^{\hat{x}_{2}^{+}\epsilon^{2}+T_{3}\hskip
0.56917pt\epsilon^{3}}=\sum_{j=0}^{\infty}b_{j}\epsilon^{j}.$ (79)

Notice that $b_{0}=1$, $b_{1}=0$, $b_{2}=\hat{x}_{2}^{+}$, and $b_{3}=T_{3}$.
In addition, $(\hat{x},y)=O(|t|^{1/2})$ from Eq. (74), and
$S_{k}(\textbf{w})=O(|t|^{k/2})$ in view of Eq. (70). Utilizing these
relations, we find that

$S_{k}(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}})=\left[S_{k}(\textbf{w})+\hat{x}_{2}^{+}S_{k-2}(\textbf{w})+T_{3}\hskip
0.56917ptS_{k-3}(\textbf{w})\right]\left[1+O\left(|t|^{-1}\right)\right],\quad\quad|t|\gg
1.$ (80)

With this formula (80), we can now determine the asymptotic expression of
$\sigma$ in Eq. (66) in the neighborhood of
$(x,y)=\left(12t+\hat{x}_{0},\hat{y}_{0}\right)$ at large $t$. The Laplace
expansion of this determinant is very similar to Eq. (60) of the previous
subsection. Using this Laplace expansion and similar techniques as in Refs.
YangYang21a ; YangYang21b , we can readily find that

$\displaystyle\sigma(x,y,t)=\left|\mu_{0}\right|^{2}\hskip
1.70709pt\left|W_{\Lambda}^{\prime}(z_{0})\right|^{2}|T_{2}|^{\rho-1}\left(\left|(x-12t-\hat{x}_{0})+2\textrm{i}(y-\hat{y}_{0})+\Delta\right|^{2}+\frac{1}{4}\right)\left[1+O\left(|t|^{-1/2}\right)\right],$
(81)

where $\mu_{0}$ is given below Eq. (72), and $\Delta=\Delta(\Lambda,z_{0})$ is
an $O(1)$ complex constant given by

$\displaystyle\Delta=\frac{1}{W_{\Lambda}^{\prime}(z_{0})}\left\\{\lambda\sum_{j=1}^{N}\det_{1\leq
i\leq N}\left[q_{n_{i}},\cdots,q_{n_{i}-(j-2)},\hskip
1.42271ptq_{n_{i}-(j-1)-2},\hskip
1.42271ptq_{n_{i}-j},\cdots,q_{n_{i}-(N-1)}\right]_{z=z_{0}}\right.$
$\displaystyle\hskip 51.21504pt\left.+\frac{4}{3}\sum_{j=1}^{N}\det_{1\leq
i\leq N}\left[q_{n_{i}},\cdots,q_{n_{i}-(j-2)},\hskip
1.42271ptq_{n_{i}-(j-1)-3},\hskip
1.42271ptq_{n_{i}-j},\cdots,q_{n_{i}-(N-1)}\right]_{z=z_{0}}\right\\},$ (82)

and

$\lambda=\left\\{\begin{array}[]{ll}\frac{1}{2}\Re(z_{0})+\textrm{i}\Im(z_{0}),&\mbox{when}\hskip
2.84544ptt<0,\\\ \Re(z_{0})+\frac{1}{2}\textrm{i}\Im(z_{0}),&\mbox{when}\hskip
2.84544ptt>0.\end{array}\right.$ (83)

The former determinant in Eq. (82) is the Wronskian-Hermit determinant in Eq.
(25) but with the $j$-th column $\\{q_{n_{i}-(j-1)}\\}$ replaced by
$\\{q_{n_{i}-(j-1)-2}\\}$, i.e., reducing the subindex value of this column by
two, while the latter determinant in (82) is the Wronskian-Hermit determinant
(25) with the $j$-th column replaced by $\\{q_{n_{i}-(j-1)-3}\\}$, i.e.,
reducing its subindex value by three.

The complex constant $\Delta$ in Eq. (81) can be absorbed into
$(\hat{x}_{0},\hat{y}_{0})$. After this absorption and rearranging terms, Eq.
(81) becomes

$\displaystyle\sigma(x,y,t)=\left|\mu_{0}\right|^{2}\hskip
1.70709pt\left|W_{\Lambda}^{\prime}(z_{0})\right|^{2}|T_{2}|^{\rho-1}\left[\left(x-12t-x_{0}\right)^{2}+4(y-y_{0})^{2}+\frac{1}{4}\right]\left[1+O\left(|t|^{-1/2}\right)\right],$
(84)

where

$x_{0}=\Re\left[z_{0}(-12t)^{1/2}\right]-\Re(\Delta),\quad
y_{0}=\frac{\Im\left[z_{0}(-12t)^{1/2}\right]}{2}-\frac{\Im(\Delta)}{2}.$ (85)

These $(x_{0},y_{0})$ formulae contain the explicit $O(1)$ corrections to the
leading $O(|t|^{1/2})$ terms, and are thus more complete than Eq. (40) in
Theorem 3. Substituting the above $\sigma$ asymptotics (84) into Eq. (7), the
asymptotics (40)-(41) for the outer region of Theorem 3 are then proved.

####  V.2.2 Proof for the inner region

In the inner region, where $\sqrt{\hat{x}^{2}+y^{2}}\leq O(|t|^{1/3})$, a
separate asymptotic analysis is needed, because the previous $S_{k}$
asymptotics (68) and (75) do not hold. In this inner region, our analysis
needs to split into two cases, depending on whether $\sqrt{\hat{x}^{2}+y^{2}}$
is $O(|t|^{1/3})$ or $O(1)$.

When $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$, it is easy to see from the
Laplace expansion of the $\sigma$ determinant (66) that, at large $|t|$, the
highest $t$-power term of $\sigma$ comes from the index choices of
$\nu_{j}=j-1$, i.e.,

$\displaystyle\sigma\sim\left|\det_{1\leq i,j\leq
N}\left[\frac{1}{2^{j-1}}S_{n_{i}+1-j}(\textbf{\emph{x}}^{+}+(j-1)\textbf{\emph{s}}+\textbf{\emph{a}})\right]\right|^{2},\qquad|t|\gg
1.$ (86)

For the determinant in the above equation, we reorganize its rows by grouping
odd-$n_{i}$ rows together (in ascending order of $n_{i}$), followed by
even-$n_{i}$ rows (also in ascending order of $n_{i}$). We also rewrite
$S_{k}(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}})$ as

$S_{k}(\textbf{\emph{x}}^{+}+\nu\textbf{\emph{s}}+\textbf{\emph{a}})=\sum_{j=0}^{\left[k/2\right]}\frac{T_{2}^{j}}{j!}S_{k-2j}(\textbf{h}+\nu\textbf{\emph{s}}+\textbf{\emph{a}}),$
(87)

where

$\textbf{h}\equiv\textbf{\emph{x}}^{+}-(0,T_{2},0,0,\cdots)=(x_{1}^{+},\hat{x}_{2}^{+},x_{3}^{+},x_{4}^{+},\cdots),$
(88)

and $[a]$ represents the largest integer less than or equal to $a$. In
addition, we notice that when $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$,

$S_{k}(\textbf{h}+\nu\textbf{\emph{s}}+\textbf{\emph{a}})\sim
S_{k}(\hat{\textbf{h}}),$ (89)

where

$\hat{\textbf{h}}=(x_{1}^{+},0,T_{3},0,0,\cdots),$ (90)

and

$S_{k}(\hat{\textbf{h}})=\left(-3T_{3}/4\right)^{k/3}p_{k}(z),$ (91)

where $z$ is as given in Eq. (59). Inserting these formulae into the above
reorganized determinant and performing row operations to eliminate certain
high powers of $T_{2}$ in lower rows of the odd-$n_{i}$-index group as well as
the even-$n_{i}$-index group, we find that the highest $t$-power term of
$\sigma$ from Eq. (86) is

$\sigma\sim\gamma_{0}\hskip
1.42271pt|T_{2}|^{N_{W}}|3T_{3}/4|^{\frac{\hat{d}(\hat{d}+1)}{3}}\left|H(z)\right|^{2},$
(92)

where $\gamma_{0}$ is a certain positive constant, $N_{W}$ is given in Eqs.
(27)-(28), $\hat{d}$ is defined in Eq. (43), $H(z)$ is the determinant

$H(z)=\det\left[\begin{array}[]{cccc}p_{1}(z)&p_{0}(z)&p_{-1}(z)&\cdots\\\
p_{3}(z)&p_{2}(z)&p_{1}(z)&\cdots\\\ \vdots&\vdots&\vdots&\vdots\\\
p_{2k_{odd}-1}(z)&p_{2k_{odd}-2}(z)&p_{2k_{odd}-3}(z)&\cdots\\\
p_{0}(z)&p_{-1}(z)&p_{-2}(z)&\cdots\\\ p_{2}(z)&p_{1}(z)&p_{0}(z)&\cdots\\\
\vdots&\vdots&\vdots&\vdots\\\
p_{2k_{even}-2}(z)&p_{2k_{even}-3}(z)&p_{2k_{even}-4}(z)&\cdots\end{array}\right],$
(93)

and $k_{odd}$, $k_{even}$ are the numbers of odd and even elements in the
index vector $(n_{1},n_{2},\dots,n_{N})$ respectively. Clearly, this $H$
determinant can be reduced to

$H(z)=Q_{\hat{d}}(z).$ (94)

Thus, when $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$,

$\sigma\sim\gamma_{0}\hskip
1.42271pt|T_{2}|^{N_{W}}|3T_{3}/4|^{\frac{\hat{d}(\hat{d}+1)}{3}}\left|Q_{\hat{d}}(z)\right|^{2}.$
(95)

In view of Eq. (7), this asymptotics shows that the solution
$u_{\Lambda}(x,y,t)$ is asymptotically zero in this region, except when
$\hat{d}>0$ (i.e., zero is a root of the Wronskian-Hermit polynomial
$W_{\Lambda}(z)$), and when $(x,y)$ is at or near the location
$\left(12t+\hat{x}_{0},\hat{y}_{0}\right)$, where

$z_{0}=(-3T_{3}/4)^{-1/3}\left(\hat{x}_{0}+2\textrm{i}\hat{y}_{0}\right)$ (96)

is a root of the Yablonskii–Vorob’ev polynomial $Q_{\hat{d}}(z)$. Solving this
equation, we get $(\hat{x}_{0},\hat{y}_{0})$ values that are the leading-order
terms in Eq. (42) of Theorem 3. Since $\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$,
root $z_{0}$ in the above equation should be nonzero.

We can further show that, near this
$(x,y)=\left(12t+\hat{x}_{0},\hat{y}_{0}\right)$ location lies a fundamental
lump. This calculation is similar to that we did in the proof of Theorem 2 and
the earlier part of this proof of Theorem 3. Specifically, we can show that in
the $O(1)$ neighborhood of this location,

$\displaystyle\sigma(x,y,t)=\gamma_{0}\hskip
1.42271pt|T_{2}|^{N_{W}}|3T_{3}/4|^{\frac{\hat{d}(\hat{d}+1)-2}{3}}\left|Q_{\hat{d}}^{\prime}(z_{0})\right|^{2}\left(\left|(x-12t-\hat{x}_{0})+2\textrm{i}(y-\hat{y}_{0})+\hat{\Delta}\right|^{2}+\frac{1}{4}\right)\left[1+O\left(|t|^{-1/3}\right)\right],$
(97)

where $\hat{\Delta}=\hat{\Delta}(\Lambda,z_{0})$ is an $O(1)$ constant. This
$\hat{\Delta}$ is the analog of a similar quantity $\Delta$ which we derived
in Eq. (82) for a fundamental lump in the outer region. It is easy to see that
the above $\sigma(x,y,t)$ gives a fundamental lump, whose position is at
$(x,y)=\left(12t+x_{0},y_{0}\right)$, where

$x_{0}=\Re(z_{0})\hskip 1.42271pt(12t)^{1/3}-\Re(\hat{\Delta}),\quad
y_{0}=\frac{\Im(z_{0})}{2}(12t)^{1/3}-\frac{1}{2}\Im(\hat{\Delta}),$ (98)

which matches (42) in Theorem 3. In addition, the error of this prediction is
$O\left(|t|^{-1/3}\right)$.

In the center region where $\hat{x}^{2}+y^{2}=O(1)$, we can use the technique
of Appendix C in Ref. YangYang21a  to show that at large time, if zero is a
root of the Yablonskii–Vorob’ev polynomial $Q_{\hat{d}}(z)$, i.e., if
$\hat{d}\equiv 1\hskip 2.84544pt\mbox{mod}\hskip 2.84544pt3$, then
$u_{\Lambda}(x,y,t)$ would approach a fundamental lump located in the $O(1)$
neighborhood of the wave center $(\hat{x},y)=(0,0)$. If zero is not a root of
$Q_{\hat{d}}(z)$, then $u_{\Lambda}(x,y,t)$ would approach zero in this center
region as $|t|\to\infty$. Details will be omitted for brevity. It may be more
illuminating for us to point out that, the leading-order term of the previous
asymptotic formula (97), which was derived for the region of
$\sqrt{\hat{x}^{2}+y^{2}}=O(|t|^{1/3})$ and nonzero roots $z_{0}$ of
$Q_{\hat{d}}(z)$, turns out to be valid for the $\hat{x}^{2}+y^{2}=O(1)$
region and the zero root $z_{0}$ of $Q_{\hat{d}}(z)$ as well [except for the
relative error term, which is now $O(|t|^{-1})$ rather than $O(|t|^{-1/3})$].
In other words, if zero is a root of $Q_{\hat{d}}(z)$, then setting $z_{0}=0$
in the leading term of (97), we would get the correct asymptotic fundamental
lump in the $\hat{x}^{2}+y^{2}=O(1)$ region. In particular, the location of
this fundamental lump would be at $(\hat{x},y)=(x_{0},y_{0})$, i.e.,
$(x,y)=(12t+x_{0},y_{0})$, where $x_{0}$ and $y_{0}$ are given by (98) with
$z_{0}=0$ and $\hat{\Delta}=\hat{\Delta}|_{z_{0}=0}$. This completes the proof
of Theorem 3.

##  VI Summary and Discussion

In this article, we have analytically studied pattern formation in higher-
order lumps of the KP-I equation at large time. For a broad class of these
higher-order lumps, we have shown that two types of solution patterns appear
at large time. The first type of patterns comprise fundamental lumps arranged
in triangular shapes, which are described analytically by root structures of
Yablonskii–Vorob’ev polynomials. As time evolves from large negative to large
positive, this triangular pattern reverses its $x$-direction. The second type
of solution patterns comprise fundamental lumps arranged in non-triangular
shapes in the outer region, which are described analytically by nonzero-root
structures of Wronskian–Hermit polynomials, together with possible fundamental
lumps arranged in triangular shapes in the inner region, which are described
analytically by root structures of Yablonskii–Vorob’ev polynomials. When time
evolves from large negative to large positive, the non-triangular pattern in
the outer region switches its $x$ and $y$ directions, while the triangular
pattern in the inner region reverses its $x$-direction. We have also compared
these predicted patterns with true solutions, and excellent agreement is
observed.

In this pattern analysis of higher-order lumps, we have set the spectral
parameter $p=1$ without any loss of generality (see Remark 4 of Sec. 2).
Because of this, lump patterns we have predicted at large time are all
$y$-symmetric (see Figs. 3, 5 and 7), since root structures of
Yablonskii–Vorob’ev and Wronskian-Hermit polynomials are symmetric with
respect to the real-$z$ axis. However, under the Galilean transformation (4),
these $y$-symmetric lump patterns can become skewed and $y$-asymmetric, and
these $y$-asymmetric patterns correspond to complex spectral parameters $p$.
Thus, $y$-asymmetric lump patterns also exist in the KP-I equation, and such
patterns can be obtained from the $y$-symmetric ones through the Galilean
transformation.

Are there other patterns of higher-order lumps at large time? The answer is
yes. Notice that in this article, we have assumed internal-parameter vectors
$\textbf{\emph{a}}_{i}$ of higher-order lumps to be equal to each other [see
Eq. (35)]. If these parameter vectors are allowed to differ from each other,
then the analytical results at large time will become different. This problem
will not be pursued in this paper, and will be left for future studies.

In a very recent preprint Ling2021 , the authors also derived higher-order
lumps in the KP-I equation and studied their large-time patterns through
Darboux transformation, and showed that their large-time patterns are
described analytically by root structures of Yablonskii–Vorob’ev polynomials.
Obviously, the higher-order lump solutions they derived are a very special
class of solutions which correspond to the index vector
$\Lambda=(1,3,5,\dots,2N-1)$ and under $\textbf{\emph{a}}_{i}$ parameter
constraints (35) in our general solutions of Theorem 1, and their large-time
pattern results are largely equivalent to our Theorem 2. However, their error
estimate of $O(|t|^{-2/3})$ for fundamental-lump predictions far away from the
wave center is different from our $O(|t|^{-1/3})$ in Theorem 2, and we have
verified numerically that their error estimate is incorrect. More importantly,
those authors have not considered the more general higher-order lumps
corresponding to the index vector $\Lambda\neq(1,3,5,\dots,2N-1)$ in our
Theorem 1, nor their large-time solution patterns. These latter patterns are
the contents of our Theorem 3 (see also our Figs. 5-8).

## Acknowledgment

This material is based upon work supported by the National Science Foundation
under award number DMS-1910282, and the Air Force Office of Scientific
Research under award number FA9550-18-1-0098.

Appendix

In this appendix, we briefly derive the bilinear higher-order lump solutions
presented in Theorem 1.

Under the variable transformation $u=2(\log\tau)_{xx}$ and notations of
$x_{1}=x,x_{2}=\textrm{i}y$ and $x_{3}=-4t$, the KP-I equation (3) is
converted to the bilinear equation

$(D_{x_{1}}^{4}-4D_{x_{1}}D_{x_{3}}+3D_{x_{2}}^{2})\,\tau\cdot\tau=0,$ (99)

where $D$ is Hirota bilinear differential operator. It is well known that if
$m_{ij}$, $\phi_{i}$ and $\psi_{j}$ are functions of $(x_{1},x_{2},x_{3})$ and
satisfy the following differential equations

$\displaystyle\partial_{x_{1}}m_{ij}=\phi_{i}\psi_{j},$ (100)
$\displaystyle\partial_{x_{n}}\phi_{i}=\partial_{x_{1}}^{n}\phi_{i},\hskip
42.67912ptn=2,3,$ (101)
$\displaystyle\partial_{x_{n}}\psi_{j}=(-1)^{n-1}\partial_{x_{1}}^{n}\psi_{j},\quad
n=2,3,$ (102)

then the $\tau$ function

$\tau=\det_{1\leq i,j\leq N}\left(m_{ij}\right)$ (103)

would satisfy the above bilinear equation Hirota_book . To derive higher-order
lump solutions, we define $m_{ij}$, $\phi_{i}$ and $\psi_{j}$ as

$m_{ij}=\mathcal{A}_{i}\mathcal{B}_{j}\frac{1}{p+q}e^{\xi_{i}+\eta_{j}},\quad\phi_{i}=\mathcal{A}_{i}e^{\xi_{i}},\quad\psi_{j}=\mathcal{B}_{j}e^{\eta_{j}},$
(104)

where

$\displaystyle\mathcal{A}_{i}=\frac{1}{n_{i}!}(p\partial_{p})^{n_{i}},\quad\mathcal{B}_{j}=\frac{1}{n_{j}!}(q\partial_{q})^{n_{j}},$
(105)
$\xi_{i}=px_{1}+p^{2}x_{2}+p^{3}x_{3}+\xi_{i,0}(p),\quad\eta_{j}=qx_{1}-q^{2}x_{2}+q^{3}x_{3}+\eta_{j,0}(q),$
(106)

$(n_{1},n_{2},\cdots,n_{N})$ is a vector of arbitrary positive integers, $p,q$
are arbitrary complex constants, and $\xi_{i,0}(p)$, $\eta_{j,0}(q)$ are
arbitrary complex functions of $p$ and $q$. It is easy to see that these
$m_{ij}$, $\phi_{i}$ and $\psi_{j}$ functions satisfy the differential
equations (100)-(102). Thus, the above $\tau$ function would satisfy the
bilinear equation (99). To guarantee that this $\tau$ function is real-valued,
we impose the parameter constraints

$q=p^{*},\quad\eta_{j,0}(q)=[\xi_{j,0}(p)]^{*}.$ (107)

Under these constraints, $\eta_{j}=\xi_{j}^{*}$,
$m_{n_{i},n_{j}}^{*}=m_{n_{j},n_{i}}$, and thus $\tau$ in (103) is real.
Following the technique of Ref. OhtaJY2012 , we can further show this $\tau$
is positive. Thus, the resulting function $u=2(\log\tau)_{xx}$ is a real-
valued solution to the KP-I equation (3).

Next, we need to simplify the matrix elements of this $\tau$ determinant and
derive their more explicit algebraic expressions. This simplification is very
similar to that we performed in OhtaJY2012 ; YangYang3wave . By expanding
$\xi_{i,0}(p)$ into a certain series containing complex parameters
$\textbf{\emph{a}}_{i}=\left(a_{i,1},a_{i,2},\cdots\right)$ and repeating the
calculations of OhtaJY2012 ; YangYang3wave , we can show that the matrix
element $m_{ij}$ in (104) can be reduced to the expression given in Eq. (9) of
Theorem 1.

We would like to make a comment here regarding the choice of differential
operators in Eq. (105). Obviously, we can also choose more general forms of
these differential operators, such as

$\displaystyle\mathcal{A}_{i}=\frac{1}{n_{i}!}\left[f(p)\partial_{p}\right]^{n_{i}},\quad\mathcal{B}_{j}=\frac{1}{n_{j}!}\left[f(q)\partial_{q}\right]^{n_{j}},$
(108)

where $f(p)$ is an arbitrary function, and the resulting $\tau$ function (103)
would still satisfy the bilinear equation (99). However, such additional
freedoms in the differential operators will not produce new higher-order lump
solutions. To see why, we can rewrite this $\mathcal{A}_{i}$ as

$\mathcal{A}_{i}=\frac{1}{n_{i}!}\left[\frac{f(p)}{p}p\partial_{p}\right]^{n_{i}}=\sum_{k=0}^{n_{i}}c_{i,k}\frac{1}{(n_{i}-k)!}(p\partial_{p})^{n_{i}-k},$
(109)

where $c_{i,k}$ are $p$-dependent complex constants. Similar treatments can be
made on $\mathcal{B}_{j}$. These differential operators in summation form are
similar to those taken in Ref. OhtaJY2012 . We can directly show that the
$m_{ij}$ matrix element with these differential operators of summation form
can be converted to one with these differential operators as a single term in
(105), after parameters $\textbf{\emph{a}}_{i}$ in the series expansion of
$\xi_{j,0}(p)$ are redefined properly. Thus, no new solutions are produced.

## References

## References

  * (1) B.B. Kadomtsev and V.I. Petviashvili, “On the stability of solitary waves in weakly dispersive media”, Sov. Phys. Dokl. 15, 539-541 (1970).
  * (2) M.J. Ablowitz and H. Segur, “On the evolution of packets of water waves”, J. Fluid Mech. 92, 691-715 (1979).
  * (3) D.E. Pelinovsky, Yu.A. Stepanyants and Yu.A. Kivshar, “Self-focusing of plane dark solitons in nonlinear defocusing media”, Phys. Rev. E 51, 5016-5026 (1995).
  * (4) S. Tsuchiya, F. Dalfovo and L.P. Pitaevskii, “Solitons in two-dimensional Bose-Einstein condensates”, Phys. Rev. A 77, 045601 (2008).
  * (5) S. Novikov, S.V. Manakov, L.P. Pitaevskii and V.E. Zakharov, _Theory of solitons: the inverse scattering method_ (Plenum, 1984).
  * (6) M.J. Ablowitz and P.A. Clarkson, _Solitons, Nonlinear Evolution Equations and Inverse Scattering_ (Cambridge University Press, 1991).
  * (7) V.I. Petviashvili, “Equation of an extraordinary soliton”, Plasma Physics 2, 469-472 (1976).
  * (8) S.V. Manakov, V.E. Zakharov, L.A. Bordag, A.R. Its and V.B. Matveev, “Two-dimensional solitons of the Kadomtsev–Petviashvili equation and their interaction”. Phys. Lett. A 63, 205-206 (1977).
  * (9) J. Satsuma and M.J. Ablowitz, “Two-dimensional lumps in nonlinear dispersive systems”, J. Math. Phys. 20, 1496 (1979).
  * (10) K. A. Gorshkov, D. E. Pelinovsky, and Yu. A. Stepanyants, “Normal and anomalous scattering, formation and decay of bound states of two-dimensional solitons described by the Kadomtsev–Petviashvili equation,” JETP 77, 237-245 (1993).
  * (11) M.J. Ablowitz and J. Villarroel, “Solutions to the time dependent Schrödinger and the Kadomtsev-Petviashvili equations,” Phys. Rev. Lett. 78, 570-573 (1997).
  * (12) M.J. Ablowitz, S. Chakravarty, A.D. Trubatch and J. Villarroel, “A novel class of solutions of the non-stationary Schrödinger and the Kadomtsev–Petviashvili I equations”, Phys. Lett. A 267, 132-146 (2000).
  * (13) D.E. Pelinovsky and Yu. A. Stepanyants, “New multisoliton solutions of the Kadomtsev–Petviashvili equation”, JETP Lett. 57, 24-28 (1993).
  * (14) D. Pelinovsky, “Rational solutions of the KP hierarchy and the dynamics of their poles. II. Construction of the degenerate polynomial solutions”, J. Math. Phys. 39, 5377-5395 (1998).
  * (15) P. Dubard, P. Gaillard, C. Klein and V.B. Matveev, “On multi-rogue wave solutions of the NLS equation and positon solutions of the KdV equation”, Eur. Phys. J. Spec. Top. 185, 247-258 (2010).
  * (16) P. Dubard and V.B. Matveev, “Multi-rogue waves solutions: from the NLS to the KP-I equation”, Nonlinearity 26, R93-R125 (2013).
  * (17) S. Chen, P. Grelu, D. Mihalache and F. Baronio, “Families of rational solution solutions of the Kadomtsev–Petviashvili I equation”, Rom. Rep. Phys. 68, 1407-1424 (2016).
  * (18) P.A. Clarkson and E. Dowie, “Rational solutions of the Boussinesq equation and applications to rogue waves”, Trans. Math. Appl. 1, 1-26 (2017).
  * (19) P. Gaillard, “Multiparametric families of solutions of the Kadomtsev–Petviashvili-I equation, the structure of their rational representations, and multi-rogue waves”, Theor. Math. Phys. 196, 1174-1199 (2018).
  * (20) J.H. Chang, “Asymptotic analysis of multilump solutions of the Kadomtsev–Petviashvili-I equation”, Theor. Math. Phys. 195, 676-689 (2018).
  * (21) W.X. Ma, “Lump solutions to the Kadomtsev–Petviashvili equation”, Phys. Lett. A 379, 1975-1978 (2015).
  * (22) C. Lester, A. Gelash, D. Zakharov and V.E. Zakharov, “Lump chains in the KP-I equation”, Stud. Appl. Math. (2021), doi.org/10.1111/sapm.12420 (see also arXiv:2102.07038).
  * (23) J. Rao, K.W. Chow, D. Mihalache and J.S. He, “Completely resonant collision of lumps and line solitons in the Kadomtsev–Petviashvili I equation”, Stud. Appl. Math. 2021 (doi.org/10.1111/sapm.12417).
  * (24) B. Yang and J. Yang, “Rogue wave patterns in the nonlinear Schrodinger equation”, Physica D 419, 132850 (2021).
  * (25) B. Yang and J. Yang, “Universal rogue wave patterns associated with the Yablonskii-Vorob’ev polynomial hierarchy”, Physica D 425, 132958 (2021).
  * (26) J. Weiss, “Modified equations, rational solutions, and the Painlevé property for the Kadomtsev–Petviashvili and Hirota–Satsuma equations”, J. Math. Phys. 26, 2174 (1985).
  * (27) Y. Ohta, J. Yang, “General high-order rogue waves and their dynamics in the nonlinear Schrödinger equation”, Proc. Roy. Soc. A 468, 1716 (2012).
  * (28) B. Yang and J. Yang, “General rogue waves in the three-wave resonant interaction systems”, IMA J. Appl. Math. 86, 378-425 (2021).
  * (29) B. Yang and J. Yang, “General rogue waves in the Boussinesq equation”, J. Phys. Soc. Jpn. 89, 024003 (2020).
  * (30) J. Chen, Y. Chen, B.F. Feng, K. Maruno and Y. Ohta, “General high-order rogue waves of the (1+1)-dimensional Yajima–Oikawa system”, J. Phys. Soc. Jpn. 87, 094007 (2018).
  * (31) A.I. Yablonskii, Vesti Akad. Navuk. BSSR Ser. Fiz. Tkh. Nauk. 3, 30 (1959) (in Russian).
  * (32) A.P. Vorob’ev, “On rational solutions of the second Painlevé equation”, Diff. Eqns. 1, 58 (1965).
  * (33) K. Kajiwara and Y. Ohta, “Determinant structure of the rational solutions for the Painlevé II equation”, J. Math. Phys. 37, 4693 (1996).
  * (34) P.A. Clarkson and E.L. Mansfield, “The second Painlevé equation, its hierarchy and associated special polynomials”, Nonlinearity 16, R1 (2003).
  * (35) S. Fukutani, K. Okamoto and H. Umemura, “Special polynomials and the Hirota bilinear relations of the second and the fourth Painlevé equations”, Nagoya Math. J. 159, 179-200 (2000).
  * (36) M. Taneda, “Remarks on the Yablonskii-Vorob’ev polynomials”, Nagoya Math. J. 159, 87-111 (2000).
  * (37) R.J. Buckingham and P.D. Miller, “Large-degree asymptotics of rational Painlevé-II functions: noncritical behaviour”, Nonlinearity 27, 2489 (2014).
  * (38) F. Balogh, M. Bertola and T. Bothner, “Hankel determinant approach to generalized Vorob’ev-Yablonski polynomials and their roots”, Constr. Approx. 44, 417 (2016).
  * (39) G. Felder, A.D. Hemery and A.P. Veselov, “Zeros of Wronskians of Hermite polynomials and Young diagrams”, Physica D 241, 2131-2137 (2012).
  * (40) M. Garca-Ferrero and D. Gmez-Ullate, “Oscillation theorems for the Wronskian of an arbitrary sequence of eigenfunctions of Schrödinger’s equation”, Lett. Math. Phys. 105, 551-573 (2015).
  * (41) N. Bonneux, C. Dunning and M. Stevens, “Coefficients of Wronskian Hermite polynomials”, Stud. Appl. Math. 144, 245–288 (2020).
  * (42) J. Dong, L. Ling and X. Zhang, “Kadomtsev-Petviashvili equation: one-constraint method and lump pattern”, arXiv:2108.09715 [nlin.SI] (2021).
  * (43) R. Hirota, _The direct method in soliton theory_ (Cambridge University Press, Cambridge, 2004).