# Relaxation and fluctuations of a mass- and dipole-conserving stochastic
lattice gas

Baruch Meerson<EMAIL_ADDRESS>Racah Institute of Physics, Hebrew
University of Jerusalem, Jerusalem 91904, Israel

###### Abstract

Han et al. [Phys. Rev. Lett. 132, 137102 (2024)] have recently introduced a
classical stochastic lattice gas model which, in addition to particle
conservation, also conserves the particles’ dipole moment. Because of its
intrinsic nonlinearity this model exhibits unusual macroscopic scaling
behaviors, different from those of lattice gases that conserve only the number
of particles. Here we investigate some basic relaxation and fluctuation
properties of this model at large scales and at long times. These properties
crucially depend on whether the total number of particles is infinite or
finite. We find similarity solutions, describing relaxation of the dipole-
conserving gas (DCG) in several standard settings. A major part of our effort
is an extension to this model of the macroscopic fluctuation theory (MFT),
previously developed for lattice gases where only the number of particles is
conserved. We apply the MFT to the calculation of the variance of
nonequilibrium fluctuations of the excess number of particles on the positive
semi-axis when starting from an (either deterministic, or random) constant
density at $t=0$. Using the MFT, we also identify the equilibrium Boltzmann-
Gibbs distribution for the DCG. Finally, based on these results, we determine
the probability distribution of, and the most probable density history leading
to, a large deviation in the form of a macroscopic void of a given size in an
initially uniform DCG at equilibrium.

##  I Introduction

Lattice gases, composed of particles undergoing symmetric random hopping,
exhibit diffusive behavior at large scales and long times Spohn91 ; Anna ;
Kipnis99 ; Schutz00 . As these models conserve the number of particles
locally, the gas density $\rho(\mathbf{x},t)$ at large scales and at long
times obeys the continuity equation

$\partial_{t}\rho+\nabla\cdot{\mathbf{J}}=0$ (1)

with the diffusion current

${\mathbf{J}}=-D(\rho)\nabla\rho\,.$ (2)

For a limited class of models – generically for the so called gradient lattice
gases Spohn91 ; Kipnis99 ; Gabrielli18  – the diffusion coefficient $D(\rho)$
can be calculated from the microscopic model exactly. In simple cases, which
include non-interacting random walkers, $D(\rho)$ is independent of the
density, and Eqs. (1) and (2) reduce to the simple (linear) diffusion equation

$\partial_{t}\rho=D\Delta\rho\,,$ (3)

where $\Delta$ is the Laplacian operator.

How is the simple diffusion model modified if the lattice gas is constrained
to satisfy additional local conservation laws (besides the conservation of the
number of particles)? Systems with additional conservation laws have attracted
much interest in recent years. Many studies concerned quantum systems with
infinitely many conservation laws, specifically integrable systems in one
spatial dimension. The corresponding generalized hydrodynamics is an active
research area, see Refs. Doyon20c ; Dubail22  for reviews. The generalized
hydrodynamics was also used for the description of the hydrodynamic behavior
in classical integrable systems in one spatial dimension Doyon19 ; Doyon20a ;
Spohn23 . Quantum higher-moment conserving models attracted much attention in
the context of fractons Chamon05 ; Pretko17 ; Gromov19 ; Williamson19 ; Lake22
; Williamson23 ; Gromov24 ; see Refs. Lucas20 ; Knap20 ; Lucas21 ; Surowka23 ;
Morningstar23  for a sample of articles about the emerging ‘fracton
hydrodynamics’. Lattice gases with additional conservation laws provide a
convenient and flexible theoretical platform for a systematic derivation
(rather than a postulation) of a ‘fracton hydrodynamics’.

Han et al. HLR2024  have recently introduced an interesting stochastic lattice
gas model which manifestly conserves, in addition to the number of particles,
their dipole moment. In one spatial dimension, the large-scale deterministic
limit of this model is governed by a nonlinear partial differential equation
(PDE) of fourth order:

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle-D\partial_{x}^{2}\left[\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}\right]$
(4) $\displaystyle=$
$\displaystyle-D\partial_{x}^{2}\left(\rho^{2}\partial_{x}^{2}\ln\rho\right)$
$\displaystyle=$
$\displaystyle-D\left[\rho\partial_{x}^{4}\rho-(\partial_{x}^{2}\rho)^{2}\right]\,.$

As one can see, here the simple first-order continuity equation (1) gives way
to the second-order equation Pretko17

$\partial_{t}\rho+\partial_{x}^{2}J=0\,,$ (5)

with the current

$J=D\rho^{2}\partial_{x}^{2}\ln\rho\,.$ (6)

Equation (4) can be generalized to an arbitrary spatial dimension HLR2024 :

$\partial_{t}\rho=-D[\rho\Delta^{2}\rho-(\Delta\rho)^{2}]\,,$ (7)

where $\Delta\equiv\nabla^{2}$ is the Laplace’s operator. The transport
coefficient $D$ comes from the microscopic model, and it has the units of
$\text{length}^{d+4}/\text{time}$, where $d$ is the dimension of space jamming
.

Because of an interplay of the nonlinearity and the fourth-order spatial
derivative, Eqs. (4) and (7) exhibit unusual macroscopic scaling behaviors
HLR2024 . Here we will further explore these scaling behaviors by deriving
some similarity solutions to these equations, which describe relaxation of the
dipole-conserving gas (DCG) in several standard settings. As we will see, the
scaling properties of the relaxation dynamics crucially depend on whether the
total number of particles is infinite or finite.

The main focus of this work, however, is on large-scale fluctuations in the
DCG. A study of fluctuations obviously requires going beyond the deterministic
limit, described by Eqs. (4) and (7). Han et al. HLR2024  have already made
this important step by deriving, from a microscopic lattice gas model, a
Langevin equation for this system, see Eq. (19) below. In addition to the
terms present in the deterministic equation (4), this stochastic PDE also
includes a noise term. A similar-in-spirit Langevin description of the mass-
only conserving lattice gases is known by the name of ‘fluctuational
hydrodynamics’ Spohn91 ; Anna ; Kipnis99 ; Schutz00 . Starting from the
Langevin equation (19), here we develop a macroscopic fluctuation theory
(MFT), which is suitable for studying large deviations of different
fluctuating quantities in the DCG. In the mass-only conserving lattice gases
the corresponding MFT was developed by Jona-Lasinio et al., see Ref. bertini
for a review, and it has been employed and further developed in numerous
subsequent works.

Here we use the MFT to establish the form of the Boltzmann-Gibbs distribution
for the DCG at equilibrium. We also apply the MFT to the calculation of the
variance of nonequilibrium fluctuations of the excess number of particles on
the positive semi-axis when starting from a (either deterministic, or random)
constant density at $t=0$. Finally, we determine the probability distribution
of, and the most probable density history leading to, a _large deviation_ in
the form of _void_ of a given size in an initially uniform DCG at equilibrium.

Here is a plan of the remainder of this paper. In Sec. II we present some
similarity solutions of the deterministic Eqs. (4) and (7), which involve
infinite and finite mass, and discuss their properties. Sections III and IV
deal with fluctuations in the DCG. Starting from fluctuational hydrodynamics,
as described by the Langevin equation (19), we introduce in Sec. III the
problem of full statistics of the excess number of particles on the positive
semi-axis. Using this setting as an example, we formulate the MFT of large
deviations in the DCG and calculate the variance of the excess number of
particles. Section IV is devoted to the MFT at equilibrium. Here we introduce
the free energy density of the DCG at equilibrium and determine the
probability distribution of, and the most probable density history leading to,
the formation of a void in an initially uniform gas. Section V presents a
brief summary and discussion of our main results. Some technical details of
the derivation of the MFT equations and boundary conditions are relegated to
Appendix A. In Appendix B we present an independent calculation of the
variance of the particle excess directly from the Langevin equation (19).

##  II Deterministic relaxation

###  II.1 Infinite-mass scaling

To start with, let us study expansion of the DCG into vacuum. Suppose that the
initial gas density has the form of a step-function:

$\rho(x,t=0)=\rho_{0}\theta(-x)\,.$ (8)

The relaxation of this system is described by the following similarity soluton
of Eq. (4):

$\rho(x,t)=\rho_{0}\,R\left[\frac{x}{(\rho_{0}Dt)^{1/4}}\right]\,.$ (9)

In this case the dynamical exponent $4$ is the same as in the linear fourth-
order equation $\partial_{t}u=-D_{0}\partial_{x}^{4}u$ originally studied in
the context of surface diffusion Mullins .

The dimensionless scaling function $R(\xi)$ obeys an ordinary differential
equation (ODE),

$\frac{1}{4}\xi
R^{\prime}(\xi)=\frac{d^{2}}{d\xi^{2}}\left[RR^{\prime\prime}-(R^{\prime})^{2}\right]\,.$
(10)

The boundary conditions are $R(-\infty)=1$, $R^{\prime}(-\infty)=0$ and
$R(+\infty)=R^{\prime}(+\infty)=0$. The scaling function $R(\xi)$ can be
obtained by solving Eqs. (10) with these boundary conditions numerically.
Alternatively, we can solve numerically the full time-dependent PDE (4) after
bringing it to a dimensionless form by rescaling $\rho_{0}x\to x$,
$\rho_{0}^{5}Dt\to t$ and $\rho/\rho_{0}\to\rho$. Figure 1 gives an example of
such a time-dependent solution for the rescaled initial condition
$\rho(x,t=0)=1-\tanh(15\,x)$. The top panel shows this initial condition and
the resulting density profiles at rescaled times $t=5$, $10$ and $15$. The
bottom panel shows the same three density profiles, but plotted against the
similarity coordinate $\xi$. As one can see, the profiles collapse into a
single curve, which describes the scaling function $R(\xi)$. Salient features
of this similarity solution are its oscillatory decay at $x\to-\infty$ and its
semi-compact support: the solution is defined for $-\infty<\xi<\xi_{*}\simeq
2.5$. The asymptotic of the solution near the edge is
$R(\xi)\simeq(\xi_{*}/48)(\xi_{*}-\xi)^{3}$, so that the first and second
derivatives of $R$ vanish at $\xi=\xi_{*}$ alongside with $R$.

Figure 1: Top panel: Numerical solution of Eq. (4) at $t=1,5$ and $15$ and
$\rho(x,0)=\rho_{0}[1-\tanh(15\,x)]$. Bottom panel: the scaling function
$R(\xi)$ obtained by plotting the same three profiles against the similarity
coordinate $\xi$.

###  II.2 Finite-mass scaling

The long-time evolution of a system with a finite number $N$ of particles is
described by the similarity solution

$\rho(x,t)=\frac{N^{4/5}}{(Dt)^{1/5}}\,R\left[\frac{x}{(NDt)^{1/5}}\right]\,,$
(11)

which exhibits a different dynamical exponent $5$ HLR2024 . The dimensionless
scaling function $R(\xi)$ obeys the normalization condition

$\int_{-\infty}^{\infty}R(\xi)\,d\xi=1\,,$ (12)

following from the conservation of the total number of particles. The form
(11) is exact if the system starts from the initial condition

$\rho(x,t=0)=N\delta(x)\,.$ (13)

Otherwise, Eq. ((11)) describes a long-time asymptotic of the solution
Barenblatt . For the scaling function $R(\xi)$ we obtain an ODE

$\frac{1}{5}\frac{d}{d\xi}\left(\xi\,R\right)=\frac{d^{2}}{d\xi^{2}}\left[RR^{\prime\prime}-(R^{\prime})^{2}\right]\,.$
(14)

Integrating once, we obtain

$\frac{1}{5}\xi\,R=\frac{d}{d\xi}\left[RR^{\prime\prime}-(R^{\prime})^{2}\right]=RR^{\prime\prime\prime}-R^{\prime}R^{\prime\prime}\,,$
(15)

while the integration constant must be zero. Since $R(\xi)$ is an even
function, we can solve Eq. (15) on the half-line $\xi>0$ with the boundary
conditions

$R(0)=a>0\,,\quad R^{\prime}(0)=0\,,\quad R(\xi\to\infty)=0\,.$ (16)

The solution must be nonnegative, and the a priori unknown constant $a$ is to
be determined from the normalization condition (12). The nonnegativity of the
solution and the boundary condition $R(\xi\to\infty)=0$ demand that the
solution have a compact support. At the edges of support both $R(\xi)$ and
$R^{\prime}(\xi)$ must vanish, thus providing continuity of the flux, see Eq.
(15).

Figure 2: The scaling function $R(\xi)$ of the similarity solution (11) found
by numerically solving the problem (15) and (16) (solid line) and the full
time-dependent PDE (4) with a localized initial condition (dashed line). Only
the positive-$\xi$ region is shown.

As one can check, the ODE (15) remains invariant under rescaling $\xi\to
C^{-1/4}\xi$ and $R\to C^{-1}R$, where $C>0$. Therefore, once we have found
the solution $R_{1}(\xi)$ of the problem (15) and (16) for $a=1$, we can find
the solution $R_{a}(\xi)$ for arbitrary $a$ by the rescaling transformation

