1. Introduction
The Kadomtsev–Petviashvili (KP) equation is a prototypical non-linear wave equation in 2+1 dimensions. It is a universal model for weakly non-linear waves in the long wavelength regime and it arises in many different fields of applications ranging from plasmas [Reference Fogaca, Navarra and Ferreira Filho26, Reference Kadomtsev and Petviashvili31] to water waves [Reference Ablowitz and Segur6], cosmology [Reference Gao and Tian27, Reference Gao and Tian28], ferromagnetics [Reference Leblond35] and elasticity [Reference Berjamin, Destrade and Saccomandi11, Reference Erbay23]. The KP equation is also the prototypical completely integrable system in 2+1 dimensions [Reference Novikov, Manakov, Pitaevskii and Zakharov42], and as such it possesses a deep mathematical structure and a variety of exact solutions. The KP equation comes in two variants, labelled respectively KPI and KPII. Both variants admit multi-soliton solutions, which are known as line solitons since the solution remains constant along certain directions at space infinity. The initial value problem for the KP equation is also amenable to exact treatment via the inverse scattering transform (IST) [Reference Ablowitz and Clarkson2, Reference Ablowitz and Fokas4, Reference Ablowitz and Fokas5, Reference Ablowitz and Segur7]. Over the last 40 years, the IST has been extended to initial conditions that are a small perturbation of an exact soliton solution [Reference Boiti, Pempinelli and Pogrebkov15, Reference Boiti, Pempinelli, Pogrebkov and Prinari16, Reference Villarroel and Ablowitz52, Reference Wu54]. The IST with nonzero boundary conditions has still not been used effectively to study the behaviour of solutions, for example, to determine the long-time dynamics. Moreover, no IST formulation exists at present for more general classes of non-decaying initial conditions.
For the KPII equation, the recent generalization of Whitham modulation theory [Reference Whitham53] to systems in 2+1 spatial dimensions [Reference Ablowitz, Biondini and Wang1] made it possible to study a variety of dynamical problems involving partial solitons and bent solitons, as well as interactions between these objects and a mean flow [Reference Ryskamp, Hoefer and Biondini46–Reference Ryskamp, Hoefer and Biondini48]. In the case of the KPI equation, the situation is different since its line solitons are unstable to transverse perturbations [Reference Alexander, Pego and Sachs10, Reference Pelinovsky and Stepanyants45]. On the other hand, the KPI equation admits stable, fully two-dimensional (2D) localized rational solutions called lumps [Reference Chakravarty17–Reference Chakravarty and Zowada19, Reference Ma38, Reference Manakov, Zakharov, Bordag, Its and Matveev40, Reference Minzoni and Smythe41, Reference Satsuma and Ablowitz49]. More general solutions such as higher-order lumps [Reference Ablowitz and Villarroel8, Reference Krichever33, Reference Pelinovsky44, Reference Yang and Yang55] and lump chains [Reference Lester, Gelash, Zakharov and Zakharov36, Reference Pelinovsky and Stepanyants45] have also been derived and studied. In fact, exact solutions of the KPI equation comprised of line solitons, lumps and lump chains can be obtained [Reference Stepanyants, Zakharov and Zakharov50], typically using Hirota’s bilinear method [Reference Hirota30]. Note that lumps have recently been observed experimentally for the first time in the realm of non-linear optics [Reference Dieli, Pierangeli, Baronio, Trillo and Conti20].
A natural question is therefore whether it is possible to formulate a modulation theory for lumps in the same spirit of the one that was developed in [Reference Ablowitz, Biondini and Wang1] for periodic solutions. In this work, we answer this question in the affirmative. Specifically, in Section 2, we derive a four-component (2+1)-dimensional hyperbolic system that governs modulations of a lump in the presence of a slowly varying mean field. In Section 3, we then demonstrate that the resulting system satisfies necessary conditions for integrability and show that it can be diagonalized in certain (1+1)-dimensional reductions. Moreover, in Section 4, we show how the modulation system can be used to effectively study quantitatively a concrete dynamical scenario, namely an interaction between a lump and a rarefaction wave. We use the modulation equations to predict the slow evolution of the amplitude and velocity of the lump during the interaction, as well as the values of the parameters of a lump that fully transmits through the rarefaction wave. Finally, Section 5 ends this work with a discussion and some concluding remarks. The appendix contains a discussion of the numerical methods used for the computer simulations.
2. Derivation of the lump modulation system
We consider the KPI equation in evolutionary form
\begin{equation}
u_{t}+uu_{x}+u_{xxx}-v_{y}=0,\qquad v_{x}=u_{y}.
\end{equation}The exact single-lump solution of (1) is given by
\begin{equation}
u(x,y,t)=\bar{u}+U(\eta,\theta),\qquad U(\eta,\theta)=24\frac{-(\eta+q\theta)^2+a^2\theta^2+3/a^2}{\big[(\eta+q\theta)^2+a^2\theta^2+3/a^2\big]^2},
\end{equation}
\begin{equation}
v(x,y,t)=\bar{v}+V(\eta,\theta),\qquad V(\eta,\theta)=\int_{-\infty}^{\eta}U_{\theta}(\eta',\theta)\,\mathrm{d}\eta',
\end{equation}
\begin{equation}
\eta=x-(\bar{u}+a^2+q^2)t + \eta_0,\qquad\theta=y+2qt+\theta_0.
\end{equation}
$(\bar{u},\bar{v})$ represents the background on which the lump is placed, and is arbitrary due to the Galilean invariance of the KP equation. The amplitude 
$A$ and velocity c of the lump are
\begin{equation}
A=U(0,0) = 8a^2,\qquad \textbf{c}=\begin{bmatrix}
\bar{u}+a^2+q^2\\
-2q
\end{bmatrix}.
\end{equation} Summarizing, the solution (1) is completely determined by four independent parameters 
$\bar u$, 
$\bar v$, 
$a$, 
$q$ and its origin 
$(\eta_0,\theta_0)$. Next, we are interested in considering solutions in which all of these parameters are slowly varying, and we aim to derive a system that governs their modulations, which will allow us to describe the dynamics of a lump subject to a slowly varying mean field. To accomplish this, we begin by deriving useful integral identities similar to one derived in [Reference Ablowitz, Cole, El, Hoefer and Luo3] for the Korteweg–de Vries (KdV) equation. Subsequently, we invoke a multiple scale ansatz consisting of a slowly varying lump and mean field, and then follow the lump as a fluid parcel is followed in Lagrangian coordinates.
We begin by introducing the splitting 
$u=\bar{u}+U$, 
$v=\bar{v}+V$, where at this stage all quantities are assumed to depend explicitly on space and time. At a later stage, we will make the assumption that the mean field 
$\bar{u},\bar{v}$ and the lump parameters 
$a$ and 
$q$ are slowly varying. Upon substitution, (1) becomes
\begin{equation}
U_{t}+UU_{x}+U_{xxx}-V_{y}=-(\bar{u} U)_{x}-F[\bar{u},\bar{v}],\qquad F[\bar{u},\bar{v}]=\bar{u}_{t}+\bar{u}\bar{u}_{x}+\bar{u}_{xxx}-\bar{v}_{y},
\end{equation}
\begin{equation}
V_{x}-U_{y}=G[\bar{u},\bar{v}],\qquad G[\bar{u},\bar{v}]=\bar{u}_{y}-\bar{v}_{x}.
\end{equation}
$U$ and 
$V$ are localized, multiplying (4a) by 
$U$ and integrating over all space gives
\begin{equation}
\frac{\mathrm{d}}{\mathrm{d} t}\iint_{\mathbb{R}^{2}}U^{2}\,\mathrm{d} x\,\mathrm{d} y=-\iint_{\mathbb{R}^{2}}U^{2}\bar{u}_{x}\,\mathrm{d} x\,\mathrm{d} y-2\iint_{\mathbb{R}^{2}}UF[\bar{u},\bar{v}]\,\mathrm{d} x\,\mathrm{d} y.
\end{equation} If 
$\overline u$ and 
$\overline v$ are constant, this simply corresponds to conservation of horizontal momentum. We now aim to derive a similar identity for the transverse momentum. To this end, from (4b), we get
\begin{equation}
\iint_{\mathbb{R}^{2}}UV_{t}\,\mathrm{d} x\,\mathrm{d} y=\iint_{\mathbb{R}^{2}}U\partial_{x}^{-1}U_{yt}\,\mathrm{d} x\,\mathrm{d} y+\iint_{\mathbb{R}^{2}}U\partial_{x}^{-1}G[\bar{u},\bar{v}]_{t}\,\mathrm{d} x\,\mathrm{d} y,
\end{equation}with a suitable definition of the antiderivative operator. Integrating by parts in 
$y$ in the first term on the right-hand side, substituting 
$U_{y}=V_{x}-G[\bar{u},\bar{v}]$, and then integrating by parts in 
$x$ in both terms yields
\begin{equation}
\iint_{\mathbb{R}^{2}}UV_{t}\,\mathrm{d} x\,\mathrm{d} y=\iint_{\mathbb{R}^{2}}VU_{t}\,\mathrm{d} x\,\mathrm{d} y+\iint_{\mathbb{R}^{2}}\big(G[\bar{u},\bar{v}]\partial_{x}^{-1}U_{t}-G[\bar{u},\bar{v}]_{t}\partial_{x}^{-1}U\big)\,\mathrm{d} x\,\mathrm{d} y.
\end{equation} Next, multiplying (4a) by 
$V$ and integrating we also find
\begin{equation}
\iint_{\mathbb{R}^{2}}VU_{t}\,\mathrm{d} x\,\mathrm{d} y+\iint_{\mathbb{R}^2}(VUU_{x}+VU_{xxx})\,\mathrm{d} x\,\mathrm{d} y=-\frac{1}{2}\iint_{\mathbb{R}^2}U^2\bar{u}_{y}\,\mathrm{d} x\,\mathrm{d} y-\iint_{\mathbb{R}^2}VF[\bar{u},\bar{v}]\,\mathrm{d} x\,\mathrm{d} y.
\end{equation}Putting (7) and (8) together and making use of
\begin{equation}
\iint_{\mathbb{R}^{2}}VUU_{x}\,\mathrm{d} x\,\mathrm{d} y=-\frac{1}{2}\iint_{\mathbb{R}^2}U^{2}G[\bar{u},\bar{v}]\,\mathrm{d} x\,\mathrm{d} y,\qquad\iint_{\mathbb{R}^{2}}VU_{xxx}\,\mathrm{d} x\,\mathrm{d} y=-\iint_{\mathbb{R}^{2}}U_{xx}G[\bar{u},\bar{v}]\,\mathrm{d} x\,\mathrm{d} y,
\end{equation}we arrive at the identity
\begin{align}
\frac{\mathrm{d}}{\mathrm{d} t}\iint_{\mathbb{R}^{2}}UV\,\mathrm{d} x\,\mathrm{d} y&=-\iint_{\mathbb{R}^{2}}U^{2}\bar{u}_{y}\,\mathrm{d} x\,\mathrm{d} y-2\iint_{\mathbb{R}^{2}}VF[\bar{u},\bar{v}]\,\mathrm{d} x\,\mathrm{d} y\\ \nonumber &\quad -\iint_{\mathbb{R}^{2}}G[\bar{u},\bar{v}]\big(\partial_{x}^{-1}V_{y}-U\bar{u}-\partial_{x}^{-1}F[\bar{u},\bar{v}]\big)\,\mathrm{d} x\,\mathrm{d} y-\iint_{\mathbb{R}^{2}}G[\bar{u},\bar{v}]_{t}\partial_{x}^{-1}U\,\mathrm{d} x\,\mathrm{d} y.\end{align} Note that for slowly varying 
$F$, 
$\partial_{x}^{-1}F$ is not necessarily well-defined, but this term will be of higher order in the multi-scale analysis we now initiate.
We now take 
$U$ and 
$V$ to be the lump solution as in (2a), (2b) and (2c) and allow the lump parameters to vary slowly with space and time. That is, introducing a small parameter 
$0 \lt \varepsilon\ll1$, we define the slow spatial and temporal scales 
$X=\varepsilon x$, 
$Y=\varepsilon y$, 
$T=\varepsilon t$ and
\begin{align}
\bar{u}&=\bar{u}(X,Y,T),\qquad \bar{v}=\bar{v}(X,Y,T),\qquad a=a(X,Y,T),\qquad q=q(X,Y,T),
\end{align}
\begin{align}
\eta&=x-\frac{1}{\varepsilon}\int_{0}^{\varepsilon t}\big[\bar{u}(X,Y,T')+a(X,Y,T')^2+q(X,Y,T')^2\big]\,\mathrm{d} T',\qquad\theta=y+\frac{2}{\varepsilon}\int_{0}^{\varepsilon t}q(X,Y,T')\,\mathrm{d} T'.
\end{align}
\begin{equation}
u(x,y,t)=\bar{u}(X,Y,T)+U(\eta,\theta,X,Y,T),\qquad v(x,y,t)=\bar{v}(X,Y,T)+V(\eta,\theta,X,Y,T),
\end{equation}which represents a modulated lump propagating on a slowly varying mean. Accordingly, we have the operator expansions
\begin{equation}
\partial_{x}=\partial_{\eta}+\varepsilon\partial_{X},\qquad\partial_{x}^{-1}=\partial_{\eta}^{-1}(1-\varepsilon\partial_{\eta}^{-1}\partial_{X}+\cdots),
\end{equation}
\begin{equation}
\partial_{y}=\partial_{\theta}+\varepsilon\partial_{Y},\qquad\partial_{t}=-(\bar{u}+a^2+q^2)\partial_{\eta}+2q\partial_{\theta}+\varepsilon\partial_{T}.
\end{equation}
$(X,Y,T)$. Thus, the total time derivative should be interpreted as the convective derivative
\begin{equation}
\frac{\mathrm{d}}{\mathrm{d} t}=\varepsilon(\partial_{T}+\mathbf{c}\cdot\nabla),\qquad\nabla=\begin{bmatrix}
\partial_{X}\\\partial_{Y}
\end{bmatrix}.
\end{equation} Inserting this ansatz and keeping terms of 
$\mathcal{O}(\varepsilon)$, we find that in this multi-scale setting the generic identities (5) and (10) provide the following relations:
\begin{equation}
(\partial_{T}+\mathbf{c}\cdot\nabla)\iint_{\mathbb{R}^2}U^2\,\mathrm{d}\eta\,\mathrm{d}\theta+\bar{u}_{X}\iint_{\mathbb{R}^{2}}U^2\,\mathrm{d}\eta\,\mathrm{d}\theta=-2(\bar{u}_{T}+\bar{u}\bar{u}_{X}-\bar{v}_{Y})\iint_{\mathbb{R}^2}U\,\mathrm{d}\eta\,\mathrm{d}\theta,
\end{equation}
\begin{align}
(\partial_{T}+\mathbf{c}\cdot\nabla)\iint_{\mathbb{R}^2}UV\,\mathrm{d}\eta\,\mathrm{d}\theta+\bar{u}_{Y}\iint_{\mathbb{R}^{2}}U^2\,\mathrm{d}\eta\,\mathrm{d}\theta&=-2(\bar{u}_{T}+\bar{u}\bar{u}_{X}-\bar{v}_{Y})\iint_{\mathbb{R}^2}V\,\mathrm{d}\eta\,\mathrm{d}\theta \\ \nonumber
&\quad +\bar{u}(\bar{u}_{Y}-\bar{v}_{X})\iint_{\mathbb{R}^{2}}U\,\mathrm{d}\eta\,\mathrm{d}\theta.\end{align}
$U$ over all space is only conditionally convergent. To remove this ambiguity, constraints must be imposed on the mean field parameters. In particular, enforcing that the coefficients of the conditionally convergent integrals on the right-hand sides of both parts of (15a) and (15b) vanish, the mean field must satisfy the dispersionless KP (dKP) equation:
\begin{equation}
\bar{u}_{T}+\bar{u}\bar{u}_{X}-\bar{v}_{Y}=0,\qquad\bar{v}_{X}=\bar{u}_{Y}.
\end{equation} Moreover, the integrals in (10) involving the antiderivative operator 
$\partial_{x}^{-1}$ do not contribute at this order. The only remaining relevant integrals are well-defined and are found to be
\begin{equation}
\iint_{\mathbb{R}^{2}}U^{2}\,\mathrm{d}\eta\,\mathrm{d}\theta=96\pi a,\qquad\iint_{\mathbb{R}^{2}}UV\,\mathrm{d}\eta\,\mathrm{d}\theta=96\pi aq.
\end{equation}Putting these into (15a) and (15b) gives the modulation equations for the lump parameters:
\begin{equation}
(\partial_{T}+\mathbf{c}\cdot\nabla)a+a\bar{u}_{X}=0,
\end{equation}
\begin{equation}
(\partial_{T}+\mathbf{c}\cdot\nabla)(aq)+a\bar{u}_{Y}=0.
\end{equation}
\begin{equation}
\mathsf{E}_3 \mathbf{v}_{T}+\mathsf{A}\mathbf{v}_{X}+\mathsf{B}\mathbf{v}_{Y}=0,\qquad
\mathbf{v}=\big[\,a\,,\, q\,,\, \bar{u}\,,\, \bar{v}\,\big]^\top,
\end{equation}where the superscript 
$\top$ denotes matrix transpose, 
$\mathsf{E}_3 = \mathop{\rm diag}\nolimits(1,1,1,0),$ and
\begin{equation}
\mathsf{A}=\begin{bmatrix}\bar{u}+a^2+q^2&0&a&0\\
0&\bar{u}+a^2+q^2&-q&0\\
0&0&\bar{u}&0\\
0&0&0&1\end{bmatrix},\qquad\mathsf{B}=\begin{bmatrix}
-2q&0&0&0\\
0&-2q&1&0\\
0&0&0&-1\\
0&0&-1&0
\end{bmatrix}.
\end{equation} Note how the system (19a and 19b) is partially decoupled, since (16) does not depend on 
$a$ and 
$q$. Conversely, (18a) and (18b) show that the dynamics of the lump is influenced by the modulation of the mean field. The system (19a and 19b) is one of the main results of this work. In the upcoming sections, we study some of the properties of the system and we also show how the system can be used to quantitatively describe the interactions between lumps and a mean field.
The characteristic speeds 
$\lambda_j$ for wave propagation in the direction 
$\lbrack\,\cos\phi\,,\,\sin\phi\,\rbrack^\top$ of the lump-mean field modulation system (19a and 19b) satisfy 
$\left|\cos\phi A+\sin\phi B-\lambda_jE_3\right|=0$ and are
\begin{equation}\lambda_1=\cos\phi\overline u-\sin\phi\tan\phi,\quad\lambda_{2,3}=\cos\phi(\overline u+a^2+q^2)-2q\sin\phi.\end{equation} The double characteristic 
$\lambda_{2,3}$ corresponds to the horizontal (
$\phi=0$) and vertical (
$\phi=\pi/2$) lump velocity 
$\mathbf{c}$ components (3). This characteristic has two linearly independent eigenvectors 
$[\,1\,,\,0\,,\,0\,,\,0\,]^\top$ and 
$[\,0\,,\,1\,,\,0\,,\,0\,]^\top$. The characteristic speed 
$\lambda_1$ corresponds to evolution of 
$(\bar{u},\bar{v})$ according to the dKP equation (16) that includes the fourth equation in (19a) and (19b), which is a constraint rather than an evolution equation. Consequently, the system (19a and 19b) is hyperbolic on the appropriate constraint manifold in the sense of [Reference Ferapontov and Khusnutdinova25].
It is worth commenting on how the modulation equations (18a) and (18b) should be interpreted, since at first sight their partial differential equation (PDE) form might seem incongruent with a lump being a localized, finite-energy object. Indeed, a single lump is described by two quantities 
$(a,q)$ evolving as a function of 
$T$ along its trajectory, and one might have therefore expected a system of ordinary differential equations for its parameters rather than the PDEs (18a and 18b). The resolution of this apparent paradox is that (18a) and (18b) should not be seen as the equations of motion of a single lump, but rather as the equations of motion of a slowly varying, non-interacting lump field defined throughout the spatial domain, of which the trajectory of any individual lump is defined by the characteristic 
$\lambda_{2,3}$. With this interpretation, the reason that 
$\lambda_{2,3}$ is a double characteristic becomes clearer. Because the lump is localized, its two degrees of freedom must propagate together, i.e., they cannot be independently modulated. The dynamics of a single lump is recovered by restricting (18a) and (18b) to a characteristic curve 
$(X(T),Y(T))$, as will be demonstrated in Sections 3 and 4. Off this curve, the lump field is to be interpreted as describing a dilute distribution of lumps with the same local parameters. This is a useful fiction for the purpose of modulation theory, but one that becomes a literal description of the solution if a lump is actually placed at the corresponding point in space.
This interpretation is analogous to the soliton limit of the Whitham modulation equations for the KdV equation, which is obtained as the 
$k \to 0$ degeneration of the genus-one KdV–Whitham system and which describes the slow evolution of the parameters of a soliton placed on a slowly varying mean flow [Reference Ablowitz, Cole, El, Hoefer and Luo3, Reference El and Hoefer22]. As noted in [Reference Ablowitz, Cole, El, Hoefer and Luo3, Reference Maiden, Anderson, Franco, El and Hoefer39], the resulting hyperbolic system has the dual interpretation of governing either a single soliton along its trajectory or a dilute (non-interacting) soliton field distributed throughout the spatial domain, the two interpretations being consistent because the soliton’s parameter propagates along a double characteristic (the soliton speed) of the modulation system. The double characteristic splits when modulation requires a nonzero cnoidal wavenumber 
$k \gt 0$, i.e., the soliton train becomes interacting. The splitting of the double characteristic due to modulation is one of the foundational results of Whitham modulation theory [Reference Whitham53] and is precisely what occurs in a dispersive shock wave [Reference El and Hoefer22]. The system (19a and 19b) plays the same role for KPI lumps, with the additional richness that both the lump and its mean-field background are now genuinely 2D. Splitting of the lump-mean field modulation system’s double characteristic (20) will occur when the lump field becomes interacting, i.e., in a lump chain or lump lattice.
At the same time, we should also point out one important difference between the KdV case and the present one: in the former, the soliton-limit modulation system arises as a degenerate reduction of an underlying Whitham system for the periodic cnoidal wavetrain [Reference El and Hoefer22], whereas, in the present setting, equations (18a) and (18b) are derived directly from integral identities for the localized lump itself, without any underlying periodic structure.
3. Integrability test, integration by characteristics and 2D reductions
Since the KP equation is completely integrable, one expects that exact asymptotic reductions of it will inherit the integrability properties. This was already shown to be true for various reductions of the Whitham modulation equations for the periodic travelling wave solutions of the KP equation in [Reference Biondini, Bivolcic, Hoefer and Moro12, Reference Biondini, Hoefer and Moro13]. In this section, we show that a similar result applies for the modulation system (19a and 19b).
3.1. Haantjes tensor test for integrability
For a (2+1)-dimensional quasilinear system of the form (19a and 19b), it was shown in [Reference Ferapontov and Khusnutdinova25] that a necessary condition for integrability is that the matrix
\begin{equation}
\mathsf{M}=(k\mathsf{E}_{3}+\mathsf{A})^{-1}(l\mathsf{E}_{3}+\mathsf{B}),
\end{equation}passes the Haantjes tensor test [Reference Haantjes29] for arbitrary 
$k$ and 
$l$. Denoting the entries of 
$\mathbf{v}$ in eq. (19a) as 
$v^{i}$ and the entries of 
$\mathsf{M}$ as 
$M_{j}^{i}$, one first defines the Nijenhuis tensor as
\begin{equation}
N_{jk}^{i}=M_{j}^{p}\partial_{v^{p}}M_{k}^{i}-M_{k}^{p}\partial_{v^{p}}M_{j}^{i}-M_{p}^{i}\left(\partial_{v^{j}}M_{k}^{p}-\partial_{v^{k}}M_{j}^{p}\right).
\end{equation} Then, the matrix 
$\mathsf{M}$ is said to pass the Haantjes tensor test if its Haantjes tensor, defined by
\begin{equation}
H_{jk}^{i}=N_{pr}^{i}M_{j}^{p}M_{k}^{r}-N_{jr}^{p}M_{p}^{i}M_{k}^{r}-N_{rk}^{p}M_{p}^{i}M_{j}^{r}+N_{jk}^{p}M_{r}^{i}M_{p}^{r},
\end{equation}is identically zero. Straightforward computations with a computer algebra system demonstrate that this is indeed the case for the system (19a and 19b). This suggests that any 2D reduction of (19a) and (19b) can be diagonalized, as will be demonstrated for the 
$XT$, 
$YT$ and 
$XY$ dependent cases in following subsections. It also suggests that the full system (19a and 19b) is integrable, as will be demonstrated in the next subsection.
3.2. Integration of the lump modulation equations by the method of characteristics
Recall that the evolution of the mean field is decoupled from the presence of the lump and is governed by the dKP equation (16), which is completely integrable. Here, we, therefore, fix 
$\bar{u}$ and 
$\bar{v}$ to be a generic solution of the dKP equation (16). We then write the remaining two modulation equations in the form
\begin{equation}
a_{T}+(\bar{u}+a^2+q^2)a_{X}-2qa_{Y}=-a\bar{u}_{X},
\end{equation}
\begin{equation}
(aq)_{T}+(\bar{u}+a^2+q^2)(aq)_{X}-2q(aq)_{Y}=-a\bar{u}_{Y}.
\end{equation}
$a(T)=a(X(T),Y(T),T)$, 
$q(T)=q(X(T),Y(T),T)$ and the characteristic trajectories 
$X(T)$, 
$Y(T)$:
\begin{equation}
\frac{\mathrm{d} a}{\mathrm{d} T}=-a\bar{u}_{X},\qquad\frac{\mathrm{d} (aq)}{\mathrm{d} T}=-a\bar{u}_{Y},\qquad \frac{\mathrm{d} X}{\mathrm{d} T}=\bar{u}+a^2+q^2,\qquad\frac{\mathrm{d} Y}{\mathrm{d} T}=-2q.
\end{equation} It is worth noting that from (25) it is clear that if the mean field is independent of 
$X$, then 
$a$ (the horizontal momentum of the lump) is a constant of motion, and if the mean field is independent of 
$Y$, then 
$aq$ (the transverse momentum of the lump) is a constant of motion. Furthermore, for any stationary mean 
$\bar{u}=\bar{u}(X,Y)$, (25) admits a conserved Hamiltonian
\begin{equation}
\mathcal{H}(X,Y,a,aq)=a\bar{u}(X,Y)+\frac{1}{3}a^{3}-\frac{(aq)^2}{a},
\end{equation}from which the system is obtained through
\begin{equation}
\frac{\mathrm{d} a}{\mathrm{d} T}=-\frac{\partial\mathcal{H}}{\partial X},\qquad\frac{\mathrm{d} (aq)}{\mathrm{d} T}=-\frac{\partial\mathcal{H}}{\partial Y},\qquad\frac{\mathrm{d} X}{\mathrm{d} T}=\frac{\partial\mathcal{H}}{\partial a},\qquad\frac{\mathrm{d} Y}{\mathrm{d} T}=\frac{\partial\mathcal{H}}{\partial (aq)}.
\end{equation}The equations Equation 25 will be used in a specific situation in Section 4 to describe the interaction of a lump with a specified mean flow.
3.3. 2D reductions of the modulation equations
In this subsection, we consider the (1+1)- and (2+0)-dimensional reductions of the full modulation system (19a and 19b) that are obtained when one can neglect variations with respect to one of the independent variables. Moreover, we show that each reduction can be diagonalized in terms of Riemann invariants.
3.3.1. 
$XT$-dependent reduction
First, in the case of modulations that are independent of 
$Y$, we have that 
$\bar{v}$ is constant, and (19a) and (19b) reduce to the three-component system
\begin{equation}
a_{T}+(\bar{u}+a^2+q^2)a_{X}+a\bar{u}_{X}=0,
\end{equation}
\begin{equation}
(aq)_{T}+(\bar{u}+a^2+q^2)(aq)_{X}=0,
\end{equation}
\begin{equation}
\bar{u}_{T}+\bar{u}\bar{u}_{X}=0.
\end{equation}
\begin{equation}
\tilde{\mathbf{v}}_{T}+\tilde{\mathsf{A}}\tilde{\mathbf{v}}_{X}=0,\qquad\tilde{\mathbf{v}}=\begin{bmatrix}
a\\q\\\bar{u}
\end{bmatrix},\qquad
\tilde{\mathsf{A}}=\begin{bmatrix}
\bar{u}+a^{2}+q^{2}&0&a\\
0&\bar{u}+a^2+q^2&-q\\
0&0&\bar{u}
\end{bmatrix}.
\end{equation} The eigenvalues of the matrix 
$\tilde{\mathsf{A}}$ are 
$\lambda_{1}=\lambda_{2}=\bar{u}+a^2+q^2$ and 
$\lambda_{3}=\bar{u}$ with respective left eigenvectors
\begin{equation}l_1=\begin{bmatrix}q\\a\\0\end{bmatrix},\qquad l_2=\begin{bmatrix}{(a^2+q^2)}/a\\0\\1\end{bmatrix},\qquad l_3=\begin{bmatrix}0\\0\\1\end{bmatrix}.\end{equation} From (28b), we have that 
$R_{1}=aq$ is a Riemann invariant with characteristic speed 
$\bar{u}+a^2+q^2$. (This is also apparent from 
$\mathrm{d} R_{1}=l_{1}\cdot \mathrm{d}\tilde{\mathbf{v}}$.) The other Riemann invariant with the same characteristic speed is found from 
$\mathrm{d} R_{2}=l_{2}\cdot \mathrm{d}\tilde{\mathbf{v}}$, which yields
\begin{equation}
R_{2}=\int\frac{a^2+q^2}{a}\,\mathrm{d} a+\bar{u}.
\end{equation} Since 
$q=R_{1}/a$, we find
\begin{equation}
R_{2}=\frac{1}{2}a^2-\frac{R_{1}^{2}}{2a^{2}}+\bar{u}=\frac{a^{2}-q^2}{2}+\bar{u}.
\end{equation} Finally, from (28c) the third Riemann invariant is 
$R_{3}=\bar{u}$. In terms of the Riemann invariants, the characteristic speeds are given by
\begin{equation}
\lambda_{1}=\lambda_{2}=R_{3}+2\sqrt{(R_{2}-R_{3})^{2}+R_{1}^{2}},\qquad\lambda_{3}=R_{3}.
\end{equation}Thus, the diagonal form of (29) can be written compactly as
\begin{equation}
\partial_{T}R_{i}+\lambda_{i}\partial_{X}R_{i}=0,\qquad i=1,2,3.
\end{equation}Note the system (29) is not genuinely non-linear nor strictly hyperbolic, since two of the characteristic speeds coincide.
YT-dependent reduction.
In the case of modulations that are independent of 
$X$, there is no reduction in the number of dependent variables, but the system (19a and 19b) reduces to
\begin{align}
a_{T}-2qa_{Y}&=0,
\end{align}
\begin{align}
q_{T}-2qq_{Y}+\bar{u}_{Y}&=0,
\end{align}
\begin{align}
\bar{u}_{T}-\bar{v}_{Y}&=0,
\end{align}
\begin{align}
\bar{u}_{Y}&=0.
\end{align}
\begin{equation}
\bar{u}=f(T),\qquad\bar{v}=f'(T)Y+g(T),
\end{equation}for some functions 
$f$ and 
$g$ of only 
$T$. Then, the remaining system shows that the lump parameters satisfy
\begin{equation}
\begin{bmatrix}
a\\q
\end{bmatrix}_{T}+\begin{bmatrix}
-2q&0\\0&-2q
\end{bmatrix}\begin{bmatrix}
a\\q
\end{bmatrix}_{Y}=0,
\end{equation}which is already in diagonal form. Similarly to the 
$XT$-reduction (29), the system (37) is not genuinely non-linear nor strictly hyperbolic, since the two characteristic speeds coincide.
XY-dependent reduction.
Finally, we consider the case of 
$T$-independent modulations, in which the system can be written for 
$q\neq0$ as
\begin{equation}
\mathbf{v}_{Y}+\mathsf{C}\mathbf{v}_{X}=0,\qquad\mathsf{C}=\mathsf{B}^{-1}\mathsf{A}=\begin{bmatrix}
-\frac{\bar{u}+a^2+q^2}{2q}&0&-\frac{a}{2q}&0\\
0&-\frac{\bar{u}+a^2+q^2}{2q}&\frac{1}{2}&-\frac{1}{2q}\\
0&0&0&-1\\
0&0&-\bar{u}&0
\end{bmatrix}.
\end{equation} The eigenvalues of 
$\mathsf{C}$ are found to be
\begin{equation}
\lambda_{1}=\lambda_{2}=-\frac{\bar{u}+a^2+q^2}{2q},\qquad\lambda_{3}=\sqrt{\bar{u}},\qquad\lambda_{4}=-\sqrt{\bar{u}},
\end{equation}with corresponding left eigenvectors
\begin{equation}l_1=\begin{bmatrix}a^2-q^2+\overline u\\-2aq\\a\\0\end{bmatrix},\qquad l_2=\begin{bmatrix}2aq\\a^2-q^2+\overline u\\q\\1\end{bmatrix},\qquad l_3=\begin{bmatrix}0\\0\\-\sqrt{\overline u}\\1\end{bmatrix},\qquad l_4=\begin{bmatrix}0\\0\\\sqrt{\overline u}\\1\end{bmatrix}.\end{equation} If 
$\bar{u} \lt 0$, the system (38) becomes elliptic. We therefore assume that 
$\bar{u}\geq0$ to ensure hyperbolicity. The left eigenvectors spanning the eigenspace of the double eigenvalue have been chosen such that it is straightforward to determine the related Riemann invariants by integration. In particular, from 
$\mathrm{d} R_{1}= l_{1}\cdot\mathrm{d} \mathbf{v}$ and 
$\mathrm{d} R_{2}= l_{2}\cdot\mathrm{d} \mathbf{v}$ we find
\begin{equation}
R_{1}=\frac{1}{3}a^3-aq^2+a\bar{u},\qquad R_{2}=a^{2}q-\frac{1}{3}q^{3}+q\bar{u}+\bar{v}.
\end{equation} Note that 
$R_{1}$ is the Hamiltonian 
$\mathcal{H}$ as in (26), which is generically a conserved quantity of the full modulation system when the mean field is time-independent. The Riemann invariants corresponding to 
$\pm\sqrt{\bar{u}}$ are found to be
\begin{equation}
R_{3}=\frac{2}{3}\bar{u}^{3/2}-\bar{v},\qquad R_{4}=-\frac{2}{3}\bar{u}^{3/2}-\bar{v}.
\end{equation}From (42), we can express the corresponding characteristic speeds as
\begin{equation}
\lambda_{3}=\left[\frac{3}{4}(R_{3}-R_{4})\right]^{1/3},\qquad\lambda_{4}=-\left[\frac{3}{4}(R_{3}-R_{4})\right]^{1/3}.
\end{equation} On the other hand, to determine the double characteristic speed 
$\lambda_{1}=\lambda_{2}$ in terms of the Riemann invariants, one needs to invert the cubic relations in (41), which is possible but involves cumbersome formulae. The diagonal system can then be expressed as
\begin{equation}
\partial_{Y}R_{i}+\lambda_{i}\partial_{X}R_{i}=0,\qquad i=1,2,3,4.
\end{equation}Note that, once again, the system (38) is not genuinely non-linear or strictly hyperbolic, because two of the characteristic speeds coincide.
We should briefly comment on how the system (38) should be interpreted, since time-independent equations may at first appear to be incongruous with the lump being inherently a time-dependent entity. One should realize, however, that (38) describes the modulations of a lump field rather than simply an individual lump (see also the comment at the end of Section 2). It is therefore entirely consistent for such modulations to be stationary even though a single lump propagating along the characteristic curves experiences changes over time. In this respect, the situation is perfectly analogous to what happens with the reductions of the soliton Whitham modulation equations studied in [Reference Ryskamp, Hoefer and Biondini46, Reference Ryskamp, Maiden, Hoefer and Biondini47], where it was demonstrated that, even though the modulation equations are independent of one of the coordinate variables, the corresponding dynamics for the soliton still involves changes over that coordinate.
4. Interactions between lumps and a rarefaction wave
The 
$Y$-independent reduction of the dKP equation is the familiar Hopf equation, 
$\bar{u}_t + \bar{u} \bar{u}_x = 0$, which admits the rarefaction wave solution
\begin{equation}
\bar{u}(X,T)=\begin{cases}
0,\quad&X \lt 0\\
X/T,\quad&0 \lt X \lt T\\
1,\quad&X \gt T
\end{cases},
\end{equation}corresponding to the initial condition 
$\bar{u}(X,Y,0)=\vartheta(X)$ for the dKP equation, where 
$\vartheta(X)$ is the Heaviside step function. Thus, (45) is also a 
$Y$-independent solution of the dKP equation (16). The bottom value can be taken to be zero without loss of generality due to the Galilean invariance of the KP equation, and the height of the step can be taken to be 
$1$ without loss of generality due to the scaling invariance of the KP equation. We should note that the case of an initial step placed along a line 
$X = c Y$ with 
$c\ne 0$ instead of 
$X=0$ (corresponding to an initial condition 
$u(X,Y,0) = \vartheta(X - c Y)$) would result in a slanted rarefaction wave. On the other hand, this scenario can always be reduced to the above scenario using the invariance of the KP equation under a slanted change of the coordinate axes (e.g., see [Reference Ablowitz, Biondini and Wang1]). We can, therefore, limit ourselves to considering 
$c =0$ without loss of generality.
If a lump is initially placed at a location 
$x_{0} \lt 0$ to the left of the step, its horizontal velocity is positive by (2c). Thus, the lump will eventually reach the ‘ramp’ region of the rarefaction wave. In this section, we are interested in studying the dynamics of such a lump as it travels up the ramp. We do so by making use of the 
$XT$-dependent modulation system (29). An illustration of this scenario is provided in Fig. 1, in which the left and right columns show temporal slices corresponding respectively to the cases of a horizontal and non-horizontal lump velocity. In both cases chronicled in Fig. 1, the lump propagates into the expanding rarefaction wave, interacts with it for some time and exits the rarefaction with altered parameters (velocity, amplitude). On the other hand, the rarefaction wave is not visibly changed by the lump post-interaction. (See the appendix for a description of the numerical methods used.) Next, we consider these two cases separately, since they give rise to qualitatively different results.
Density plots of the numerically computed solution 
$u(x,y,t)$ of the KPI equation as a function of 
$x$ and 
$y$ at the times 
$t = 0$, 
$t = 20$, 
$t = 40$, 
$t = 60$ and 
$t=80$ (top to bottom) with 
$a_0=1.75$ and 
$q_0=0$ (left) and 
$q_0=0.1$ (right). The rarefaction wave region lies between the two dashed vertical lines. (See Appendix for details of the numerical methods.).
1 Long description
The composite image contains eight heatmaps arranged in two columns and four rows. Each heatmap shares the same axis structure. The horizontal axis spans from 0 to 250 and the vertical axis spans from 0 to 2. A vertical dashed line appears near the horizontal position of 100 in each heatmap, dividing the grid into a left portion and a right portion. Each heatmap has a gradient scale bar positioned to its right, with values ranging from 0 to 2, where lower values correspond to a lighter appearance and higher values correspond to a more intense concentrated appearance. Across all eight heatmaps, the field values are predominantly low and uniform throughout most of the grid area. A small localized region of elevated intensity is visible in each heatmap, appearing as a distinct concentrated spot. In the top row, the localized spot appears near the horizontal position of 100 to 110 and is positioned near the vertical center of the grid. Moving down through subsequent rows, the position of the localized high-intensity spot shifts progressively toward higher horizontal axis values, reaching near the horizontal position of 200 to 225 in the bottom row. The left column and right column heatmaps at each row level show closely similar spatial distributions, with the localized spot appearing at comparable horizontal and vertical positions within each paired row. The surrounding field values remain near 0 throughout, making the localized elevated region the primary visible feature in each heatmap.
4.1. Interactions between a horizontally travelling lump and a rarefaction wave
We begin with the simpler scenario 
$q=0$, in which case the velocity of the lump is oriented purely in the horizontal direction. In particular, (2c) yields the initial horizontal velocity of the lump as 
$a_{0}^{2}$, where 
$a_{0}$ is the initial value of the amplitude parameter 
$a$. If the lump is transmitted through the rarefaction wave, it will emerge after some time from the ramp with an amplitude 
$a_1$ to be determined and propagate (on top of the unit mean) with a constant velocity of 
$a_{1}^2+1$.
To determine 
$a_{1}$ in terms of 
$a_{0}$, one can use the Riemann invariant of the 
$XT$-dependent modulation system 
$R_{2}=a^2/2+\bar{u}$. In particular, from (34), we know that this Riemann invariant is constant along the characteristics defined by 
$\mathrm{d} X/\mathrm{d} T=\lambda_{2}$ where 
$\lambda_{2}=\bar{u}+a^2$ is exactly the lump velocity. That is, 
$R_{2}$ is constant on the trajectory of the lump. As such, the values of 
$R_{2}$ before and after the interaction with the rarefaction ramp must be equal, which implies that
\begin{equation}
a_{1}=\sqrt{a_{0}^{2}-2}.
\end{equation} Note that this relationship suggests that if 
$a_{0} \lt \sqrt{2}$, the lump will never emerge from the ramp. We therefore have the following dichotomy:
(1)
$a_{0} \gt \sqrt{2}$ results in lump transmission,(2)
$a_{0} \lt \sqrt{2}$ results in lump trapping.
This result is similar to the transmission condition found for the KdV soliton-rarefaction wave interaction in [Reference Ablowitz, Cole, El, Hoefer and Luo3]. As we show next, further information about the dynamics of the lump within the ramp region, including the precise temporal dependence of the lump amplitude as it travels through the rarefaction wave, can be obtained by explicitly solving the modulation system (25) or, equivalently in this case, the reduction (28a).
In the 
$XT$-dependent reduction with 
$q=0$, (28b) is satisfied trivially. Moreover, once a solution for (28c) is given, (28a) yields, along characteristics,
\begin{equation}
\frac{\mathrm{d} a}{\mathrm{d} T}=-a\bar{u}_{X},\qquad\frac{\mathrm{d} X}{\mathrm{d} T}=\bar{u}+a^2.
\end{equation} Since the mean field at the bottom of the ramp is zero and the lump is initially centred at the location 
$x_0 \lt 0$, (2c) yields the time at which the lump reaches the base of the ramp as
\begin{equation}
T_{0}=-\frac{x_{0}}{a_{0}^2}.
\end{equation} In the ramp region, we can substitute the self-similar solution 
$\bar{u}=X/T$ from (45), in which case the solution of (47) subject to the initial conditions 
$a(T_{0})=a_{0}$, 
$X(T_{0})=0$ is given by
\begin{equation}
a(T)=\frac{a_{0}T_{0}}{T},\qquad X(T)=\frac{a_{0}^{2}}{2}T\left(1-\frac{T_{0}^{2}}{T^{2}}\right).
\end{equation} Now, if the lump is to transmit through the rarefaction wave, there must exist a time 
$T_{1} \gt T_{0}$ such that 
$X(T_{1})=T_{1}$ is satisfied, which yields a condition that can be used to solve for 
$T_1$. With the form of 
$X(T)$ given above, the only positive solution of this equation is
\begin{equation}
T_{1}=\frac{a_{0}T_{0}}{\sqrt{a_{0}^{2}-2}}.
\end{equation} Note that 
$T_{1}$ is not real-valued when 
$a_{0} \lt \sqrt{2}$, consistent with the previously found trapping condition. Conversely, if 
$T_{1}$ is real-valued, then the final value of the amplitude parameter after the lump emerges from the ramp is
\begin{equation}
a_{1}=a(T_{1})=\frac{a_{0}T_{0}}{T_{1}},
\end{equation}which coincides with (46) after substituting (50).
Recall that the physical amplitude 
$A$ of the lump is related to the parameter 
$a$ via 
$A=8a^2$. Summarizing the present results in terms of the amplitude, we then have that if the initial amplitude of the lump satisfies the condition 
$A_{0} \gt 16$, then the lump is transmitted through the rarefaction wave with a final amplitude of 
$A_{1}=A_{0}-16$ on top of the unit mean field. Conversely, in the case where 
$A_{0} \lt 16$, the lump will remain trapped in the ramp forever, never reaching the top of the rarefaction wave, and its amplitude will asymptotically decay to zero as 
$T\to\infty$, similarly to what happens with soliton trapping in the KdV equation [Reference Ablowitz, Cole, El, Hoefer and Luo3].
Note that the full trajectory of a transmitting lump for all time is described by
\begin{equation}
X_{\text{lump}}(T)=\begin{cases}
x_{0}+a_{0}^{2}T,\quad&T\leq T_{0}\\
\frac{a_{0}^{2}}{2}T\left(1-\frac{T_{0}^{2}}{T^{2}}\right),\quad& T_{0} \lt T \lt T_{1}\\
x_{1}+(a_{1}^{2}+1)T,\quad&T\geq T_{1}
\end{cases}.
\end{equation} Using the explicit expressions (49) shows that 
$x_{1}=x_{0}a_{1}/a_{0}$, i.e., the lump emerges from the interaction with a phase shift given by
\begin{equation}
\delta=x_{1}-x_{0}=\left(\frac{a_{1}}{a_{0}}-1\right)x_{0}.
\end{equation} This phase shift can also be determined without knowledge of the full trajectory. To this end, we introduce an equally spaced array of horizontally travelling lumps with wavenumber 
$k$, such that 
$k/(2\pi)\ll1$, and couple an approximate conservation of waves equation expressing the conservation of the number of lumps in the train to the modulation system,
\begin{equation}
k_{T}+\big[\big(2R_{2}-\bar{u}\big)k\big]_{X}=0,
\end{equation}where 
$2R_{2}-\bar{u}=a^2+\bar{u}$ is the lump velocity. We restrict to the trajectory of the lumps, on which 
$R_{2}$ is constant. In this case, we find that (54) is diagonalized by the Riemann invariant 
$R_{k}$ given by
\begin{equation}
\partial_{T}R_{k}+(2R_{2}-\bar{u})\partial_{X}R_{k}=0,\qquad R_{k}=k\sqrt{R_{2}-\bar{u}}.
\end{equation}The constancy of this Riemann invariant along the lump trajectory yields
\begin{equation}
k_{1}\sqrt{R_{2}-1}=k_{0}\sqrt{R_{2}},
\end{equation}where 
$k_{0}$ and 
$k_{1}$ are the wavenumber (i.e., lump spacing) before and after the interaction, respectively. The phase shift can then be predicted by examining the ratio of the spacing before and after the interaction in the limit as 
$k\rightarrow0$,
\begin{equation}
\frac{x_{1}}{x_{0}}\sim\frac{k_{0}}{k_{1}}\sim\sqrt{\frac{R_{2}-1}{R_{2}}}=\frac{a_{1}}{a_{0}},
\end{equation}which is consistent with (53).
Figure 2 displays a comparison between the analytical prediction (49) for the amplitude parameter and the numerically measured amplitude of the lump throughout the simulation. In particular, the left panel corresponds to a case of lump transmission (i.e., 
$a_0 \gt \sqrt2$), while the right panel illustrates the phenomenon of lump trapping (i.e., 
$a_0 \lt \sqrt2$), in which the lump never emerges from the ramp and its amplitude decays for all time, and the centre panel corresponds to the critical case (i.e., 
$a_0 = \sqrt2$). Here and in subsequent comparisons with numerics, we set 
$\varepsilon=1$ so that 
$T=t$, 
$X=x$, 
$Y=y$ and set 
$-x_0 \gg 1$. Note how the analytical prediction for the final lump amplitude is in excellent agreement with the results from the numerical simulations. (See Fig. 5 for further confirmation of this agreement.) The only small disagreements between the theory and the numerics are the small oscillations of the lump amplitude before it enters the rarefaction wave and the smooth transition to the asymptotic value of the lump amplitude. Both of these phenomena are due to dispersive effects that are not captured by the ramp solution (45), which results from the reduction from the dispersive KP equation to the dispersionless Hopf equation. On the other hand, the general features of the temporal dependence of the lump amplitude, and its asymptotic value are very well captured by the solution of the modulation equations. (The smaller noise-like oscillations of the lump amplitude towards the end of simulation are likely due to corruption of the lump due to interference with the rest of the numerical solution field. See the Appendix for a discussion of the numerical methods.)
Numerically measured lump amplitude parameter, 
$a(t)$ as a function of 
$t$ (solid blue line) with theoretical prediction (dashed black line) in a case of lump transmission, 
$a_0=1.75$, 
$q_0 = 0$ (left), critical trapping threshold 
$a_0 = \sqrt{2}$, 
$q=0$ (centre) and below the critical threshold 
$a_0 = 1.25$, 
$q=0$ (right).
Figure 2 Long description
Null
4.2. Interactions between an obliquely travelling lump and a rarefaction wave
We now move on to study the interaction between the rarefaction wave and a full two-parameter lump with nonzero transverse velocity. In this case, from (2c), we have that if the lump is initially placed at a location 
$(x_{0},y_{0})$ with 
$x_0 \lt 0$ and 
$y_{0}$ arbitrary, the time at which the lump reaches the base of the rarefaction ramp is
\begin{equation}
T_{0}=-\frac{x_{0}}{a_{0}^2+q_{0}^2}.
\end{equation} Interestingly, we find that in this case, the lump will always transmit through the ramp as long as 
$q_{0}\neq0$. We first show that this is the case by considering the Riemann invariants of the 
$XT$-dependent modulation system (29) given by 
$R_{1}=aq$ and 
$R_{2}=(a^2-q^2)/2+\bar{u}$. Using similar logic as in the previous subsection, since 
$R_{1}$ and 
$R_{2}$ are both constant on the trajectory of the lump, it must be the case that
\begin{equation}
a_{1}q_{1}=a_{0}q_{0},
\end{equation}
\begin{equation}
a_{1}^2-q_{1}^2=a_{0}^2-q_{0}^{2}-2.
\end{equation}
$q_{1}=a_{0}q_{0}/a_{1}$ in (59b) to get the biquadratic equation
\begin{equation}
a_{1}^4-\mu a_{1}^2-a_{0}^2q_{0}^2=0,\qquad\mu=a_{0}^2-q_{0}^2-2,
\end{equation}whose only positive solution is
\begin{equation}
a_{1}^2=\frac{1}{2}\left(\mu+\sqrt{\mu^2+4a_{0}^2q_{0}^2}\right).
\end{equation} Note that in the one-parameter case 
$q_{0}=0$, the sign of 
$\mu$ distinguishes between the cases of trapping and transmission. However, in the two-parameter case 
$q_{0}\neq0$, the final amplitude parameter 
$a_{1}$ is well-defined regardless of the sign of 
$\mu$.
The full trajectory of the lump as it interacts with the rarefaction wave can be found by solving (25). The evolution of the amplitude parameter 
$a$ is the same as in the previous subsection, since the first equation in (25) applies in either case. Moreover, since 
$aq$ is a constant of motion for a 
$Y$-independent mean field, we have that 
$q$ grows with time as 
$a$ decays. That is, explicitly,
\begin{equation}
a(T)=\frac{a_{0}T_{0}}{T},\qquad q(T)=\frac{q_{0}T}{T_{0}}.
\end{equation} The time-independence of the product 
$aq$ throughout the interaction is confirmed by the numerical results, as shown in Fig. 3. The temporal dependence of 
$q(T)$ is also shown in Fig. 3 in two representative cases and compared to the theoretical predictions from (62), again showing very good agreement. (As before, we suspect that the small noise-like oscillations of 
$q$ towards the end of simulation are due to corruption of the lump due to interference with the rest of the numerical solution field.) Note that the growth of 
$q$ is the reason that a lump with nonzero transverse momentum is always transmitted through a rarefaction wave: Since the horizontal component of the lump velocity is 
$\bar{u}+a^2+q^2$, the fact that 
$q$ grows with time shows that the lump experiences acceleration in the horizontal direction as a result of the interaction.
The lump parameter 
$q(t)$ as a function of time as the lump propagates through the rarefaction wave, with 
$a_0=1.75$ and with 
$q_0=1$ (left panel) and 
$q_0=0.5$ (centre panel) as reconstructed from numerical simulations (blue and gold curves), compared to the theoretical predictions (dashed black curves). Right panel: The product 
$aq$ as function of time in three simulations with 
$a_0=1.75$ and 
$q_0=0$ (green curve), 
$q_0=0.5$ (gold curve) and 
$q_0=1$ (blue curve).
Figure 3 Long description
Null
Finally, we can solve the modulation equations for the characteristic trajectories to find
\begin{equation}
X(T)=\frac{a_{0}^{2}}{2}T\left(1-\frac{T_{0}^{2}}{T^{2}}\right)-\frac{q_{0}^{2}}{2}T\left(1-\frac{T^{2}}{T_{0}^{2}}\right),\qquad Y(T)=q_{0}T_{0}\left(1-\frac{T^{2}}{T_{0}^{2}}\right)+Y(T_{0}).
\end{equation} The time at which the lump exits the ramp is 
$T_{1}$ such that 
$X(T_{1})=T_{1}$. This is equivalent to
\begin{equation}
q_{0}^{2}T_{1}^{4}+\mu T_{0}^{2}T_{1}^2-a_{0}^{2}T_{0}^{4}=0.
\end{equation}The relevant solution is
\begin{equation}
T_{1}=\frac{T_{0}}{\sqrt{2}|q_{0}|}\sqrt{-\mu+\sqrt{\mu^2+4a_{0}^{2}q_{0}^2}},
\end{equation}which is defined for either sign of 
$\mu$ as long as 
$q_{0}\neq0$. Note that taking the limit of the above expression as 
$q_{0}\rightarrow0$ recovers (50). With this expression for 
$T_{1}$, 
$a_{1}=a(T_{1})$ is identical to the previously obtained prediction (61). Figure 4 shows the theoretical predictions in (63) for the lump coordinates in a number of representative cases, compared with the results of corresponding numerical simulations. The agreement between theory and numerical results is excellent in all cases considered.
The 
$x$-coordinate of the lump (left panel), the 
$y$-coordinate of the lump (centre panel) and the lump trajectory 
$(x(t),y(t))$ (right panel) for 
$a=1.75$ and 
$q=0$ (green), 
$q=0.5$ (gold) and 
$q=1$ (blue) compared with theory (dashed black line). The squares mark the time at which the lump enters the rarefaction wave region, while the circles mark the time at which the lump exits the rarefaction wave.
Figure 4 Long description
A multi-line line graph with the vertical axis label, x of t and the horizontal axis label, t. The horizontal axis tick labels are 0, 15, 30. The vertical axis tick labels are 0, 40, 80. Three straight lines start at (0, 0). Series 1 uses square and circle markers: (12, 50), (15, 65). Series 2 uses square and circle markers: (15, 50), (22, 70). Series 3 uses square and circle markers: (16, 50), (28, 78). Trend summary for image A: all three series increase with t; Series 1 is steepest among the listed points, Series 3 is least steep among the listed points. The image B showing a multi-line line graph with the vertical axis label, y of t and the horizontal axis label, t. The horizontal axis tick labels are 0, 15, 30. The vertical axis tick labels are 0, 40, 80. Three straight lines start at (0, 0). Series 1 uses square and circle markers: (10, 25), (15, 35). Series 2 uses square and circle markers: (15, 15), (22, 25). Series 3 lies on y equals 0 with markers: (15, 0), (28, 0). Trend summary for image B: Series 1 increases faster than Series 2 across the listed points; Series 3 remains at y equals 0. The image C showing a multi-line line graph with the vertical axis label, y of t and the horizontal axis label, x of t. The horizontal axis tick labels are 0, 50, 100. The vertical axis tick labels are negative 10, 0, 20, 40. Three straight lines start at (0, 0). Series 1 uses square and circle markers: (50, 25), (60, 32). Series 2 uses square and circle markers: (50, 15), (70, 25). Series 3 lies on y equals 0 with markers: (50, 0), (80, 0). Trend summary for image C: Series 1 increases faster than Series 2 across the listed points; Series 3 remains at y equals 0.
Summarizing, in terms of the physical lump amplitude 
$A=8a^{2}$, we have that any lump with initial amplitude 
$A_{0}$ and 
$q_{0}\neq0$ is expected to emerge from the rarefaction wave with new amplitude given by
\begin{equation}
A_{1}=4\left(\mu+\sqrt{\mu^2+\frac{1}{2}A_{0}q_{0}^2}\right),\qquad\mu=\frac{1}{8}A_{0}-q_{0}^2-2.
\end{equation}The trajectory of the lump for all time is
\begin{equation}
X_{\text{lump}}(T)=\begin{cases}
x_{0}+(a_{0}^{2}+q_{0}^{2})T,\quad&T\leq T_{0}\\
\frac{a_{0}^{2}}{2}T\left(1-\frac{T_{0}^{2}}{T^{2}}\right)-\frac{q_{0}^{2}}{2}T\left(1-\frac{T^{2}}{T_{0}^{2}}\right),\quad& T_{0} \lt T \lt T_{1}\\
x_{1}+(a_{1}^{2}+q_{1}^{2}+1)T,\quad&T\geq T_{1}
\end{cases},
\end{equation}
\begin{equation}
Y_{\text{lump}}(T)=\begin{cases}
y_{0}-2q_{0}T,\quad&T\leq T_{0}\\
q_{0}T_{0}\left(1-\frac{T^{2}}{T_{0}^{2}}\right)+y_{0}-2q_{0}T_{0},\quad& T_{0} \lt T \lt T_{1}\\
y_{1}-2q_{1}T,\quad&T\geq T_{1}
\end{cases}.
\end{equation}
$T=T_{1}$ we find
\begin{equation}
x_{1}=\frac{x_{0}T_{1}}{T_{0}}\frac{a_{1}^2+q_{1}^2}{a_{0}^{2}+q_{0}^{2}},\qquad y_{1}=y_{0}-q_{0}T_{0}\left(1-\frac{T_{1}^{2}}{T_{0}^{2}}\right).
\end{equation} This first equation above shows that the phase shift in 
$x$ as a result of the interaction is
\begin{equation}
\delta=x_{1}-x_{0}=\left(\frac{T_{1}}{T_{0}}\frac{a_{1}^{2}+q_{1}^{2}}{a_{0}^{2}+q_{0}^{2}}-1\right)x_{0},
\end{equation}which reduces to (53) when 
$q_{0}=0$. Similarly to the previous case, this phase shift can also be obtained by considering the conservation of waves type equation
\begin{equation}
k_{T}+\left[\left(\bar{u}+2\sqrt{(R_{2}-\bar{u})^{2}+R_{1}^{2}}\right)k\right]_{X}=0,
\end{equation}with wavenumber 
$k/(2\pi)\ll1$. Along the lump trajectory, 
$R_{1}$ and 
$R_{2}$ are both constant and (70) is diagonalized by the Riemann invariant
\begin{equation}
R_{k}=kp(\bar{u}),\qquad p(\bar{u})=\sqrt{\frac{(R_{2}-\bar{u})^2+R_{1}^2}{R_{2}-\bar{u}+\sqrt{(R_{2}-\bar{u})^2+R_{1}^2}}}.
\end{equation} Using the fact that 
$R_{k}$ remains constant throughout the interaction yields
\begin{equation}
\frac{x_{1}}{x_{0}}\sim\frac{k_{0}}{k_{1}}\sim\frac{p(1)}{p(0)}.
\end{equation}Direct calculations show that this provides the same phase shift as in (69).
In Fig. 5, the predicted values of the parameters of the transmitted lump are plotted as a function of the amplitude of the initial lump for several values of 
$q_{0}$. The theoretical predictions are then compared with the numerically measured values from the corresponding numerical simulations. In all cases, the agreement between theoretical predictions and numerical results is excellent.
Left: Analytical prediction for the final lump amplitude 
$A_\mathrm{out} = 8a_1^2$ as a function of the initial lump amplitude 
$A_\mathrm{in} = 8a_0^2$ for 
$q_0 = 0$ (blue curve), 
$q_0=0.5$ (gold curve) and 
$q_0=1$ (green curve) versus the numerically recovered values, marked by circles and triangles. Right panel: Analytical prediction for the final slope parameter 
$q_1$ as a function of the initial lump amplitude 
$A_\mathrm{in} = 8a_0^2$ (for the same values of 
$q_0$ as in the left panel) versus the numerically recovered values.
Figure 5 Long description
Null
5. Discussion and outlook
In this work, we derived a four-component (2+1)-dimensional hyperbolic system of equations governing the modulations of a lump solution of the KPI equation in the presence of a slowly varying mean field. We showed that the system passes the Haantjes tensor test for integrability, admits a Hamiltonian formulation when the mean field is stationary, and we diagonalized it in terms of Riemann invariants in each of its three 2D reductions. As a concrete application, we used the modulation system to derive sharp transmission and trapping criteria for a lump traversing a mean-field rarefaction wave: a horizontally travelling lump is transmitted if and only if its initial amplitude 
$A_{0}$ exceeds the critical value 
$16$, while a lump with any nonzero transverse parameter 
$q_0$ is always transmitted, with the transverse momentum 
$aq$ conserved throughout the interaction. The analytic predictions were compared with direct numerical simulations of the KPI equation and found to be in excellent agreement.
Although this work was motivated by recent extensions of Whitham modulation theory to the line-soliton and one-phase periodic solutions of the KP equation [Reference Ablowitz, Biondini and Wang1], it differs from [Reference Ablowitz, Biondini and Wang1] in one essential respect. Line solitons and one-phase periodic solutions are co-dimension-one objects (i.e., the solution is constant along the wave front), and the corresponding modulation theory is built by averaging over a one-dimensional (1D) family of phases. Lumps, by contrast, are localized in both spatial directions and decay only algebraically at infinity, so the conventional phase average is unavailable. The derivation in Section 2 instead proceeds from integral identities for the lump itself, taken in the moving reference frame of the lump and at leading order in the modulation parameter 
$\varepsilon$. To the best of our knowledge, the system (19a and 19b) is therefore the first Whitham-type modulation theory for a genuinely 2D, fully localized coherent structure of an integrable PDE.
The trapping/transmission dichotomy obtained in Section 4 places this work within the broader research on soliton–mean field interactions that has emerged over the past decade in the context of dispersive hydrodynamics [Reference El and Hoefer22]. For the KdV equation, the analogous problem of a soliton traversing a rarefaction was analysed in [Reference Ablowitz, Cole, El, Hoefer and Luo3], and the trapping condition 
$A_0 \lt 2$ was obtained there; closely related analyses of soliton–dispersive shock and soliton–periodic mean interactions have since been carried out for KdV and other weakly dispersive hydrodynamic models in [Reference Biondini and Lottes14, Reference El and Hoefer22, Reference Maiden, Anderson, Franco, El and Hoefer39]. Indeed, the extent to which the modulation theory of a genuinely 2D integrable coherent structure parallels the 1D KdV theory [Reference Ablowitz, Cole, El, Hoefer and Luo3] is perhaps striking: a Hamiltonian system on the parameter space, Riemann invariants in each coordinate reduction and a sharp threshold separating transmission from trapping. At the same time, our results extend this programme in an essential way: here the coherent structure is genuinely 2D, so that nontrivial transverse dynamics (notably the conservation of 
$aq$ that drives the unconditional transmission of obliquely travelling lumps) becomes available for the first time.
We also emphasize that the quantitative predictions for the final lump amplitude are in principle testable in the same optical setting in which lumps were recently observed [Reference Dieli, Pierangeli, Baronio, Trillo and Conti20], and we hope the present work will stimulate such experiments.
The fact that the modulation system (19a and 19b) passes the Haantjes test, admits a Hamiltonian form for stationary means and reduces to diagonal (1+1)-systems with explicit Riemann invariants in all three coordinate-plane reductions strongly suggests that the system should be amenable to the full machinery of Tsarev’s generalized hodograph method [Reference Tsarev51] and to the broader theory of integrable systems of hydrodynamic type [Reference Dubrovin and Novikov21, Reference Ferapontov and Khusnutdinova24, Reference Pavlov43]. Two questions arise naturally. First, can the full 
$(2{+}1)$-system be integrated by hydrodynamic-type techniques – for example, by exhibiting four functionally independent commuting flows, or by embedding its reductions in a hydrodynamic chain in the sense of Pavlov? Second, all three (1+1)-dimensional reductions analysed in Section 3 fail to be strictly hyperbolic, with two characteristic speeds always coinciding. The geometric meaning of this degeneracy and its consequences for shock formation in the modulation variables deserve a more thorough analysis. (We note that analogous degeneracies appear in the line-soliton modulations of the KPII equation [Reference Ablowitz, Biondini and Wang1, Reference Ryskamp, Hoefer and Biondini46].)
A complementary perspective is provided by Krichever’s algebro-geometric approach to averaging on integrable systems [Reference Krichever34]. Since lumps arise as further rational degenerations of singular-curve theta-function solutions of the KP equation [Reference Agostini, Celik and Little9,Reference Little37], the modulation system (19a and 19b) may admit an interpretation as a singular limit of the corresponding Krichever–Whitham equations for the periodic case. Establishing this connection rigorously, and using it to inform the modulations of higher-order rational solutions, is an attractive open problem. At the same time, we should point out that the same problem is also still open for the one-phase KP–Whitham system derived in [Reference Ablowitz, Biondini and Wang1].
Finally, a long-standing obstacle in the analysis of the KP equation is the absence of an effective inverse scattering theory for genuinely non-decaying initial data: the IST has been developed for decaying data [Reference Ablowitz and Fokas4, Reference Ablowitz and Fokas5, Reference Ablowitz and Segur7] and for perturbations of an exact lump or line soliton [Reference Boiti, Pempinelli and Pogrebkov15, Reference Boiti, Pempinelli, Pogrebkov and Prinari16, Reference Villarroel and Ablowitz52, Reference Wu54], but not for backgrounds of the type considered here. We believe that the development of a viable IST for non-decaying initial conditions and the formulation of a suitable semiclassical limit that can derive the system (19a and 19b) in a rigorous way is one of the most challenging and at the same time intriguing open problems suggested by this work.
We end the discussion by mentioning that we view this work as the initiation of a new research programme dealing with modulation theory for 2D localized structures. In addition to the open problems discussed above, there are many directions for further investigation, including modulations of multi-lump configurations and lump chains, lump interactions with more complex mean fields, such as dispersive shock waves, and modulation theory for lump-like solutions of other physically relevant models. We hope that the results of this work and the present discussion will stimulate further work on these topics.
Acknowledgements
We thank Gennady El for many interesting discussions on related topics.
Funding statement
G.B. was partially supported by the Simons Foundation under grant number SFI-MPS-TSM-00013369. The work of M.H. was supported by NSF Grant No. DMS-2306319.
Competing interests
The authors declare that they have no conflict of interest.
Appendix
In this appendix, we provide a concise description of the numerical methods used to perform the numerical simulations discussed in the main text.
Numerical integration of the KP equation. The KP equation (1) is solved numerically by a pseudo-spectral method combined with a fourth-order exponential time-differencing in time [Reference Kassam and Trefethen32] by recasting the Fourier-transformed equation into the form
\begin{equation}
\partial_t \hat u_k = L_k \hat u_k + \hat N(u)
\end{equation}where 
$\hat u_k(t)$ is the Fourier coefficient of 
$u(x,y,t)$ with 
$k = (k_x,k_y)$, with
\begin{equation}
L_k =
\begin{cases} \displaystyle
ik_x^3 + i{k_y^2}/{k_x}\,\,&k_x\neq 0 \\
0\,\, &\mbox{otherwise}
\end{cases}
\end{equation}and where the non-linear term is
\begin{equation}
\hat N(u) = -ik_x \mathcal{F} \left[\frac{u^2}{2}\right],
\end{equation} 
$\mathcal{F}$ denoting the Fourier transform. Equation (A.1) is recast as an integral equation,
\begin{equation}{\widehat u}_k(t+h)=e^{hL_k}{\widehat u}_k+\int_0^he^{L(h-s)}\widehat N(u(t+s))\,\mathrm ds,\end{equation}and the non-linear term is approximated by a fourth order Runge-Kutta (RK4) method. The resulting method is often referred to as the ETDRK4 [Reference Kassam and Trefethen32]. The accuracy of the time integration is verified by monitoring the four conserved integrals
\begin{equation}
M = \iint_{\mathbb{R}^{2}} u \,\mathrm{d} x\,\mathrm{d} y, \qquad
P = \iint_{\mathbb{R}^{2}} u^2 \,\mathrm{d} x\,\mathrm{d} y, \qquad
Q = \iint_{\mathbb{R}^{2}} u v \,\mathrm{d} x\,\mathrm{d} y,
\end{equation}
\begin{equation}
E = \iint_{\mathbb{R}^{2}} \left[u_x^2 - v^2 -\frac{1}{3}u^3\right] \,\mathrm{d} x\,\mathrm{d} y,
\end{equation}
$v_x = u_y$. Over the course of all simulations reported, the largest numerical error was observed in 
$E$, with 
$\mathop{\rm max}\nolimits(\Delta E/E) \sim 10^{-4}$.
Initial and boundary conditions, simulation parameters. The theoretical problem is posed on the infinite plane 
$(x,y)\in\mathbb{R}^2$. For the purposes of numerical simulation, we restrict the spatial variables to a computational box 
$(x,y)\in[0,L_x]\times[0,L_y]$ with periodic boundary conditions in both directions. Note that the Heaviside initial condition considered in Section 4 is not periodic with respect to 
$X$. Moreover, the discontinuous transition between the two step values would result in the generation of spurious high-wavenumber components that would pollute the solution field. To eliminate both of these problems, we choose numerical initial conditions consisting of a symmetrized, regularized step function,
\begin{equation}
u_s(x;x_s,w,\delta) = \tfrac{1}{2}\left(1 - \tanh \left(f(x)/\delta\right) \right),
\qquad
f(x) = (x-x_s)^2 - w^2,
\end{equation}where 
$x_s$ and 
$w$ control the step centre and its half-width, and 
$\delta$ gives the characteristic width of the step transition smoothness, a discontinuous step corresponding to the limit 
$\delta\to0$. The evenness of 
$f(\cdot)$ with respect to 
$x_s$ ensures that the limits of (A.6) as 
$x\to\pm\infty$ coincide. A lump soliton 
$U(\cdot)$ with parameters 
$a$ and 
$q$ given by equations (2a), (2b) and (2c) is then added to the above step-like initial condition, so that the full initial condition is
\begin{equation}
u(x,y,t=0) = u_s(x; x_s, w, \delta) + U(x-x_l,y-y_l; a_o,q_o),
\end{equation}where the initial lump parameters are 
$x_l = 50$, 
$y_l = 93.75$ and the step parameters are 
$x_s=275$, 
$w = 175$ and 
$\delta=100$ on a computational domain up to 
$L_x = 500$ by 
$L_y = 125$. (Note that the slow, algebraic decay of the lump solution leads to the need to use a large spatial domain.) The initial lump parameters 
$a_0$ and 
$q_0$ differ among the various simulations (see Figs 2–5). Up to 
$32768\times8192$ Fourier modes were used to resolve the solution on the numerical grid, even though in some cases a much smaller number of Fourier modes was sufficient to accurately capture the solution. The computational were performed on the computers of the Center for Computational Research at the University at Buffalo.
Lump detection algorithm. In the course of the simulations, the lump propagates into the rarefaction wave and needs to be detected, and its parameters 
$a,q$ recovered. A postprocessing script is run on the numerical data for 
$u(x,y,t_n)$ for a set of time slices 
$t_n = n\tau$ where 
$\tau = T_{max}/M$ with 
$M=500$ and 
$T_\mathrm{max} = 100$. At each time slice, a rough estimate for the 
$\mathop{\rm max}\nolimits u(x,y,t_n)$ is determined by scanning the maximum on the numerical grid,
\begin{equation}
(x_{i_\mathrm{max}}, y_{j_\mathrm{max}}) = \arg \mathop{\rm max}\nolimits \limits_{(x_i,y_j)} u(x_i,y_j,t_n).
\end{equation}Next, a local grid patch at 
$(x_i,y_j)$ with 
$i = i_\mathrm{max}-l,\ldots,i_\mathrm{max}+l$ and 
$j = j_\mathrm{max}-l,\ldots,j_\mathrm{max}+l$ is considered, with the parameter 
$l$ controlling the size of this local patch. Once the data points 
$(x_i,y_j,u_{ij})$ are determined, a non-linear least-squares fit to the lump soliton (2a, 2b and 2c) is performed to determine five parameters, 
$(x_\mathrm{max},y_\mathrm{max})$, 
$a$,
$q$ and 
$\bar u$. The resulting lump parameters are then compared to the theoretical predictions in Figs 2–5.
Open access