$R_{a}(\xi)=aR_{1}\left(\frac{\xi}{a^{1/4}}\right)\,.$ (17)

Using Eqs. (12) and (17), we obtain

$a=\left[2\int_{0}^{\infty}R_{1}(\xi)\,d\xi\right]^{-4/5}\,,$ (18)

so what remains is to find $R_{1}(\xi)$. This can be achieved numerically by
the shooting method. We set $a=1$ and trade the boundary condition at infinity
$R_{1}(\infty)=0$ for the condition $R^{\prime\prime}(0)=\gamma$, where
$\gamma<0$ serves as the shooting parameter. Having found $R_{1}(\xi)$ and
employing Eqs. (17) and (18), we obtained the numerical solution shown in Fig.
2. Here $a=0.356\dots$, while the edges of support are at
$|\xi|=\xi_{*}=3.140\dots$. (The proximity of the latter number to $\pi$
raises curiosity but most likely coincidental.) The asymptotic of $R(\xi)$
near the edges is $R(\xi)\simeq(\xi_{*}/60)(\xi_{*}-|\xi|)^{3}$, so that the
first and second derivatives of $R$ vanish alongside with $R$ at
$|\xi|=\xi_{*}$ similarly to the step-like solution of the previous
subsection.

##  III Macroscopic fluctuation theory. Fluctuations of excess number of
particles

Fluctuation hydrodynamics of the DCG is described by the Langevin equation
which has been recently derived in Ref. HLR2024 :

$\partial_{t}\rho=-D\partial_{x}^{2}\left[\rho\partial_{x}^{2}\rho-(\partial_{x}\rho)^{2}\right]+\sqrt{2D}\,\partial_{x}^{2}\left[\rho\eta(x,t)\right]\,,$
(19)

where $\eta(x,t)$ is a white Gaussian noise,
$\langle\eta(x_{1},t_{1})\eta(x_{2},t_{2})\rangle=\delta(x_{1}-x_{2})\delta(t_{1}-t_{2})$,
and we confine ourselves to one spatial dimension.

As we will show in Sec. IV, the equilibrium state of the DCG can be described
by the Boltzmann-Gibbs distribution with a well-defined free-energy density
$F(\rho)$. That is, when the gas is in equilibrium at density $\rho_{0}$, the
probability density of observing an arbitrary density profile $\rho_{*}(x)$ is
given by $-\ln{\mathcal{P}}[\rho_{*}(x)]\simeq S_{\text{eq}}$, where

$S_{\text{eq}}\\!\simeq\\!\int_{-\infty}^{\infty}dx\left[F(\rho_{*}(x))-F(\rho_{0})-F^{\prime}(\rho_{0})(\rho_{*}(x)-\rho_{0})\right].$
(20)

Remarkably, the free energy density of this gas,

$F(\rho)=\rho\ln\rho-\rho\,,$ (21)

coincides with that of the lattice gas of noninteracting random walkers (RWs).
The latter is described by the more familiar second-order Langevin equation
Spohn91

$\partial_{t}\rho=\partial_{x}\left[D\partial_{x}\rho+\sqrt{2D\rho}\,\eta(x,t)\right]\,.$
(22)

Plugging Eq. (21) for $F(\rho)$ into Eq. (20), we obtain

$S_{\text{eq}}=\int_{-\infty}^{\infty}dx\left[\rho_{*}(x)\ln\frac{\rho_{*}(x)}{\rho_{0}}+\rho_{0}-\rho_{*}(x)\right]\,,$
(23)

which coincides with the corresponding expression for the gas of
noninteracting RWs.

Large deviations in the DCG can be described by the MFT bertini . The MFT is a
weak-noise theory which is based on the Langevin equation (19) and relies on a
problem-specific small parameter, that we identify below. When such a
parameter is present, the probability distribution of interest can be
approximately determined via a saddle-point evaluation of the exact path
integral corresponding to Eq. (19) and the problem-specific constraints.

We will introduce the MFT for the DCG on the example of the statistics of the
excess number of particles $K$ on the positive semi-axis at specified time
$T$. One can consider two different settings. In the first of them we start
from a constant gas density $\rho_{0}$ and condition the process on the value
of the integral

$\int_{0}^{\infty}\left[\rho(x,T)-\rho_{0}\right]\,dx=K.$ (24)

at time $T$. This setting corresponds to infinite total mass of the system.

In the second setting (a finite total mass) there is a large but finite number
of particles, $N\gg 1$, in the system. Here the particle excess condition at
$t=T$ is

$\int_{0}^{\infty}\rho(x,T)\,dx-\frac{N}{2}=K\,,$ (25)

where $-N/2\leq K\leq N/2$.

As we will see shortly, the scaling behavior of the particle excess statistics
in the two settings is quite different. In the infinite-mass setting the
result also depends on whether the initial condition is quenched: that is,
deterministically prepared, or annealed: that is, randomly sampled from the
equilibrium distribution of the gas at density $\rho_{0}$.

###  III.1 Infinite-mass scaling

Let us rescale the variables: $t/T\to t$, $x/(\rho_{0}DT)^{1/4}\to x$, and
$\rho/\rho_{0}\to\rho$. In the new variables Eqs. (19) and (24) become

$\partial_{t}\rho=-\partial_{x}^{2}\left(\rho^{2}\partial_{x}^{2}\ln\rho\right)+\sqrt{2}\,\left(\rho_{0}^{5}DT\right)^{-1/8}\partial_{x}^{2}\left[\rho\eta(x,t)\right]$
(26)

and

$\int_{0}^{\infty}\left[\rho(x,1)-1\right]\,dx=j\equiv\frac{K}{(\rho_{0}^{5}DT)^{1/4}},$
(27)

respectively. At long time the noise becomes effectively weak. The presence of
a large parameter $(\rho_{0}^{5}DT)^{1/4}\gg 1$, which is a typical number of
particles within a region with the length scale $(\rho_{0}DT)^{1/4}$, makes it
possible to develop the MFT, that is to perform a saddle-point evaluation of
the exact path integral for Eq. (26) with account of the constraint (27). The
calculation boils down to minimization of the action functional (see Appendix
A) and leads to the following Hamilton’s equations that describe the optimal
(that is, the most likely) path of the system conditioned on Eq. (27):

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle\partial_{x}^{2}\left[-\rho\partial_{x}^{2}\rho+(\partial_{x}\rho)^{2}+2\rho^{2}v\right]\,,$
(28) $\displaystyle\partial_{t}v$ $\displaystyle=$
$\displaystyle\partial_{x}^{2}\left[v\partial_{x}^{2}\rho+\partial_{x}^{2}(v\rho)+2\partial_{x}(v\partial_{x}\rho)-2v^{2}\rho\right]\,.$
(29)

The boundary condition at $t=T$ is

$v(x,1)=\lambda\,\delta^{\prime}(x)\,,$ (30)

where $\lambda$ is a Lagrange multiplier, introduced when accommodating the
constraint (27), and $\delta^{\prime}(x)$ is the $x$-derivative of the delta
function. The quenched initial condition is

$\rho(x,0)=1\,.$ (31)

The annealed initial condition is introduced shortly.

Once the optimal path is found, the probability distribution is given, up to a
pre-exponent, by the action along the optimal path:
$-\ln\mathcal{P}(K,T,\rho_{0})\simeq S$, where

$S=\left(\rho_{0}^{5}DT\right)^{1/4}s(j)\,,$ (32)

where

$s(j)=\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\,\rho^{2}(x,t)v^{2}(x,t)\,.$
(33)

The rescaled excess of the number of particles $j$ is defined in Eq. (27).

For the annealed initial condition with _average_ rescaled density $1$ the
full action includes the cost of creating the optimal initial condition
$\rho(x,0)$: $S_{\text{annealed}}=S+S_{0}$, where $S$ is the dynamical action,
described by Eq. (32), and

$S_{0}=\int_{-\infty}^{\infty}dx\left\\{F[\rho(x,0)]-F(1)-F^{\prime}(1)[\rho(x,0)-1]\right\\},$
(34)

As a result, Eq. (31) gives way to a different condition DG2009b , which
describes a relation between the a priori unknown $\rho(x,0)$ and $v(x,0)$:

$v(x,0)-\frac{d^{2}}{dx^{2}}\ln\rho(x,0)=\lambda\delta^{\prime}(x)\,.$ (35)

###  III.2 Finite-mass scaling

In this case the parameter $\rho_{0}$ is absent, and the rescaling of
variables is different:

$\frac{t}{T}\to t,\quad\frac{x}{(NDT)^{1/5}}\to
x\quad\text{and}\quad\frac{(DT)^{1/5}\rho}{N^{4/5}}\to\rho\,.$ (36)

Notice that the infinite-mass dynamical exponent $4$ gives way to the finite-
mass exponent $5$ HLR2024 . The rescaled Eqs. (19) and (25) become

$\partial_{t}\rho=-\partial_{x}^{2}\left(\rho^{2}\partial_{x}^{2}\ln\rho\right)+\sqrt{\frac{2}{N}}\,\partial_{x}^{2}\left[\rho\eta(x,t)\right]\,,$
(37)

and

$\int_{0}^{\infty}\rho(x,1)\,dx-\frac{1}{2}=\frac{K}{N}\equiv j\,,$ (38)

where $|j|\leq 1/2$. Here the saddle-point expansion, leading to the MFT,
relies on the large parameter $N\gg 1$.

The rescaled MFT equations coincide with Eqs. (28) and (29). The boundary
condition (30) is also the same, but Eq. (31) gives way to the condition

$\rho(x,0)=\delta(x)\,.$ (39)

Here the probability distribution of the number of particles, transferred to
the right, is independent of the observation time $T$:
$-\ln\mathcal{P}(K,N)\simeq S$, where

$S=Ns(j)\,,$ (40)

and $s(j)$ is again described by Eq. (33).

###  III.3 Variance of particle excess

The statistics of typical (that is small) fluctuations of the excess number of
particles $K$ is Gaussian, and its variance scales as a characteristic number
of particles involved. In the infinite-mass case this is the typical number of
particles over the dynamical length scale $(\rho_{0}DT)^{1/4}$, that is
$\text{var}_{K}\sim(\rho_{0}^{5}DT)^{1/4}$. In the finite-mass case it is
simply $\text{var}_{K}\sim N$. Essentially, the role of theory (a first-order
perturbation theory in $|\lambda|\ll 1$, developed in Ref. KrMe ) is to
provide the numerical coefficients $O(1)$ in these expressions.

For the infinite-mass case the calculations are straightforward. Upon
linearization with respect to $\lambda$, the MFT equations (28) and (29)
become

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle-\rho_{0}\partial_{x}^{4}\delta\rho+\rho_{0}^{2}\partial_{x}^{2}v\,,$
(41) $\displaystyle\partial_{t}v$ $\displaystyle=$
$\displaystyle\partial_{x}^{4}v\,.$ (42)

We can solve Eq. (42) backward in time with the “initial condition” (30). The
solution can be obtained by differentiating with respect to $x$ the previously
known solution for the initial condition in the form of a delta-function (see
e.g. Ref. MV2016  and Eq. (74) in Appendix B below). The result has the
similarity form

$v(x,t)=\frac{\lambda}{(1-t)^{1/2}}\,V\left[\frac{x}{(1-t)^{1/4}}\right]\,.$
(43)

The scaling function $V(z)$ can be expressed via the hypergeometric function
${}_{0}F_{2}$:

$V(z)=-\frac{z\Gamma\left(\frac{3}{4}\right)\,_{0}F_{2}\left(\frac{5}{4},\frac{3}{2};\frac{z^{4}}{256}\right)}{4\pi}+\frac{z^{3}\Gamma\left(\frac{5}{4}\right)\,_{0}F_{2}\left(\frac{3}{2},\frac{7}{4};\frac{z^{4}}{256}\right)}{24\pi}-\frac{z^{5}\Gamma\left(\frac{3}{4}\right)\,_{0}F_{2}\left(\frac{9}{4},\frac{5}{2};\frac{z^{4}}{256}\right)}{960\pi}\,.$
(44)

Now we can evaluate the rescaled dynamical action (33) in terms of $\lambda$.
Within the linear theory in $\lambda$, we should replace the density
$\rho(x,t)$ in Eq. (33) by $1$. We obtain

$\displaystyle s(\lambda)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\,v^{2}(x,t)$ (45)
$\displaystyle=$
$\displaystyle\lambda^{2}\int_{0}^{1}\frac{dt}{(1-t)^{3/4}}\int_{-\infty}^{\infty}dz\,V^{2}(z)\,.$

The integral over $t$ gives $4$. The integral over $z$ can be evaluated
numerically: $\int_{-\infty}^{\infty}V^{2}(z)\,dz=\alpha=0.05798\dots$.
Overall, we obtain $s=4\alpha\lambda^{2}$.

####  III.3.1 Quenched initial condition

For the quenched setting, the formula $s=4\alpha\lambda^{2}$ suffices for
expressing the action in terms of $j$. Indeed, using the “shortcut relation”
$ds/dj=\lambda$ (see, e.g. Ref. Vivo ), we obtain

$\frac{ds}{d\lambda}\frac{d\lambda}{dj}=2\alpha\lambda\,\frac{d\lambda}{dj}=\lambda\,,$
(46)

which gives $\lambda=j/(8\alpha)$. As a result,

$s=4\alpha\lambda^{2}=\frac{j^{2}}{16\alpha}\,.$ (47)

Back to the dimensional variables, we obtain the variance of the typical
fluctuations of $K$:

$\text{var}K=8\alpha(\rho_{0}^{5}DT)^{1/4}\,.$ (48)

Exactly the same result (48) for the quenched initial condition can be
obtained directly from the linearized version of the Langevin equation (26).
This calculation is presented in Appendix B. The MFT, however, also enables
one to calculate the variance for the annealed initial condition, where a
direct calculation with the linearized Langevin equation does not seem to be
available.

####  III.3.2 Annealed initial condition

In the annealed setting one should also take into account (the small-$\lambda$
expansions of) Eqs. (34) and (35), which are the following:

$s_{0}=\frac{1}{2}\int_{-\infty}^{\infty}\delta\rho^{2}(x,t=0)\,dx$ (49)

and

$v(x,t=0)-\frac{d^{2}}{dx^{2}}\delta\rho(x,t=0)=\lambda\delta^{\prime}(x)\,,$
(50)

respectively, where $\delta\rho(x,t=0)$ is the a priori unknown small
perturbation on the background of the constant density $\rho=\rho_{0}$, and
$v(x,t=0)$ is also small.

Evaluating Eq. (43) at $t=0$, plugging the result into Eq. (50) and
integrating the latter equation twice over $x$, we determine the optimal
initial density field $\rho(x,0)\simeq 1+\delta\rho(x,0)$, where

$\delta\rho(x,0)=\lambda\left[\int_{-\infty}^{x}dy\int_{-\infty}^{y}dz\,F(z)-\theta(x)\right]\,,$
(51)

and $\theta(x)$ is the step-function. The double integral in this expression
can be evaluated with “Mathematica” analytically, but the result is too
cumbersome to present it here. Plugging it into Eq. (49) and evaluating the
resulting integral numerically, we obtain

$s_{0}=\beta\,\lambda^{2}\,,$ (52)

where $\beta=0.1581\dots$. The total annealed action is, therefore,

$s(\lambda)=(4\alpha+\beta)\lambda^{2}=\frac{j^{2}}{4(4\alpha+\beta)}\,,$ (53)

where $\alpha=0.05798\dots$, as obtained above, and we have again used the
shortcut relation $ds/dj=\lambda$ to express the action through $j$.

At given $j$, the annealed action (53) is smaller than the quenched action
(47). Back in the original variables, the variance of the typical fluctuations
of $K$ in the annealed case,

$\text{var}K=(8\alpha+2\beta)(\rho_{0}^{5}DT)^{1/4}\,,$ (54)

is larger than that in the quenched case [see Eq. (48)], as to be expected on
physical grounds.

##  IV MFT of large deviations at equilibrium

Now we suppose that the DCG is at equilibrium and, using the MFT, evaluate the
probability density of observing a specified density profile $\rho_{*}(x)$.
Importantly, the Hamiltonian MFT equations are still Eqs. (28) and (29) (where
we return to the dimensional variables except setting $D=1$). The boundary
conditions in time, however, become $\rho(x,t=-\infty)=\rho_{0}$ and
$\rho(x,t=0)=\rho_{*}(x)$ KrMe ; KMSvoid .

At the microscopic level, the DCG obeys detailed balance. At the macroscopic
level the detailed balance manifests itself as the Onsager principle of
reversibility. In particular, the optimal activation path $\rho(x,t)$, leading
from $\rho(x,t=-\infty)=\rho_{0}$ to $\rho(x,t=0)=\rho_{*}(x)$, must coincide
with the _time-reversed_ relaxation path from $\rho(x,t=-\infty)=\rho_{*}(x)$
to $\rho(x,t=\infty)=\rho_{0}$. That is, the optimal path must obey the
equation

$\partial_{t}\rho=\partial_{x}^{2}\left(\rho^{2}\partial_{x}^{2}\ln\rho\right)\,.$
(55)

Combining this equation with Eq. (28), we obtain the important relation

$v(x,t)=\partial_{x}^{2}\ln\rho(x,t)\,,$ (56)

which describes the equilibrium manifold of this system. The relation (56)
brings about two important consequences. Firstly, as one can check by a direct
substitution, Eq. (29) is now obeyed automatically. Secondly, the mechanical
action can be calculated as follows:

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dt\int_{-\infty}^{\infty}dx\,\rho^{2}(x,t)v^{2}(x,t)$
(57) $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dt\int_{-\infty}^{\infty}dx\,\rho^{2}(x,t)\left(\partial_{x}^{2}\ln\rho\right)^{2}\,.$

After two integrations by part over $x$, this expression becomes

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dt\int_{-\infty}^{\infty}dx\ln\rho\,\partial_{x}^{2}\left(\rho^{2}\partial_{x}^{2}\ln\rho\right)\,.$
(58)

By virtue of Eq. (55), this expression can be recast as

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx\int_{-\infty}^{0}dt\,\ln\rho\,\partial_{t}\rho$
(59) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx\int_{-\infty}^{0}dt\,\partial_{t}\left(\rho\ln\rho-\rho\right)\,.$

Performing the integration over time, we arrive at the announced Boltzmann-
Gibbs relation (23).

As a simple but instructive example, let us evaluate the probability
distribution $\mathcal{P}_{\text{void}}$ of observing a void of size $2L$ in a
uniform gas at equilibrium:

$\rho_{*}(x)=\begin{cases}0,&|x|<L,\\\ \rho_{0},&|x|>L.\end{cases}$ (60)

Using Eq. (23), we obtain $\mathcal{P}_{\text{void}}\sim\exp(-2\rho_{0}L)$, as
in the gas of noninteracting RWs KMSvoid . The optimal path of the system
toward the void formation, $\rho(x,t)$ is obtained by the time reversal of the
relaxation dynamics of the void, see Eq. (55). We computed this optimal path
numerically, and the results are shown in Fig. 3. Noticeable are spatial
oscillations of the density, which are absent in the optimal path of the void
formation in the gas of noninteracting RWs.

Figure 3: The optimal path of void formation at $t=0$ and $\rho_{0}=L=D=1$.
Shown is the optimal density profile at times $t=-\infty$ (black), $-6$
(blue), $-2$ (brown), $-1/5$ (magenta) and $0$ (red).

##  V Summary and Discussion

We have considered some basic macroscopic relaxation and fluctuation
properties of the mass- and dipole-conserving stochastic lattice gas. We have
extended the MFT approach to this system. Using some carefully selected
examples, we have demonstrated how one can apply the MFT both to typical
fluctuations of the DPG, and to its large deviations. We hope that the MFT
formalism will find applications in particular settings relevant to
experiments.

We would like to conclude this work with a fascinating example. Let us briefly
recall one basic property of diffusive lattice gases that conserve only the
number of particles. The fluctuational hydrodynamics of such gases is
described by the Langevin equation Spohn91

$\partial_{t}\rho=\partial_{x}\left[\mathcal{D}(\rho)\partial_{x}\rho+\sqrt{\sigma(\rho)}\,\eta(x,t)\right]\,,$
(61)

which generalizes Eq. (22) to a broad class of diffusive lattice gases. The
free energy density $F(\rho)$ of this class of gases is determined by the
diffusion coefficient $\mathcal{D}(\rho)$ and the mobility $\sigma(\rho)$ via
the Einstein relation Spohn91

$F^{\prime\prime}(\rho)=\frac{2\mathcal{D}(\rho)}{\sigma(\rho)}\,.$ (62)

In the particular case of noninteracting RWs, see Eq. (22), one has
$\mathcal{D}(\rho)=D=\text{const}$, and $\sigma(\rho)=2D\rho$. Then Eq. (62)
yields the free energy density $F(\rho)=\rho\ln\rho-\rho$. The remarkable
coincidence of the free energy density (21) for the DCG with this expression
may suggest the validity of the Einstein relation for the DCG at the level of
the Langevin equation (26). Indeed, if, by analogy with Eq. (61), we identify
the functions $\mathcal{D}(\rho)=D\rho$ (which enters the highest-derivative
term) and $\sigma(\rho)=2D\rho^{2}$, then Eq. (62) correctly reproduces the
free energy density (21). This observation hints at possible additional
surprises in the studies of the dipole-conserving lattice gases.

## Acknowledgments

The author is very grateful to P. L. Krapivsky for the stimulating discussions
and comments. This research was supported by the Israel Science Foundation
(Grant No. 1499/20).

##  Appendix A Derivation of Eqs. (28) and (29) and boundary conditions

Our starting point is Eq. (26). One way to deal with the second derivative in
the noise term is to define the “second-order potential” $\psi(x,t)$ by the
relation $\rho(x,t)=\partial_{x}^{2}\psi(x,t)$ and rewrite Eq. (26) in terms
of $\psi(x,t)$:

$\partial_{t}\psi=-\partial_{x}^{2}\psi\,\partial_{x}^{4}\psi+(\partial_{x}^{3}\psi)^{2}+\sqrt{2\epsilon}\,\partial_{x}^{2}\psi\,\eta(x,t)\,,$
(63)

where $\epsilon=\left(\rho_{0}^{5}DT\right)^{-1/4}\ll 1$. The probability
density of a realization of the (dimensionless) Gaussian white noise
$\eta(x,t)$ is

$\mathcal{P}\left[\eta\right]\sim\exp\left(-\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\frac{\eta^{2}}{2}\right).$
(64)

Expressing $\xi$ through $\psi$ and its derivatives from Eq. (63), we obtain

$\mathcal{P}[\psi(x,t)]\sim\exp\left\\{-\frac{1}{4\epsilon}\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\frac{\left[\partial_{t}\psi+\partial_{x}^{2}\psi\,\partial_{x}^{4}\psi-(\partial_{x}^{3}\psi)^{2}\right]^{2}}{(\partial_{x}^{2}\psi)^{2}}\right\\}\,.$
(65)

Employing the small parameter $\epsilon\ll 1$, we can apply the saddle-point
approximation. The calculations boil down to a minimization of the action
functional

$s_{\lambda}[\psi(x,t)]=\frac{1}{4}\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\frac{\left[\partial_{t}\psi+\partial_{x}^{2}\psi\,\partial_{x}^{4}\psi-(\partial_{x}^{3}\psi)^{2}\right]^{2}}{(\partial_{x}^{2}\psi)^{2}}-\lambda\int_{-\infty}^{\infty}\theta(x)\left[\partial_{x}^{2}\psi(x,1)-\partial_{x}^{2}\psi(x,0)\right]\,dx\,,$
(66)

where the constraint (27) is incorporated with the help of a Lagrange
multiplier $\lambda$. Let us introduce the second derivative
$\partial_{x}^{2}p$ of the canonical momentum density $p(x,t)$ by
differentiating the Lagrangian of the action functional (66) with respect to
$\psi_{t}$ LL . We obtain

$\partial_{x}^{2}p=\frac{\partial_{t}\psi+\partial_{x}^{2}\psi\,\partial_{x}^{4}\psi-(\partial_{x}^{3}\psi)^{2}}{2(\partial_{x}^{2}\psi)^{2}}\,.$
(67)

Denoting $\partial_{x}^{2}p(x,t)\equiv v(x,t)$ and going back from $\psi(x,t)$
to the density $\rho(x,t)$, we obtain Eq. (28).

The Hamiltonian $H$, corresponding to the first term of the action functional
(66), is the following LL :

$H=\int_{-\infty}^{\infty}dx\left[v^{2}\partial_{x}^{2}\psi-v\partial_{x}^{2}\psi\,\partial_{x}^{4}\psi+v(\partial_{x}^{3}\psi)^{2}\right]\,.$
(68)

Equation (29) can be obtained by taking the minus variational derivative of
this Hamiltonian with respect to $\psi$, and then going back from $\psi$ to
$\rho$.

This derivation, however, ignores possible boundary terms at $t=0$ and $t=1$,
arising when calculating the linear variation of the constrained action
$s_{\lambda}[\psi(x,t)]$. One boundary term at $t=1$ comes from the term
proportional to $\lambda$ in Eq. (66). An additional term at $t=1$ results
from integration by parts of the term that includes
$\partial_{t}\delta\psi(x,t)$, where $\delta\psi(x,t)$ is the linear variation
of $\psi(x,t)$. These two terms yield the boundary condition (30).

For the quenched initial condition (31) there are no boundary terms at $t=0$.
In the annealed setting there are two boundary terms at $t=0$. One of them
arises from integration by parts of the term that includes
$\partial_{t}\delta\psi(x,t)$, while the other comes from the variation of the
term proportional to $\lambda$ in Eq. (66). (Indeed, the initial density
$\rho(x,t=0)$ in this case is a priori unknown.) These two term lead to the
annealed initial condition (35).

##  Appendix B Variance of particle excess number from the Langevin equation

Typical, small fluctuations of $K$ are caused by typical, small fluctuations
of the density profile $\rho(x,t)$. In order to account for such fluctuation,
one can linearize the Langevin equation (26) around the equilibrium state
$\rho=1$: $\rho(x,t)=1+u(x,t)$, where $|u|\ll 1$ KrMe . The linearized
equation is

$\partial_{t}u(x,t)=-\partial_{x}^{4}u(x,t)+\sqrt{2\epsilon}\,\partial_{x}^{2}\eta(x,t)\,,$
(69)

where $\epsilon=\left(\rho_{0}^{5}DT\right)^{-1/4}\ll 1$. The rescaled
particle excess number at $t=1$, see Eq. (27), becomes

$j\equiv\epsilon K=\int_{0}^{\infty}u(x,1)\,dx\,.$ (70)

Using the new variable,

$\psi(x,t)=\int_{-\infty}^{x}dy\int_{-\infty}^{y}dz\,u(z,t)\,,$ (71)

we can rewrite Eq. (69) as

$\partial_{t}\psi(x,t)=-\partial_{x}^{4}\psi(x,t)+\sqrt{2\epsilon}\,\eta(x,t)\,,$
(72)

and formally solve it for a given realization of noise $\eta(x,t)$ and the
quenched initial condition $\psi(x,t=0)=0$. In particular, at rescaled time
$t=1$ we obtain

$\psi(x,t=1)=\sqrt{2\epsilon}\int_{0}^{1}dt\int_{-\infty}^{\infty}dy\,\eta(y,t)G(x-y,1-t)\,,$
(73)

where $G(z,\tau)$ is the Green’s function of the homogeneous equation,
corresponding to Eq. (69). That is, $G(z,\tau)$ solves the equation
$\partial_{\tau}G(z,\tau)=-\partial_{z}^{4}G(z,\tau)$ with the initial
condition $G(z,\tau=0)=\delta(z)$. The explicit form of $G(z,\tau)$ is the
following:

$G(z,\tau)=\frac{\Gamma\left(\frac{5}{4}\right)\,_{0}F_{2}\left(\frac{1}{2},\frac{3}{4};\frac{z^{4}}{256\tau}\right)}{\pi\tau^{1/4}}-\frac{z^{2}\Gamma\left(\frac{3}{4}\right)\,_{0}F_{2}\left(\frac{5}{4},\frac{3}{2};\frac{z^{4}}{256\tau}\right)}{8\pi\tau^{3/4}}\,.$
(74)

Plugging Eq. (73) into the constraint (70) and performing the integration over
$x$, we obtain the fluctuating particle excess number at $t=1$:
$j=-\partial_{x}\psi(x=0,t=1)$. The variance of these fluctuations is

$\text{Var}_{K}=\langle\left[\partial_{x}\psi(0,1)\right]^{2}\rangle=2\epsilon\int_{0}^{1}dt_{1}\int_{0}^{1}dt_{2}\int_{-\infty}^{\infty}dy_{1}\int_{-\infty}^{\infty}dy_{2}\left\langle\eta(y_{1},t_{1})\eta(y_{2},t_{2})\right\rangle\frac{\partial
G(y_{1},1-t_{1})}{\partial y_{1}}\times\frac{\partial
G(y_{2},1-t_{2})}{\partial y_{2}}\,.$ (75)

Since
$\left\langle\eta(y_{1},t_{1})\eta(y_{2},t_{2})\right\rangle=\delta(y_{1}-y_{2})\delta(t_{1}-t_{2})$,
Eq. (75) simplifies to

$\text{Var}_{K}=2\epsilon\int_{0}^{1}dt\int_{-\infty}^{\infty}dx\left[\frac{\partial
G(x,1-t)}{\partial x}\right]^{2}\,.$ (76)

Up to a constant factor, this double integral is exactly the same as in the
linearized MFT calculation, see Eq. (45). Equation (76) gives
$\text{Var}_{K}=8\alpha\epsilon$ (where $\alpha=0.05798\dots$), which
perfectly coincides with our MFT result (48).

## References

  * (1) H. Spohn, Large Scale Dynamics of Interacting Particles (Springer, Berlin, 1991).
  * (2) A. De Masi and E. Presutti, Mathematical methods for hydrodynamic limits (Springer, Berlin, 1991).
  * (3) C. Kipnis and C. Landim, Scaling Limits of Interacting Particle Systems (Springer, Berlin, 1999).
  * (4) G. M. Schütz, in Phase Transitions and Critical Phenomena, Vol. 19, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 2001), pp. 1–251.
  * (5) D. Gabrielli and P. L. Krapivsky, J. Stat. Mech 04, 043212 (2018).
  * (6) B. Doyon, SciPost Phys. Lect. Notes 18, 1 (2020).
  * (7) I. Bouchoule and J. Dubail, J. Stat. Mech. 2022, 014003 (2022).
  * (8) B Doyon, J. Math. Phys. 60, 073302 (2019).
  * (9) B. Doyon and J. Myers, Ann. Henri Poincaré 21, 255 (2020).
  * (10) H. Spohn, arXiv:2301.08504 (2023).
  * (11) C. Chamon, Phys. Rev. Lett. 94, 040402 (2005).
  * (12) M. Pretko, Phys. Rev. B 95, 115139 (2017).
  * (13) A. Gromov, Phys. Rev. X 9, 031035 (2019).
  * (14) D. J. Williamson, Z. Bi, and M. Cheng, Phys. Rev. B 100, 125150 (2019).
  * (15) E. Lake, M. Hermele, and T. Senthil, Phys. Rev. B 106, 064511 (2022).
  * (16) B. C. Rayhaun and D. J. Williamson, SciPost Phys. 15, 017 (2023).
  * (17) A. Gromov and L. Radzihovsky, Rev. Mod. Phys. 96, 011001 (2024).
  * (18) A. Gromov, A. Lucas, and R. M. Nandkishore, Phys. Rev. Res. 2, 033124 (2020).
  * (19) J. Feldmeier, P. Sala, G. De Tomasi, F. Pollmann, and M. Knap, Phys. Rev. Lett. 125, 245303 (2020).
  * (20) J. Iaconis, A. Lucas, and R. Nandkishore, Phys. Rev. E 103, 022142 (2021).
  * (21) A. Głódkowski, F. Peña Benítez, and P. Surówka, Phys. Rev. E 107, 034142 (2023).
  * (22) A. Morningstar, Ni. O’Dea, and J. Richter, Phys. Rev. B 108, L020304 (2023).
  * (23) J. H. Han, E. Lake, and S. Ro, Phys. Rev. Lett. 132, 137102 (2024).
  * (24) For a coarse-grained description to apply the particle density must be sufficiently large so as to avoid jamming, see e.g. Ref. Arita .
  * (25) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Rev. Mod. Phys. 87, 593 (2015).
  * (26) W. W. Mullins, J. Appl. Phys. 28, 333 (1957).
  * (27) G. I. Barenblatt, Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics (Cambridge Texts in Applied Mathematics, Series Number 14) (Cambridge University Press, Cambridge, UK, 1996).
  * (28) B. Derrida and A. Gerschenfeld, J. Stat. Phys. 137, 978 (2009).
  * (29) P. L. Krapivsky and B. Meerson, Phys. Rev. E 86, 031106 (2012).
  * (30) B. Meerson and A. Vilenkin, Phys. Rev. E (R) 93, 020102 (2016); Phys. Rev. E 93, 059901 (2016).
  * (31) F. D. Cunden, P. Facchi, and P. Vivo, J. Phys. A 49, 135202 (2016).
  * (32) P. L. Krapivsky, B. Meerson, and P. V. Sasorov, J. Stat. Mech. (2012) P12014.
  * (33) L. D. Landau and E. M. Lifshitz, Mechanics (Reed, Oxford, UK, 1981).
  * (34) C. Arita, P. L. Krapivsky, and K. Mallick, Phys. Rev. E 95, 032121 (2017).


# Propagation of wave packets along large-scale background waves

D. V. Shaykin  Institute of Spectroscopy, Russian Academy of Sciences,
Troitsk, Moscow, 108840, Russia  Moscow Institute of Physics and Technology,
Institutsky lane 9, Dolgoprudny, Moscow region, 141700, Russia  A. M.
Kamchatnov<EMAIL_ADDRESS>Institute of Spectroscopy, Russian Academy
of Sciences, Troitsk, Moscow, 108840, Russia  Moscow Institute of Physics and
Technology, Institutsky lane 9, Dolgoprudny, Moscow region, 141700, Russia

###### Abstract

We study propagation of high-frequency wave packets along a large-scale
background wave which evolves according to dispersionless hydrodynamic
equations for two variables (fluid density and flow velocity). Influence of
the wave packet on evolution of the background wave is neglected, so the
large-scale evolution can be found independently of the wave packet’s motion.
At the same time, propagation of the packet depends in essential way on the
background wave and it can be considered in framework of geometric optics
approximation with the use of Hamilton equations for the carrier wave number
and the mean coordinate of the packet. We derive equations for the carrier
wave number as a function of the parameters which describe the background
wave. When they are solved, the path of the packet can be found by simple
integration of the Hamilton equation. The theory is illustrated by its
application to the problem of propagation of wave packets along expanding
large-scale wave which evolution is described by the shallow water equations.
In particular, they correspond to the dispersionless limit of the defocusing
nonlinear Schrödinger equation, and then the expanding wave can be considered
as an expanding cloud of the Bose-Einstein condensate. Reflection of wave
packets from upstream flows and their propagation along stationary flows are
also discussed. The analytical solutions found for these particular cases
agree very well with exact numerical solution of the nonlinear Schrödinger
equation.

###### pacs:

43.20.Bi, 47.35.+i

##  I Introduction

Propagation of high-frequency wave packets through non-uniform and non-
stationary media is a well-developed field of wave physics. In this theory, it
is usually assumed that there exist three different scales of length: the
wavelength $\lambda=2\pi/k$, corresponding to the modulus $k=|{\bf k}|$ of the
carrier wave vector ${\bf k}$, is much smaller that the size $d$ of the
packet, and $d$ is much smaller than the characteristic length $l$ at which
parameters of the medium considerably change,

$\lambda\ll d\ll l.$ (1)

The difference of scales $d\ll l$ allows one to introduce a “mean” coordinate
${\bf r}$ of the packet with accuracy $d$, and the carrier wave vector ${\bf
k}$ is defined with accuracy $2\pi/d\ll 2\pi/\lambda=k$. Thus, in this
geometrical optics approximation the packet’s dynamics is described by two
variables ${\bf r}$ and ${\bf k}$ which obey the Hamilton equations

$\frac{d{\bf r}}{dt}=\frac{\partial\omega}{\partial{\bf k}},\qquad\frac{d{\bf
k}}{dt}=-\frac{\partial\omega}{\partial{\bf r}},$ (2)

where $\omega=\omega({\bf k},{\bf r},t)$ is the dispersion function for linear
waves propagating through the medium whose parameters have the values
corresponding to the point ${\bf r}$ at the moment $t$. From an elementary
physics point of view, the first equation (2) is just a definition of the
group velocity of a packet and the second equation (2) describes “refraction”
of waves and can be considered as a generalization of Snell’s refraction law.
Coincidence of Eqs. (2) with the Hamilton equations of classical mechanics of
point particles is the essence of the optical-mechanical analogy which, on one
hand, played so important role in development of quantum mechanics and, on the
other hand, allows one to use for solving wave propagation problems the well-
developed methods of analytical mechanics (see, e.g., [1, 2] and references
therein).

In practice, the possibility of analytical treatment of Eqs. (2) depends on
specification of properties of the medium and, consequently, on the form of
the dispersion relation $\omega=\omega({\bf k},{\bf r},t)$. As has recently
been noticed [3], one such a fruitful specification appears when non-
uniformity and time-dependence are caused by a large-scale wave propagating
through the medium. If such a large-scale wave is described by two variables
$\rho({\bf r},t)$ (“density” of the medium) and ${\bf u}({\bf r},t)$ (its
“flow velocity”), then the packet’s frequency $\omega$ depends on ${\bf r}$
and $t$ only via the variables $\rho$ and ${\bf u}$: $\omega=\omega({\bf
k};\rho({\bf r},t),{\bf u}({\bf r},t))$. The problem of propagation of wave
packets was solved in Ref. [3] for a particular case of the large-scale
background wave in the form of a rarefaction wave in a unidirectional small-
amplitude (KdV) approximation with the use of the small-amplitude limit of the
Whitham modulation equations [4, 5] instead of Hamilton’s equations (2). This
approach was generalized in Ref. [6] to the case of general simple waves in
one-dimensional geometry (1D) when $\rho$ and $u$ are functionally related
with each other so that $\omega$ becomes a function of only one medium’s
parameter, $\omega=\omega(k,u(x,t))$, where $u(x,t)$ obeys in the simple wave
case to the Hopf equation

$u_{t}+V_{0}(u)u_{x}=0.$ (3)

Combining Eqs. (2) and (3), one readily gets [6] the differential equation for
the function $k=k(u)$,

$\frac{dk}{du}=\frac{\partial\omega(k,u)/\partial
u}{V_{0}(u)-\partial\omega(k,u)/\partial k}$ (4)

(it was first found in Ref. [7] as a consequence of Whitham’s equations at the
small-amplitude edge of dispersive shock waves). Its solution $k=k(u)$ for a
given initial condition $k=k_{0}$ at $u=u_{0}$ together with the solution of
the Hopf equation (3) can be used for solving the first Hamilton equation (2)
(see Ref. [8]), so that one gets the path $x=x(t)$ of the packet along
evolving background simple wave.

Although simple waves represent an important case of nonlinear waves appearing
in fluid dynamics problems, there are many situations described by general
solutions of compressible fluid dynamics equations, when there is no any
functional relationship between $\rho$ and $u$. For example, such general
solutions appear in problems of expansion of gas clouds, problems of
reflection of nonlinear waves from solid walls, etc. (see, e.g., Ref. [9]). In
this case, the wave number $k$ is affected by two phenomena—refraction due to
change of density and Doppler shift due to change of flow velocity, which are
independent of each other on the contrary to the simple wave case. The
resulting value $k=k(x,t)$ is to be obtained by means of integration of these
changes along the packet’s path which itself depends on the carrier wave
number $k$ as well as on the initial conditions $\rho=\rho_{0}(x),u=u_{0}(x)$
for the background flow. We will look for situations when the wave number $k$
is the same function of the local values of $\rho$ and $u$,
$k=k(\rho(x,t),u(x,t))$, for any choice of the initial conditions for the
background flow. This greatly simplifies the analysis but imposes essential
conditions on the form of the function $k=k(\rho,u)$. Nevertheless, the
inequalities (1) allow one to distinguish situations when such a function
$k=k(\rho,u)$ does exist in the most physically important limit of high-
frequency wave packets, i.e. in the limit of large $k$.

The aim of this paper is to develop the theory of propagation of wave packets
along large-scale background waves represented by general solutions of fluid
dynamics equations. We will formulate the applicability condition of our
approach and will show in a particular case of generalized nonlinear
Schrödinger (gNLS) equation that in the limit of large wave numbers the
function $k=k(\rho,u)$ is given by a simple analytical expression. The theory
is illustrated by several examples of large-scale background flows and
nonlinearity types and our analytical results agree very well with exact
numerical solutions of the gNLS equation for appropriate initial conditions.

##  II General theory

As was noticed in Ref. [10], dynamics of many physical systems depending on
two parameters can be written in hydrodynamics-like form

$\rho_{t}+(\rho u)_{x}=0,\qquad(\rho u)_{t}+(\rho u^{2}+P)_{x}=0,$ (5)

where $\rho$ and $u$ are considered as ‘density’ and ‘flow velocity’ of the
‘fluid’. Consequently, Eqs. (5) can be treated as conservation laws of mass
and momentum, respectively. The function $P=P(\rho,u,\rho_{x},u_{x},...)$
describes the action of pressure as well as effects of dispersion and/or
viscosity which contribution into the momentum density depends on gradients
and higher order derivatives of $\rho$ and $u$. Equations (5) have uniform
solutions $\rho=\overline{\rho}=\mathrm{const}$,
$u=\overline{u}=\mathrm{const}$. If we linearize these equations with respect
to small deviations $\rho^{\prime}=\rho-\overline{\rho}$,
$u^{\prime}=u-\overline{u}$ and look for harmonic wave solutions
$\rho^{\prime},u^{\prime}\propto\exp[i(kx-\omega t)]$, then we get the
dispersion relation for linear waves

$\omega=\omega(k,\overline{\rho},\overline{u}).$ (6)

It is worth noticing that usually one obtains two such expressions
corresponding to linear waves which propagate upstream or downstream the flow
$\overline{u}$, but we do not distinguish them here in our notation. In the
theory of propagation of high-frequency wave packets the variables
$\overline{\rho}$ and $\overline{u}$ can be replaced by the local background
values $\rho$ and $u$ as long as they change slowly compared with the packet’s
size (see Introduction). Of course, this long wave dynamics obeys the
nonlinear equations obtained from Eqs. (5) in dispersionless limit.

Equations for evolution of the background variables are obtained from Eqs. (5)
when we neglect the terms with derivatives in the function $P$ and assume that
the limiting equation of state has the form $P=P(\rho)$, so that the resulting
equations can be written in the form

$\rho_{t}+(\rho u)_{x}=0,\qquad u_{t}+uu_{x}+\frac{c^{2}}{\rho}\rho_{x}=0,$
(7)

where $c=c(\rho)$ denotes the sound velocity ($c^{2}=dP/d\rho$), that is in
the long wavelength limit $k\to 0$ the dispersion relation (6) takes one of
the forms

$\frac{\omega}{k}=v_{+}=u+c\quad\text{or}\quad\frac{\omega}{k}=v_{-}=u-c$ (8)

for sound waves propagating upstream or downstream the flow $u$. The sound
velocity $c=c(\rho)$ depends on the local value of the density $\rho$ and can
also serve as a local wave variable instead of $\rho$.

In the linear limit the variables $u^{\prime}$ and $\rho^{\prime}$ are related
by the formulas $u^{\prime}=\pm(c/\rho)\rho^{\prime}$ for these two ‘right’
and ‘left’ propagating waves with velocities $v_{+}$ and $v_{-}$,
respectively. In a large scale wave the variables $u^{\prime}$ and
$\rho^{\prime}$ can be considered as differentials along small segments of
such a wave, and integration of the equation $du=\pm(c/\rho)d\rho$ yields two
important variables

$r_{\pm}=\frac{u}{2}\pm\frac{1}{2}\int_{0}^{\rho}\frac{cd\rho}{\rho}$ (9)

called Riemann invariants (factor 1/2 is introduced for further convenience).
If we invert the function $c=c(\rho)$ and substitute $\rho=\rho(c)$ into Eq.
(9), then we get the Riemann invariants in the form

$r_{\pm}=\frac{u}{2}\pm\sigma(c),\qquad\sigma(c)=\frac{1}{2}\int_{0}^{c}\frac{c}{\rho(c)}\frac{d\rho}{dc}dc.$
(10)

Evidently, all physical variables $u,\rho,c$ can be expressed in terms of the
Riemann invariants,

$u=r_{+}+r_{-},\quad\rho=\rho(r_{+},r_{-}),\quad c=c(r_{+},r_{-}).$ (11)

Equations (7) written in terms of $r_{\pm}$ take very simple diagonal form

$\frac{\partial r_{+}}{\partial t}+v_{+}\frac{\partial r_{+}}{\partial
x}=0,\quad\frac{\partial r_{-}}{\partial t}+v_{-}\frac{\partial
r_{-}}{\partial x}=0,$ (12)

where the characteristic velocities are equal to

$v_{\pm}=v_{\pm}(r_{+},r_{-})=r_{+}+r_{-}\pm c(r_{+},r_{-}).$ (13)

Simple waves correspond to unidirectional flows with one of the Riemann
invariants constant. Then one of Eqs. (12) is satisfied identically and the
other one reduces to the Hopf equation, so we return to situations considered
in Ref. [6]. In general solutions both Riemann invariants $r_{\pm}$ change
with space coordinate $x$ and time $t$, so they can serve as new coordinates
in the region of the general solution. This hodograph transform makes $x$ and
$t$ functions of $r_{+}$ and $r_{-}$,

$x=x(r_{+},r_{-}),\qquad t=t(r_{+},r_{-}),$ (14)

and it casts Eqs. (12) to linear equations for $x$ and $t$,

$\frac{\partial x}{\partial r_{+}}-v_{-}\frac{\partial t}{\partial
r_{+}}=0,\quad\frac{\partial x}{\partial r_{-}}-v_{+}\frac{\partial
t}{\partial r_{-}}=0,$ (15)

which can be solved, for example, by the hodograph method of Ref. [11]. If the
solution (14) is known, then either it can be inverted to give explicit
formulas for the Riemann invariants

$r_{+}=r_{+}(x,t),\qquad r_{-}=r_{-}(x,t)$ (16)

and, consequently, for the physical variables (11), or it gives these physical
variables in a parametric form.

When a high-frequency wave packet propagates along the background wave given
by the general solution of hydrodynamic equations in the Riemann diagonal form
(12), then both the dispersion relation and the group velocity become
functions of the Riemann invariants,

$\omega=\omega(k,r_{+},r_{-}),\qquad v_{g}=\frac{\partial\omega}{\partial
k}=v_{g}(k,r_{+},r_{-}).$ (17)

Let the packet enter into the region of the general solution at some point
$x_{0}=x(r_{+}^{0},r_{-}^{0})$ at the moment $t_{0}=t(r_{+}^{0},r_{-}^{0})$
with the carrier wave number $k_{0}$. Then its values at later moments of time
depend on the path in the $(r_{+},r_{-})$ plane from the initial point
$(r_{+}^{0},r_{-}^{0})$ to the final point $(r_{+},r_{-})$. Motion of the
packet is governed by the 1D version of Hamilton’s equations (2),

$\frac{dx}{dt}=\frac{\partial\omega}{\partial
k},\qquad\frac{dk}{dt}=-\frac{\partial\omega}{\partial x}.$ (18)

Then the second equation gives

$\frac{dk}{dt}=-\frac{\partial\omega}{\partial r_{+}}\frac{\partial
r_{+}}{\partial x}-\frac{\partial\omega}{\partial r_{-}}\frac{\partial
r_{-}}{\partial x}.$

In Eqs. (17) the frequency $\omega=\omega(k,r_{+},r_{-})$ depends, besides
$k$, only on local values $r_{\pm}$ of the Riemann invariants. We assume here
that the wave number $k$ is also a function of the Riemann invariants,
$k=k(r_{+},r_{-})$. Then along the packet’s path we have

$\begin{split}\frac{dk}{dt}&=\frac{\partial k}{\partial
r_{+}}\frac{dr_{+}}{dt}+\frac{\partial k}{\partial r_{-}}\frac{dr_{-}}{dt}\\\
&=\frac{\partial k}{\partial r_{+}}\left(\frac{\partial r_{+}}{\partial
t}+v_{g}\frac{\partial r_{+}}{\partial x}\right)+\frac{\partial k}{\partial
r_{-}}\left(\frac{\partial r_{-}}{\partial t}+v_{g}\frac{\partial
r_{-}}{\partial x}\right)\\\ &=-\frac{\partial k}{\partial
r_{+}}(v_{+}-v_{g})\frac{\partial r_{+}}{\partial x}-\frac{\partial
k}{\partial r_{-}}(v_{-}-v_{g})\frac{\partial r_{-}}{\partial x},\end{split}$
(19)

where we have used Eqs. (12). Comparison of these two expressions shows that
they are consistent with each other, if the condition

$\begin{split}&\left[\frac{\partial k}{\partial
r_{+}}(v_{+}-v_{g})-\frac{\partial\omega}{\partial r_{+}}\right]\frac{\partial
r_{+}}{\partial x}\\\ &+\left[\frac{\partial k}{\partial
r_{-}}(v_{-}-v_{g})-\frac{\partial\omega}{\partial r_{-}}\right]\frac{\partial
r_{-}}{\partial x}=0\end{split}$ (20)

is fulfilled. As we assumed above, the function $k=k(x,t)$ reduces in
situations under consideration to the function $k=k(r_{+}(x,t),r_{-}(x,t))$
independent of the initial conditions, so the expressions in square brackets
are only functions of $r_{+}$ and $r_{-}$. Let us denote them as
$X_{+}(r_{+},r_{-})$ and $X_{-}(r_{+},r_{-})$. At the same time $r_{+}(x,t)$
and $r_{-}(x,t)$ depend on the initial conditions $r_{+}^{0}(x),r_{-}^{0}(x)$,
so the derivatives $\partial r_{\pm}/\partial x$ at the points with the same
values of $r_{\pm}$ can be different for different choices of the initial
conditions. Let us choose two such initial conditions
$r_{\pm}^{0(1)}(x),r_{\pm}^{0(2)}(x)$, that the corresponding solutions
$r_{\pm}^{(1)}(x,t),r_{\pm}^{(2)}(x,t)$ of Eqs. (12) have at the points with
$r_{\pm}^{(1)}=r_{\pm}^{(2)}\equiv r_{\pm}$ the derivatives with respect to
$x$ satisfying the inequality

$\frac{\partial r_{+}^{(1)}}{\partial x}\cdot\frac{\partial
r_{-}^{(2)}}{\partial x}-\frac{\partial r_{-}^{(1)}}{\partial
x}\cdot\frac{\partial r_{+}^{(2)}}{\partial x}\neq 0.$

Then the system

$\begin{split}&X_{+}(r_{+},r_{-})\frac{\partial r_{+}^{(1)}}{\partial
x}+X_{-}(r_{+},r_{-})\frac{\partial r_{-}^{(1)}}{\partial x}=0,\\\
&X_{+}(r_{+},r_{-})\frac{\partial r_{+}^{(2)}}{\partial
x}+X_{-}(r_{+},r_{-})\frac{\partial r_{-}^{(2)}}{\partial x}=0\end{split}$
(21)

gives $X_{+}(r_{+},r_{-})=0,X_{-}(r_{+},r_{-})=0$, and the function
$k=k(r_{+},r_{-})$ must satisfy the equations

$\frac{\partial k}{\partial r_{+}}=\frac{\partial\omega/\partial
r_{+}}{v_{+}-v_{g}},\qquad\frac{\partial k}{\partial
r_{-}}=\frac{\partial\omega/\partial r_{-}}{v_{-}-v_{g}}.$ (22)

For existence of such a function $k=k(r_{+},r_{-})$ these derivatives must
commute

$\frac{\partial}{\partial r_{+}}\left(\frac{\partial k}{\partial
r_{-}}\right)=\frac{\partial}{\partial r_{-}}\left(\frac{\partial k}{\partial
r_{+}}\right).$ (23)

If Eqs. (22) satisfy this condition, then we look for their solution
$k=k(r_{+},r_{-},q)$, where $q$ is an integration constant which value is
determined by the initial value $k_{0}$ of the wave number at the point
$x_{0}$ and moment of time $t_{0}$ with local values of the Riemann invariants
$r_{\pm}^{0}=r_{\pm}(x_{0},t_{0})$. As a result, the solution
$k=k(r_{+},r_{-},q)$ defines the carrier wave number $k$ in the whole region
of the solution $r_{\pm}=r_{\pm}(x,t)$ of hydrodynamic equations (12).
Substitution of these functions into the first Hamilton equation (18) gives
the equation

$\frac{dx}{dt}=\left.\frac{\partial}{\partial
k}\omega\big{(}k,r_{+}(x,t),r_{-}(x,t)\big{)}\right|_{k=k[r_{+}(x,t),r_{-}(x,t)]}$
(24)

for the packet’s path $x=x(t)$.

If the condition (23) is not fulfilled exactly, then we should confine
ourselves to an approximate solution correct in the limit of large $k$ and
such a solution would be enough for our treatment of propagation of high-
frequency wave packets. To formulate the conditions of applicability of such
an approximation, we assume that the flow and sound velocities have the same
order of magnitude, $|u|\sim c$, and therefore $r_{\pm}\sim c$. Now we look
for an approximate solution for $k\gg c$, where $k\sim k_{0}$, $k_{0}$ being
the initial value of the wave number. If we denote the right-hand sides of
Eqs. (22) as $R_{\pm}(k,r_{+},r_{-})$, then we expand them in series with
respect to small parameters $|r_{\pm}|/k$ and keep only the terms
$R_{\pm}^{\text{asymp}}(k,r_{+},r_{-})$ satisfying the commutativity condition
(23). The solution $k=k^{\text{asymp}}(r_{+},r_{-},q)$ of the resulting
equations

$\frac{\partial k}{\partial r_{\pm}}=R_{\pm}^{\text{asymp}}(k,r_{+},r_{-})$
(25)

should also be looked for with the same accuracy.

If the solution of hydrodynamic equations (7) for the background wave is found
directly for the physical variables $\rho$ and $u$, then it is convenient to
write equations for $k=k(\rho,u)$ without transition to the Riemann invariants
although they may be known. Simple calculation similar to derivation of Eqs.
(22) gives

$\begin{split}&\frac{\partial
k}{\partial\rho}=\frac{(v_{g}-u)\frac{\partial\omega}{\partial\rho}+\frac{c^{2}}{\rho}\frac{\partial\omega}{\partial
u}}{c^{2}-(v_{g}-u)^{2}},\\\ &\frac{\partial k}{\partial
u}=\frac{(v_{g}-u)\frac{\partial\omega}{\partial
u}+\rho\frac{\partial\omega}{\partial\rho}}{c^{2}-(v_{g}-u)^{2}}.\end{split}$
(26)

Naturally, these derivatives should satisfy the commutativity condition
$\frac{\partial}{\partial u}\left(\frac{\partial
k}{\partial\rho}\right)=\frac{\partial}{\partial\rho}\left(\frac{\partial
k}{\partial u}\right)$ at least in the limit of large $k$, when the asymptotic
solution $k=k^{\text{asymp}}(\rho,u,q)$ can be found.

Let us illustrate the general theory by its application to systems whose
evolution is governed by the gNLS equation.

##  III Application to gNLS equation

Here we shall apply the developed above approach to the gNLS equation

$i\psi_{t}+\frac{1}{2}\psi_{xx}-f(|\psi|^{2})\psi=0,$ (27)

which has a number of physical applications and is written here in standard
dimensionless variables. To be definite, we shall imply here that this
equation describes dynamics of Bose-Einstein condensates (BECs) with different
kinds of repulsive interaction between atoms (or quasiparticles in case of
polariton condensates). To cast this equation to the form (5), we make the
Madelung transformation of the BEC’s wave function

$\psi(x,t)=\sqrt{\rho(x,t)}\exp\left(i\int^{x}u(q,t)dq\right)$ (28)

and obtain the system

$\begin{split}&\rho_{t}+(\rho u)_{x}=0,\\\
&u_{t}+uu_{x}+\frac{c^{2}}{\rho}\rho_{x}+\bigg{(}\frac{\rho_{x}^{2}}{8\rho^{2}}-\frac{\rho_{xx}}{4\rho}\bigg{)}_{x}=0,\end{split}$
(29)

where

$c^{2}=\rho f^{\prime}(\rho)$ (30)

is a function of the condensate’s density $\rho=|\psi|^{2}$ and $u$ has a
meaning of its flow velocity. The last term in the second equation (29) is
responsible for the dispersive effects, so linearization of these equations
with respect to small deviations from a uniform flow with constant density and
flow velocity yields the Bogoliubov dispersion relation

$\omega=k\left(u\pm\sqrt{c^{2}+\frac{k^{2}}{4}}\right).$ (31)

In the limit $k\to 0$ we obtain the characteristic velocities (8),
consequently $c$ equals to the sound velocity of long waves. Dynamics of a
large scale wave obeys the dispersionless ‘shallow water’ equations (7) which
are obtained from (29) by neglecting the dispersion terms.

Formula (30) determines $c$ as a function of $\rho$, $c=c(\rho)$. It is
convenient to invert this function and then the function $\rho=\rho(c)$
characterizes the nonlinear properties of the system under consideration.
Thus, in this notation the large scale flow is described by the variables
$c(x,t),u(x,t)$ which are equivalent to the previous variables
$\rho(x,t),u(x,t)$ but they are more convenient in some calculations. Then we
easily obtain

$\begin{split}&v_{g}=\frac{\partial\omega}{\partial
k}=u+\frac{k^{2}+2c^{2}}{\sqrt{k^{2}+4c^{2}}},\qquad\frac{\partial\omega}{\partial
u}=k,\\\
&\frac{\partial\omega}{\partial\rho}=\frac{2ck/\rho^{\prime}}{\sqrt{k^{2}+4c^{2}}},\end{split}$
(32)

and substitution of these formulas into Eqs. (26) casts them to the form

$\begin{split}&\frac{\partial k}{\partial
c}=-\frac{c[(2+c\rho^{\prime}/\rho)k^{2}+4(1+c\rho^{\prime}/\rho)c^{2}]}{k(k^{2}+3c^{2})},\\\
&\frac{\partial k}{\partial
u}=-\frac{\sqrt{k^{2}+4c^{2}}[k^{2}+2(1+\rho/(c\rho^{\prime}))c^{2}]}{k(k^{2}+3c^{2})}.\end{split}$
(33)

The difference of cross-derivatives equals to

$\begin{split}&\frac{\partial}{\partial c}\left(\frac{\partial k}{\partial
u}\right)-\frac{\partial}{\partial u}\left(\frac{\partial k}{\partial
c}\right)=\frac{\sqrt{k^{2}+4c^{2}}}{k\rho\rho^{\prime}(k^{2}+3c^{2})^{2}}\\\
&\times\left[(k^{2}+6c^{2})\rho^{\prime
2}(\rho^{\prime}-2c)+2(k^{2}+3c^{2})\rho^{2}(c\rho^{\prime\prime}-\rho^{\prime})\right].\end{split}$
(34)

As one can see, this difference only vanishes in the case of $\rho=c^{2}$
which corresponds to the standard NLS equation (27) with $f(\rho)=\rho$. In
the limit of large $k$ this difference vanishes as $\propto k^{-2}$. Thus, in
the limit of large $k$ the condition (23) is fulfilled and the wave number $k$
becomes a function of two variables, $k=k(c,u)$. In this limit the
corresponding partial derivatives are given by the asymptotic expressions

$\frac{\partial k^{2}}{\partial
c^{2}}=-2\left(1+\frac{c^{2}}{\rho}\frac{d\rho}{dc^{2}}\right),\qquad\frac{\partial
k}{\partial u}=-1,$ (35)

obtained from Eqs. (33) for large $k$. To solve these equations in the
asymptotic limit, we notice that the second equation gives $k=q-u+F(c^{2})$,
where in the main approximation $k\approx q\sim k_{0}\gg|u|,c$. We suppose
that $F(c^{2})\sim c^{2}$ (this will be confirmed by the final result), so
$k^{2}\approx(q-u)^{2}+2qF(c^{2})$, where we have neglected small terms
$uF\sim c^{3}$ and $F^{2}\sim c^{4}$. Then the first equation (35) gives at
once

$F(c^{2})=-\frac{1}{q}\left(c^{2}+\int_{0}^{\rho(c)}\frac{c^{2}}{\rho}d\rho\right)\sim\frac{c^{2}}{q}$

and we arrive at the asymptotic solution

$k^{2}=(q-u)^{2}-2\left(c^{2}+\int_{0}^{\rho(c)}\frac{c^{2}}{\rho}d\rho\right).$
(36)

Here $q$ is an integration constant determined by the value of $k=k_{0}$ at
some initial point in the $(c,u)$-plane and, hence, $q=k_{0}\gg|u|,c$ in this
asymptotic solution. For an important case of the nonlinearity function

$f(\rho)=\frac{1}{\gamma-1}\rho^{\gamma-1},\quad\text{when}\quad
c^{2}=\rho^{\gamma-1},\quad\rho=c^{2/(\gamma-1)},$ (37)

we get

$k^{2}=(q-u)^{2}-\frac{2\gamma}{\gamma-1}c^{2},$ (38)

where $\gamma$ is an effective ‘adiabatic constant’ (in gNLS equation case it
does not have a meaning of ratio of specific heat capacities, of course). If
$\gamma=2$, then the formula

$k^{2}=(q-u)^{2}-4c^{2}$ (39)

gives the exact solution of Eqs. (33) with $\rho=c^{2}$, when these
derivatives commute. Expression (36) for the function $k=k(c,u)$ and its
particular cases (38) and (39) can be used for solving various problems on
propagation of high-frequency wave packets.

##  IV Propagation along expanding BEC

###  IV.1 NLS equation $(\gamma=2)$

At first we shall consider the case of the standard NLS equation with the
nonlinearity function $f(\rho)=\rho$, $(\gamma=2)$. Let at the initial moment
of time the BEC be confined in a standard trap with harmonic potential, so in
the Thomas-Fermi approximation distribution of the density has a form of
inverted parabola, and $u$ is equal to zero everywhere

$\begin{cases}\rho(x,0)=a^{2}\big{[}1-\big{(}\frac{x}{l}\big{)}^{2}\big{]},\quad|x|\leq
l,\\\ u(x,0)=0,\end{cases}$ (40)

where $l\gg 1$ is the size of the BEC cloud and $a^{2}$ is the maximal density
at the center $x=0$ of the trap. At the moment $t=0$ the trap is switched off
and the condensate starts its expansion. Equations (7) with the initial
conditions (40) were solved in Ref. [14] with the use of ansatz suggested in
Ref. [15] for description of self-focusing of light beams propagating through
nonlinear medium. The exact solution has the form

$\begin{cases}\rho(x,t)=\frac{a^{2}}{f(t)}\big{[}1-\big{(}\frac{x}{lf(t)}\big{)}^{2}\big{]},\quad|x|\leq
lf(t),\\\ u(x,t)=\phi(t)x,\end{cases}$ (41)

where $f(0)=1,\phi(0)=0$ according to the initial conditions (40).
Substitution of these expressions into (7) leads to the system for the
functions $f(t),\phi(t)$:

$\phi=\frac{1}{f}\frac{df}{dt},\qquad\frac{d^{2}f}{dt^{2}}=\frac{2a^{2}}{l^{2}}\frac{1}{f^{2}},$
(42)

These equations can be readily solved to give the function $f(t)$ in an
implicit form

$\frac{2a}{l}t=\sqrt{f(f-1)}+\ln\left(\sqrt{f-1}+\sqrt{f}\right),$ (43)

and then the first expression in Eqs. (42) defines the function $\phi(t)$.
Equations (41)-(43) describe in a parametric form the large scale evolution of
the condensate cloud.

(a)

(b)

Figure 1: (a) Profiles $\rho(x,t)$ of the density obtained by numerical
solution of the Gross-Pitaevskii equation with the initial distribution (40)
for $a=5$, $l=35$, $k_{0}=5$ and the wave packet added at the center; the
curves correspond to the moments of time $t=0$ (red), $t=35$ (blue). (b)
Profiles of density in the linear wave packet obtained by subtraction of the
result of evolution of the smooth background pulse from the evolution of the
whole condensate with account of contribution of the wave packet at the same
moments of time.

Now we assume that a high-frequency disturbance of condensate’s density is
located in vicinity of its center at the moment $t=0$ of switching off the
trap. This disturbance can be represented as a wave packet made from the
Bogoliubov harmonics with both signs in Eq. (31) and some carrier wave vector
$k_{0}$. Consequently, this disturbance splits to two wave packets propagating
symmetrically along the expanding cloud of BEC (see Fig. 1) and the carrier
wave vector evolves according to the equations derived in Section III. As was
noticed above, in this case Eqs. (33) take simple form

$\begin{split}\frac{\partial k}{\partial c}=-4\frac{c}{k},\qquad\frac{\partial
k}{\partial u}=-\frac{\sqrt{k^{2}+4c^{2}}}{k},\end{split}$ (44)

and these derivatives commute with each other. It is remarkable that this case
corresponds to the complete integrability of the NLS equation, so one may
suppose that these two properties can be related with each other. The exact
solution of these equations coincides with the solution (39) of their large
$k$ approximations. Equation (24) for the packet’s path takes with account of
Eqs. (31), (39) the form (we assume here $q-u>0$)

$\frac{dx}{dt}=q-\frac{2\rho(x,t)}{q-u(x,t)},$ (45)

and it can be easily solved numerically for the given distributions (41)-(43)
of $\rho(x,t)$ and $u(x,t)$. We compared this theory with the exact numerical
solution of the Gross-Pitaevskii equation. To this end, we have chosen the
initial distributions (40) with the parameters $a=5$, $l=35$ and added at the
center $x=0$ a small perturbation in the form of the wave packet with
$k_{0}=5$ or $k_{0}=10$. In numerical solution of Eq. (27) the initial
disturbance splits into two packets (see Fig. 1) and we will consider the
right-propagating packet. Since splitting of the initial wave packet to the
right and left propagating packets takes some time, their coordinates may be
not well enough defined at small values of time and in these cases the
numerical corresponding points are not depicted in our figures. The paths of
the wave packet for these two initial values of $k_{0}$ are shown in Fig. 2,
where the solid line depicts the solution of Eq. (45) and the dots correspond
to locations of the packet extracted from the full numerical solution of the
Gross-Pitaevskii equation (27). As one can see, our analytical approach agrees
very well with the exact numerical solution even for relatively small values
of $k$.

Figure 2: Paths $x(t)$ of the wave packets along expanding BEC $(\gamma=2)$
for two values of the initial carrier wave vector $k_{0}=5$ or $k_{0}=10$.

###  IV.2 gNLS equation $(\gamma=3)$

Now we shall consider the case of the gNLS equation (27) with
$f(\rho)=\rho^{2}/2$, so that $c=\rho$ (see Eq. (37)). In this case it is
convenient to redefine the Riemann invariants (10) by multiplying them by 2 so
we get $r_{\pm}=v_{\pm}=u\pm c$ and

$u=\frac{1}{2}(r_{+}+r_{-}),\quad c=\rho=\frac{1}{2}(r_{+}-r_{-}).$ (46)

As a result, the hydrodynamic equations (12) for the background flow split
into two independent Hopf equations

$\frac{\partial r_{\pm}}{\partial t}+r_{\pm}\frac{\partial r_{\pm}}{\partial
x}=0.$ (47)

Such a splitting to two separate Hopf equations can be related with a specific
symmetry of the gNLS equation for the case $\gamma=3$. This symmetry was first
found in Ref. [16] for NLS with two spatial dimensions and generalized in Ref.
[17] to gNLS with any number of spatial dimensions so that $\gamma=3$
corresponds to our case with one spatial dimension. If in the initial state
the flow velocity equals everywhere to zero, then the initial distributions of
the Riemann invariants only differ by signs, $r_{+}(x,0)=c(x,0)$,
$r_{-}(x,0)=-c(x,0)$, so the solutions of the Hopf equations are given by the
formulas

$x-r_{+}t=\overline{x}(r_{+}),\qquad x-r_{-}t=\overline{x}(-r_{-}),$ (48)

where $\overline{x}(r_{+})$ is the function inverse to $r_{+}=r_{+}(x,0)$.

For simplicity we take again the initial distribution of $\rho=c$ in a
parabolic form (40), so that $\overline{x}(r)=l\sqrt{1-r/a^{2}}$. Then the
solution (48) can be written as

$\begin{split}&x-(u+c)t=l\sqrt{1-(u+c)/a^{2}},\\\
&x-(u-c)t=l\sqrt{1+(u-c)/a^{2}},\end{split}$ (49)

and the time-dependent distributions $c=c(x,t),u=u(x,t)$ can be easily
expressed by relatively simple analytical formulas

$\begin{split}\rho(x,t)&=c(x,t)=\frac{l}{4a^{2}t^{2}}\sqrt{l^{2}+4a^{4}l^{2}+4a^{2}tx}+\\\
&+\frac{l}{4a^{2}t^{2}}\big{(}-2l+\sqrt{l^{2}+4a^{4}l^{2}-4a^{2}tx}\big{)},\end{split}$
(50)
$\begin{split}u(x,t)&=\frac{x}{t}+\frac{l}{4a^{2}t^{2}}\sqrt{l^{2}+4a^{4}l^{2}-4a^{2}tx}-\\\
&-\frac{l}{4a^{2}t^{2}}\sqrt{l^{2}+4a^{4}l^{2}+4a^{2}tx}\end{split}$ (51)

Equation (24) takes (with account of $k^{2}=(q-u)^{2}-3c^{2}$, see Eq. (38))
the form

$\frac{dx}{dt}=u+\frac{(q-u)^{2}-c^{2}}{\sqrt{(q-u)^{2}+c^{2}}},$ (52)

where $u=u(x,t)$ and $c=c(x,t)$ are known functions. Again we put a wave
packet at the top of the density distribution with some initial value of the
carrier wave number $k_{0}$, which determines the value of $q$, and solve Eq.
(52) numerically. The results are compared with the exact solution of the full
gNLS equation, see Fig. 3, and one can see that the approximate theory using
the asymptotic formula (38) agrees very well with the exact solutions.

Figure 3: Paths $x(t)$ of the wave packets along expanding BEC $(\gamma=3)$
for two values of the initial carrier wave vector $k_{0}=9$ or $k_{0}=25$.

##  V Packet’s reflection from an upstream flow

So far we have considered propagation of a packet along a single background
wave with coinciding directions of the background flow and the packet’s
propagation. It is of considerable interest to discuss situations when the
packet with the fixed value of $q$ propagates through the region of the
background flow where the group velocity can vanish. In the optical-mechanical
analogy this corresponds to the ‘turning points’ where the packet changes
direction of its propagation. Due to presence of the flow velocity this does
not mean that the wave vector vanishes, so the asymptotic theory can be still
applicable. In this section, we will check the validity of Eq. (39)
($\gamma=2$) with the same value of $q$ after such a reflection of a wave
packet from the turning point.

Figure 4: Initial distributions of two BEC clouds in the form of (40) with
parameters indicated in the picture. The wave packet at $t=0$ is drawn in
brown, it splits for $t>0$ to two modes, one of which propagates to the left,
the other to the right. Our aim is to investigate the left moving packet with
$q=-\sqrt{14^{2}+4}$.

To this end, we launch a relatively slow packet towards a large expanding
condensate in the following way. At the initial moment of time we form two
clouds with parabolic distributions of density: a large cloud occupies the
region with $l=100$, its right edge is located at $x=-20$, and the field
amplitude equals to $a=30$, a small cloud occupies the region $-10<x<10$, so
its left edge is located at $x=-10$, and the amplitude equals to $a=1$. The
wave packet is generated initially at the top of the small cloud around the
point $x=0$ with the wave number $k_{0}=14$, see Fig. 4. Then the clouds start
their expansions, the initial packet splits to two packets, and the left
packet with $q=-\sqrt{14^{2}+4}\approx-14.14$ propagates first towards the
empty region between two clouds, propagates through it, and reaches finally
the large cloud with upstream background flow. In this case the first Hamilton
equation (with sign “–” in the dispersion relation and $q-u<0$) has the same
form as (45). Interestingly enough, propagation through the empty region does
not look in numerical calculations like motion of a well-defined packet: it
disperses at the left edge of the small cloud into a wide distribution of
small-amplitude waves which are collected again into a narrow packet at the
right edge of the large cloud. The essential point is that at the edges of the
both clouds, where $\rho=0$, we find from Eqs. (45) that the following
estimate holds at the edges of the clouds:

$\left.v_{g}\right|_{\text{edges}}=q.$ (53)

This allows us to find analytically trajectories of the packet along three
characteristic regions: the small cloud, the empty region, and the large
cloud; they are shown in Fig. 5 by a solid line. Positions of the packet at
different moments of time are shown by dots; as was indicated above, they
cannot be identified within the empty region. The packet enters the big cloud
with negative group velocity, at some point the group velocity vanishes due to
increase of the upstream background flow velocity, and later vanishes again
because the flow velocity is small near the center of the cloud. The
boundaries of the gray region are obtained by equating the right-hand side of
Eq. (45) to zero for given distributions of $\rho(x,t)$ and $u(x,t)$. Fig. 5
shows good agreement of our analytical approach with numerical solutions of
the NLS equation.

Figure 5: Comparison of the trajectories obtained by numerical solution of the
NLS equation (circles) and by solutions of Eq. 45. The dashed line inside the
ellipse shows the transition through the empty region with group velocity
$v_{g}=q$. In the grey and white areas the packets with the fixed value of
$q=-\sqrt{14^{2}+4}$ have positive and negative velocities, respectively.

##  VI Conservation of wave action

The example of the preceding section demonstrates that though the packet
undergoes drastic transformations during its propagation (in particular, its
amplitude becomes extremely small in the region between the two large scale
flows), nevertheless the value of the constant $q$ in Eq. (39) remains the
same at all stages of evolution. So one can predict the value of the carrier
wave number $k$ at any moment of time.

Equation (39) does not depend on the amplitude of the packet. Therefore it
would be desirable to have another conservation law which could provide some
information about the packet’s amplitude. As is known, the packet’s energy is
not conserved in situations with time-dependent external parameters of the
large scale flow. Instead, as was shown by Whitham [18, 4], the wave action is
preserved provided the external parameters change slowly enough. Here we want
to demonstrate that the packet’s wave action remains constant in spite of
drastic packet’s transformations shown the preceding section.

Figure 6: Blue line denotes analytical prediction for the time dependence of
the carrier wave number. Circles correspond to numerical solutions of the NLS
equation. Black dashed line denotes transition between two expanding BEC
clouds.

At first we shall derive by the Whitham method [18, 4] the expression for the
wave action in case of gNLS equation (27). Linearization of the hydrodynamic-
like equations (29) with respect to small packet’s variables
$\rho^{\prime}=\rho-\overline{\rho},u^{\prime}=u-\overline{u}$, where
$\overline{\rho},\overline{u}$ obey the dispersionless equations (7), yields

$\rho^{\prime}_{t}+\overline{u}\rho^{\prime}_{x}+\overline{\rho}\phi_{xx}=0,\qquad\phi_{t}+\overline{u}\phi_{x}+\frac{c^{2}}{\overline{\rho}}\rho^{\prime}-\frac{\rho^{\prime}_{xx}}{4\overline{\rho}}=0,$
(54)

where we have also introduced the potential $\phi$ of the packet’s
contribution to the flow velocity $u^{\prime}=\phi_{x}$. These equations can
be formally derived from the principle of minimal action with the Lagrangian
density

$\mathop{\mathcal{L}}=\phi_{t}\rho^{\prime}+\overline{u}\rho^{\prime}\phi_{x}+\frac{c^{2}}{2\overline{\rho}}\rho^{\prime
2}+\frac{\overline{\rho}}{2}\phi_{x}^{2}+\frac{(\rho^{\prime}_{x})^{2}}{8\overline{\rho}},$
(55)

where $\overline{\rho},\overline{u}$ are known functions and they should not
be varied in our variational problem (actually they can be considered constant
at the packet’s length scale). Now we define the packet’s envelope variables
by the formulas

$\rho^{\prime}=A(x,t)\cos{\theta(x,t)},\quad\phi=\frac{B(x,t)}{\theta_{x}}\sin{\theta(x,t)},$
(56)

where $\theta(x,t)$ is the high-frequency phase with slowly changing wave
vector $k=\theta_{x}$ and frequency $\omega=-\theta_{t}$. Substitution of Eqs.
(56) into first Eq. (54) gives the relationship

$B=\frac{A}{\overline{\rho}k}\big{(}\omega-\overline{u}k\big{)}.$

Consequently, the Lagrangian density becomes proportional to $A^{2}$ and
averaging of $\mathcal{L}$ over fast oscillations at the packet’s wavelength
scale yields

$\overline{\mathop{\mathcal{L}}}=\frac{A^{2}}{\overline{\rho}}\left\\{c^{2}+\frac{1}{4}\theta_{x}^{2}-\left(\frac{\theta_{t}}{\theta_{x}}+\overline{u}\right)^{2}\right\\},$
(57)

where we have omitted an inessential numerical factor 1/4.

Figure 7: Dependence of the packet’s wave action (63) during propagation of
the packet along the path shown in Fig. 5.

According to Whitham, equations for the envelope variables are to be obtained
from the variational problem for the averaged action

$S=\int\overline{\mathop{\mathcal{L}}}\,dxdt.$ (58)

Variation of $A$ gives

$c^{2}+\frac{k^{2}}{4}=\left(\frac{\omega}{k}-\overline{u}\right)^{2}$ (59)

and this equation reproduces the dispersion relation (31). Variation of
$\theta$ leads to the conservation law

$\frac{\partial}{\partial
t}\frac{\partial\mathop{\mathcal{L}}}{\partial\theta_{t}}+\frac{\partial}{\partial
x}\frac{\partial\mathop{\mathcal{L}}}{\partial\theta_{x}}=0,$ (60)

where

$\begin{split}&\frac{\partial\mathop{\mathcal{L}}}{\partial\theta_{t}}=\frac{2A^{2}}{\overline{\rho}k}\left(\frac{\omega}{k}-\overline{u}\right)=\frac{A^{2}}{\overline{\rho}k}\sqrt{k^{2}+4c^{2}}\equiv
w,\\\
&\frac{\partial\mathop{\mathcal{L}}}{\partial\theta_{x}}=\frac{A^{2}}{\overline{\rho}k}\sqrt{k^{2}+4c^{2}}\left(\overline{u}+\frac{k^{2}+2c^{2}}{\sqrt{k^{2}+4c^{2}}}\right)\equiv
v_{g}w.\end{split}$ (61)

Thus, we arrive at the wave action conservation law

$\frac{\partial w}{\partial t}+\frac{\partial(v_{g}w)}{\partial x}=0,$ (62)

where $w$ is defined in Eq. (61). If the packet is localized, then its total
action

$\mathcal{A}=\int\limits_{-\infty}^{+\infty}\frac{A^{2}}{\overline{\rho}k}\sqrt{4c^{2}+k^{2}}\,dx\int\limits_{-\infty}^{+\infty}\frac{A^{2}}{\overline{\rho}k}|q-\overline{u}|dx$
(63)

should be preserved during propagation through a slowly changing background
(we have used Eq. (39) in the last expression for the wave action for the case
$\gamma=2$). The wave action $\mathcal{A}$ depends on the amplitude $A$ of the
packet and can be used for estimates of $A$, if the packet’s size can be
measured independently.

We have checked conservation of the wave action (63) during propagation of the
wave packet along the path discussed in the preceding section (see Figs. 4-6).
As one can see in Fig. 7, the wave action calculated numerically keeps the
constant value in spite of quite steep turns in the packet’s trajectory.

##  VII Stationary background flow

Our approach assumes that the background flow is described by Eqs. (7) which
do not include external forces: if we add a term $-U_{x}$ corresponding to the
external potential $U(x)$ acting on the gas, then we cannot derive Eqs. (22)
or (26) anymore, so $k$ cannot be a function only of $\rho$ and $u$. However,
in an important particular case of a stationary background flow described by
the equations

$(\rho u)_{x}=0,\qquad uu_{x}+\frac{c^{2}}{\rho}\rho_{x}=-U_{x},$ (64)

Eqs. (22) or (26) are not necessary for finding a path of the high-frequency
wave packet. Indeed, in this case the distributions $\rho=\rho(x)$, $u=u(x)$
do not depend on time, consequently the frequency $\omega=\omega(k,\rho,u)$
obtained from linearized equations does not depend on time either, in
agreement with Hamilton equations (18), i.e., $\omega$ is their ‘energy
integral’. This means that the dependence $k=k(x)$ can be found from the
frequency conservation law

$\omega(k,\rho(x),u(x))=\omega_{0}=\mathrm{const},$ (65)

and then the packet’s path can be found by integration of the equation

$\frac{dx}{dt}=\left.\frac{\partial}{\partial
k}\omega(k,\rho(x),u(x))\right|_{k=k(x)}$ (66)

with a given initial condition.

Figure 8: Wave packet’s path along a stationary flow.

We illustrate here this situation by a simple example for the NLS equation
case with $c^{2}=\rho$ ($\gamma=2$). Equations (64) can be integrated once to
give

$\rho
u=\rho_{0}u_{0},\qquad\frac{1}{2}u^{2}+\rho+U=\frac{1}{2}u_{0}^{2}+\rho_{0},$
(67)

where $\rho_{0},u_{0}$ are the flow variables far from the obstacle whose
action is described by the potential $U(x)$, $U\to 0$ as $|x|\to\infty$.
Elimination of $\rho=\rho_{0}u_{0}/u$ yields the equation

$\frac{1}{2}\left(u_{0}^{2}-u^{2}(x)\right)+\rho_{0}\left(1-\frac{u_{0}}{u(x)}\right)=U(x)$
(68)

which defines the function $u(x)$ for all values of $x$, $-\infty<x<\infty$,
if the maximal value $U_{m}>0$ of the potential $U(x)$ satisfies the condition
(see [19])

$U_{m}\leq\frac{1}{2}u_{0}^{2}-\frac{3}{2}(\rho_{0}u_{0})^{2/3}+\rho_{0}.$
(69)

This inequality becomes equality for two values of $u_{0}$, $u_{0}=u_{\pm}$,
$u_{-}<u_{+}$, so the smooth solutions of Eq. (68) exist for subcritical
$u_{0}<u_{-}$ and supercritical $u_{0}>u_{+}$ flows.

Figure 9: Dependence of the carrier wave number during packet’s propagation
along a stationary flow.

Let the potential be given by the expression

$U(x)=\frac{U_{m}}{\cosh(x/\sigma)}$ (70)

with $U_{m}=15,\sigma=15$ and the flow at infinity have the values of the
parameters $\rho_{0}=1,u_{0}=8$. We launch the high-frequency wave packet at
$x_{0}=50$ with the initial carrier wave number $k_{0}=12$. The functions
$u(x)$ and $\rho(x)=\rho_{0}u_{0}/u(x)$ are found by numerical solution of Eq.
(68) and then the function $k(x)$ is defined by Eq. (65). The path $x=x(t)$
obtained by solving Eq. (66) is shown in Fig. 8 by a solid line and it agrees
very well with the exact numerical solution of the NLS equation whose results
are shown by dots. When $x=x(t)$ is known, then the time dependence of the
carrier wave number $k=k(x(t))$ is known, too; it is shown in Fig. 9 and also
agrees with the numerical solution of the NLS equation. The final value
$k=k_{f}$ at $t\to\infty$ can be predicted from the conservation law (65),

$k_{f}\left(u_{0}-\sqrt{\rho_{0}+\frac{k_{f}^{2}}{4}}\right)=k_{0}\left(u_{0}-\sqrt{\rho_{0}+\frac{k_{0}^{2}}{4}}\right),$
(71)

where the sign before the square root is chosen in such a way that the initial
group velocity (17) is negative. Then we obtain $k_{f}\approx 3.988$, that is
the carrier wave vector remains positive during propagation along this
stationary flow, but the group velocity changes sign at the turning point, so
the packet changes direction of its propagation (see Fig. 8).

##  VIII Conclusion

In this paper, we have extended the approach proposed in Ref. [6] for
description of propagation of high-frequency wave packets along a large scale
simple wave to motion of such wave packets along large-scale background pulses
described by general solutions of hydrodynamic type systems for two variables
$\rho$ and $u$. Due to large difference of scales inherent in our task, the
wave packet’s motion is separated from the background pulse evolution.
Consequently, the background wave variables $\rho,u$ evolve according to
dispersionless equations, whereas propagation of the wave packet is governed
by the Hamilton equations with the carrier wave frequency $\omega$ and wave
vector $k$ playing the roles of the Hamiltonian and canonical momentum,
respectively. We have formulated conditions under which fulfillment the wave
vector $k$ is a function of two variables $\rho,u$. Typically, these
conditions are fulfilled for high-frequency packets in the limit of large $k$.
In this case, the combined system of Hamilton equations and dispersionless
hydrodynamic equations can be reduced to equations for $k$ as a function of
the background wave variables $\rho,u$. Asymptotic solution of these equations
for the gNLS equation case has a quite simple general form (38) and it allows
one to find the packet’s path and the dependence $k(t)$ of the carrier wave
vector on time. Application of this theory to concrete examples demonstrated
its good agreement with the exact numerical solution of the full system.

###### Acknowledgements.

We thank S. K. Ivanov for useful discussions. This research was funded by the
research project FFUU-2021-0003 of the Institute of Spectroscopy of the
Russian Academy of Sciences (sections I-IV) and by the Foundation for the
Advancement of Theoretical Physics and Mathematics “BASIS” (sections V-VII).

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] J. L. Synge, Geometrical Optics. An Introduction into Hamilton’s method, (Cambridge University Press, 1937).
  * [2] Yu. A. Kravtsov, Yu. I. Orlov, Geometrical Optics of Inhomogeneous Media, (Springer, Berlin, 1990).
  * [3] T. Congy, G. A. El, M. A. Hoefer, Interaction of linear modulated waves and unsteady dispersive hydrodynamic states with application to shallow water waves, J. Fluid Mech. 875, 1145-1174 (2019).
  * [4] G. B. Whitham, Non-linear dispersive waves. Proc. R. Soc. Lond. A 283, 238–261 (1965).
  * [5] G. B. Whitham, Linear and Nonlinear Waves, (Wiley Interscience, New York, 1974).
  * [6] A. M. Kamchatnov and D. V. Shaykin, Propagation of wave packets along intensive simple waves, Phys. Fluids, 33, 052120 (2021).
  * [7] G. A. El, Resolution of a shock in hyperbolic systems modified by weak dispersion, Chaos, 15, 037103 (2005).
  * [8] A. M. Kamchatnov, Dispersive shock wave theory for nonintegrable equations, Phys. Rev. E 99, 012203 (2019).
  * [9] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics, Vol. 6: Fluid Mechanics, (Pergamon, New York, 1987).
  * [10] C. H. Su and C. S. Gardner, Korteweg–de Vries equation and generalizations. III. Derivation of the Korteweg–de Vries equation and Burgers equation, J. Math. Phys., 10, 536 (1969).
  * [11] S. P. Tsarev, Geometry of Hamiltonian systems of hydrodynamic type. Generalized hodograph method, Math. USSR-Izvestiya, 37 (2), 397 (1991).
  * [12] B. L. Rozdestvenskii, N. N. Yanenko, Systems of quasilinear equations and their applications to gas dynamics, (AMS, Transl. Math. Monographs; Vol. 55) (Providence, Rhode Island, 1983).
  * [13] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation, (Clarendon Press, Oxford, 2003).
  * [14] V. A. Brazhnyi, A. M. Kamchatnov, V. V. Konotop, Hydrodynamic flow of expanding Bose-Einstein condensates, Phys. Rev., A 68, 035603 (2003).
  * [15] V. I. Talanov, Self-focusing of light beams in nonlinear media, Pisma Zh. Eksp. Teor. Fiz. 2, 218 (1965) [JETP Lett. 2, 138 (1965)].
  * [16] V. I. Talanov, Focusing of light in cubic media, Pisma Zh. Eksp. Teor. Fiz., 11, 303 (1970) [JETP Lett., 11, 199 (1970)].
  * [17] E. A. Kuznetsov, S. K. Turitsyn, Talanov transformations in self-focusing problems and instability of stationary waveguides, Phys. Lett. A 112, 273-275 (1985).
  * [18] G. B. Whitham, A general approach to linear and non-linear dispersive waves using a Lagrangian. J. Fluid Mech., 22, 273-283 (1965).
  * [19] A. M. Leszczyszyn, G. A. El, Yu. G. Gladush, and A. M. Kamchatnov, Transcritical flow of a Bose-Einstein condensate through a penetrable barrier, Phys. Rev., A 79, 063608 (2009).