# Whitham modulation theory for the defocusing nonlinear Schrödinger equation
in two and three spatial dimensions

Asela Abeya1, Gino Biondini1 and Mark A. Hoefer2 1 State University of New
York, Department of Mathematics, Buffalo, NY 14260, USA  2 University of
Colorado, Department of Applied Mathematics, Boulder, CO 80309, USA

###### Abstract

The Whitham modulation equations for the defocusing nonlinear Schrödinger
(NLS) equation in two, three and higher spatial dimensions are derived using a
two-phase ansatz for the periodic traveling wave solutions and by period-
averaging the conservation laws of the NLS equation. The resulting Whitham
modulation equations are written in vector form, which allows one to show that
they preserve the rotational invariance of the NLS equation, as well as the
invariance with respect to scaling and Galilean transformations, and to
immediately generalize the calculations from two spatial dimensions to three.
The transformation to Riemann-type variables is described in detail; the
harmonic and soliton limits of the Whitham modulation equations are explicitly
written down; and the reduction of the Whitham equations to those for the
radial NLS equation is explicitly carried out. Finally, the extension of the
theory to higher spatial dimensions is briefly outlined. The multidimensional
NLS-Whitham equations obtained here may be used to study large amplitude
wavetrains in a variety of applications including nonlinear photonics and
matter waves.

##  1 Introduction

The nonlinear Schrödinger (NLS) equation in one, two and three spatial
dimensions is a ubiquitous model in nonlinear science. One reason is its
universality as a model for the evolution of weakly nonlinear dispersive wave
trains [9, 18, 53]. The NLS equation arises as the governing equation in a
broad variety of physical contexts, ranging from water waves to optics,
acoustics, Bose-Einstein condensates and beyond [6, 36, 39, 44]. As a result,
enormous attention has been devoted over the last half century to the study of
its solutions. It is also the case that in many physical situations,
dispersive effects are much weaker than nonlinear ones and these scenarios,
which can be formulated as small dispersion limits of the governing equations,
give rise to a variety of interesting physical phenomena [25]. In particular,
the small dispersion limits often lead to the formation of dispersive shock
waves, a coherent, slowly modulated and expanding train of nonlinear
oscillations.

A powerful tool in the study of small dispersion limits is Whitham modulation
theory (also simply called Whitham theory) [56, 57]. Whitham theory is an
asymptotic framework within which one can derive the Whitham modulation
equations or Whitham equations for brevity. The Whitham equations are a system
of first-order, quasi-linear partial differential equations (PDEs) that govern
the evolution of the periodic traveling wave solutions of the original PDE
over spatial and temporal scales that are larger than the traveling wave
solution’s wavelength and period, respectively. Whitham theory does not
require integrability of the original PDE, and therefore it can also be
applied to non-integrable PDEs. Thanks to Whitham theory and, when applicable,
the inverse scattering transform (IST), much is known about small dispersion
limits for (1+1)-dimensional nonlinear wave equations (e.g., see [13, 21, 25,
30, 38, 45] and references therein). On the other hand, small dispersion
limits for (2+1)-dimensional systems have been much less studied and
(3+1)-dimensional systems apparently have not been studied at all. Recently,
the Whitham modulation equations for the Kadomtsev-Petviashvili (KP) and two-
dimensional Benjamin-Ono equations and, more generally, a class of
(2+1)-dimensional equations of KP type were derived [2, 3, 1]. The properties
of the resulting KP-Whitham equations were then studied in [14] and the
soliton limit of these equations was used in [48, 50, 49] to study the time
evolution of a variety of piecewise-constant initial conditions in the
modulation equations and, in the process, characterize the resulting dynamics
of the solutions of the KP equation. Recently, the Whitham equations for the
radial NLS equation [4] and those for focusing and defocusing two-dimensional
nonlinear Schrödinger (2DNLS) equations [5] were also derived using a multiple
scales approach.

The goal of this work is to derive and study the Whitham modulation equations
for the defocusing multi-dimensional nonlinear Schrödinger equation, which we
write in the semiclassical scaling as

$i\varepsilon\psi_{t}+\varepsilon^{2}\nabla^{2}\psi-2|\psi|^{2}\psi=0\,$ (1.1)

for a complex-valued field $\psi(\mathbf{x},t)$, where
$\mathbf{x}=(x_{1},\dots,x_{N})^{T}$ and
$\nabla^{2}\psi=\psi_{x_{1}x_{1}}+\cdots+\psi_{x_{N}x_{N}}$ is the spatial
Laplacian, and subscripts $x_{j}$ and $t$ denote partial differentiation
throughout. Equation (1.1) arises as a governing equation in water waves [6],
optics [44], plasmas [36], Bose-Einstein condensates [39], magnetic materials
[59] and beyond. The small parameter $0<\varepsilon\ll 1$ quantifies the
relative strength of dispersive effects compared to nonlinear ones and sets a
spatial and temporal scale for oscillatory solutions. In the (1+1)-dimensional
case, the Whitham modulation equations have been shown to provide quantitative
predictions for experiments in ultracold quantum fluids [34, 35] and nonlinear
optics [55, 58, 10, 8].

While the Whitham equations for the two-dimensional version of (1.1)
(hereafter referred to as the 2DNLS equation) were obtained in [5], this work
differs from [5] in several important respects. First, our derivation employs
a two-phase ansatz for the periodic solutions of the 2DNLS equation, which has
several practical advantages. For one thing, it immediately yields a second
conservation of waves equation in vector form that was missed in [5]. It is
well known that several methods can be used to derive the Whitham equations:
averaged conservation laws, averaged Lagrangian, and multiple scales
perturbation theory. Our derivation employs averaged conservation laws which
are directly tied to the physical symmetries of the NLS equation, rather than
secularity conditions as used in [5]. Moreover, the ability to take advantage
of the second conservation of waves equation also simplifies the calculations.
In contrast, one of the secularity conditions obtained in [5] is equivalent to
the averaged energy equation, which is more complicated and requires more
significant manipulation than the second conservation of waves equation. Our
approach dramatically simplifies the calculations and enables us to carry out
the whole derivation in vector form. Consequently, the resulting NLS-Whitham
equations are obtained in a simpler way, which lays the groundwork for
generalizations to other NLS-type equations and higher dimensions.

In this work, we also show how our approach allows one to easily generalize
the derivation of the Whitham equations to the NLS equation in an arbitrary
number of spatial dimensions. We primarily concentrate on the two and three
dimensional cases, though some of our results apply to an arbitrary number of
spatial dimensions. This generalization to higher dimensions is particularly
relevant because the NLS equation in three spatial dimensions is the zero-
potential version of the Gross-Pitaevski equation, and is therefore of
fundamental importance in describing the dynamics of Bose-Einstein condensates
[39], so we expect our results to be directly applicable in that context.

We use our representation of the NLS-Whitham equations to identify several
symmetries and reductions of the Whitham equations. For example, we verify
that the Whitham equations preserve the invariance of the (N+1)-dimensional
NLS equation with respect to scaling and Galilean transformations, and we take
advantage of the vector formulation of the modulation equations, which we use
to show that they preserve the rotation symmetry of the multidimensional NLS
equation. We also explicitly write down both the harmonic and soliton limits
of the Whitham equations in a mathematically convenient set of independent
variables (which we refer to as Riemann-type variables) and in physical
variables. We identify the self-consistent reduction of the 2DNLS-Whitham
equations to the Whitham equations for the radial NLS equation.

The outline of this work is as follows. In section 2 we write the NLS equation
in hydrodynamic form, write down its conservation laws, and obtain a
representation for the periodic solutions. In section 3 we average the
conservation laws to obtain the Whitham equations in physical variables. In
section 4 we begin to study the reductions of the Whitham equations in
physical variables, including one-dimensional reductions as well as the
harmonic and soliton limits. In section 5 we discuss two different
transformations to Riemann-type variables. In section 6 we derive further
symmetries and reductions of the Whitham equations, including the reduction to
the Whitham equations of the radial NLS equation and the harmonic and soliton
limits of the Whitham equations in Riemann-type variables. In section 7 we
present the generalization of the results to the NLS equation in three spatial
dimensions, and in section 8 we end this work with a discussion of the results
and some final remarks. The details of various calculations are relegated to
the Appendix.

##  2 Hydrodynamic form, conservation laws and periodic solutions of the NLS
equation

###  2.1 Madelung form of the NLS equation and its conservation laws

We begin by writing down the first few conservation laws of the NLS equation
(1.1) in an arbitrary number of dimensions. It is convenient to introduce the
Madelung transformation

$\displaystyle\psi(\mathbf{x},t)=\sqrt{\rho(\mathbf{x},t)}\,e^{i\Phi(\mathbf{x},t)}\,,$
(2.0a)
$\displaystyle\mathbf{u}(\mathbf{x},t)=\varepsilon{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\Phi(\mathbf{x},t)\,.$ (2.0b)

where $\mathbf{u}=(u_{1},\dots,u_{N})^{T}$,
$\mathbf{x}=(x_{1},\dots,x_{N})^{T}$ and ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}=(\partial_{x_{1}},\dots,\partial_{x_{N}})^{T}$.
Substituting (2.1) into the NLS equation (1.1), separating into real and
imaginary parts, and differentiating the real part with respect to each of the
spatial variables yields the following dispersive hydrodynamic system of PDEs:

$\displaystyle\rho_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\rho\mathbf{u})=0\,,$
(2.0a) $\displaystyle\mathbf{u}_{t}+2(\mathbf{u}\cdot{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}})\,\mathbf{u}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\rho-{\textstyle\frac{1}{4}}\varepsilon^{2}{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\Big{(}\nabla^{2}\ln\rho+\frac{1}{\rho}\nabla^{2}\rho\Big{)}=0\,.$
(2.0b)

The conservation laws for (1.1) for the mass $E$, momentum $\mathbf{P}$ and
energy $H$ in integrated form are:

${\frac{dE}{dt}}=0\,,\qquad{\frac{d\mathbf{P}}{dt}}=0\,,\qquad{\frac{dH}{dt}}=0\,,$
(2.0a)

where

$\displaystyle
E=\int_{\mathbb{R}^{N}}|\psi|^{2}\,({\mathrm{d}}\mathbf{x})\,,\quad\mathbf{P}={\textstyle\frac{i}{2}}\varepsilon\int_{\mathbb{R}^{N}}(\psi{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi^{*}-\psi^{*}{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi)\,({\mathrm{d}}\mathbf{x})\,,\quad
H=\int_{\mathbb{R}^{N}}\big{(}\varepsilon^{2}\|{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi\|^{2}+|\psi|^{4}\big{)}\,({\mathrm{d}}\mathbf{x})\,,$

$\|\mathbf{v}\|^{2}=|v_{1}|^{2}+\cdots+|v_{N}|^{2}$ is the Euclidean vector
norm and $({\mathrm{d}}\mathbf{x})={\mathrm{d}}x_{1}\cdots{\mathrm{d}}x_{N}$
is the volume element in $\mathbb{R}^{N}$. These conservation laws correspond,
via Noether’s theorem, to the invariance of the NLS equation (1.1) with
respect to phase rotations, space and time translations, respectively [53]. In
differential form, and in terms of the Madelung variables, these conservation
laws become [37]

$\displaystyle\rho_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\rho\mathbf{u})=0\,,$
(2.0a) $\displaystyle(\rho\mathbf{u})_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\rho\,\mathbf{u}\otimes\mathbf{u})+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\rho^{2})={\textstyle\frac{1}{2}}\varepsilon^{2}\Big{(}{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\nabla^{2}\rho)-{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\Big{(}\frac{1}{\rho}\nabla\rho\otimes\nabla\rho\Big{)}\Big{)}\,,$
(2.0b) $\displaystyle h_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\big{(}(h+\rho^{2})\mathbf{u}\big{)}=\varepsilon^{2}{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\Big{(}\mathbf{u}\nabla^{2}\rho-\frac{1}{\rho}({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\rho\mathbf{u})\nabla\rho\Big{)}\,,$
(2.0c)

where $\otimes$ denotes the dyadic [namely,
$\mathbf{v}\otimes\mathbf{w}=\mathbf{v}\mathbf{w}^{T}$, so that
$(\mathbf{v}\otimes\mathbf{w})_{i,j}=v_{i}w_{j}$] and the mass density,
momentum density and energy density of (1.1) are, respectively

$\displaystyle\rho=|\psi|^{2},\quad\rho\mathbf{u}=i2\varepsilon(\psi{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi^{*}-\psi^{*}{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi),\qquad
h=\varepsilon^{2}\|{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\psi\|^{2}+|\psi|^{4}=\rho\|\mathbf{u}\|^{2}+\rho^{2}+\frac{\varepsilon^{2}}{4\rho}\|{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\rho\|^{2}\,.$

The first two of the conservation laws (2.1) are equivalent to the real and
imaginary parts of NLS equation in hydrodynamic form (2.1), but only up to an
extra differentiation, an issue that we will return to later.

###  2.2 Periodic solutions of the NLS equation via a two-phase ansatz

The Whitham modulation equations govern the slow dynamics of the parameters of
the periodic solutions of the PDE of interest. Next, we therefore write down
the periodic solutions of the hydrodynamic system (2.1) in arbitrary
dimensions. We begin by looking for solutions in the form of the following
two-phase ansatz:

$\rho(\mathbf{x},t)=\rho(Z)\,,\qquad\Phi(\mathbf{x},t)=\phi(Z)+S\,,$ (2.0a)

where $\rho(Z)$ and $\phi(Z)$ are periodic function of $Z$ with period one,
and the “fast phases” $Z$ and $S$ are

$Z(\mathbf{x},t)=(\mathbf{k}\cdot\mathbf{x}-\omega t)/\varepsilon\,,\qquad
S(\mathbf{x},t)=(\mathbf{v}\cdot\mathbf{x}-\mu t)/\varepsilon\,.$ (2.0b)

where $\mathbf{k}=(k_{1},\dots,k_{N})^{T}$ and
$\mathbf{v}=(v_{1},\dots,v_{N})^{T}$. The reason for using a two-phase ansatz
is the fact that the solution $\psi(\mathbf{x},t)$ of the NLS equation (1.1)
is complex-valued, unlike that of the Korteweg-deVries (KdV) equation (of
which the KP equation mentioned earlier is a two-dimensional generalization),
which is real-valued. Therefore, a one-phase ansatz (e.g., as in [5]) leads
only to a subclass of all periodic solutions, and one would need to apply a
Galilean boost a posteriori in order to capture the most general family of
periodic solutions of the NLS equation. Two-phase ansatzes are standard when
deriving the Whitham equations using Lagrangian averaging (e.g., see [57]);
the novelty here is that such a two-phase ansatz is combined with the use of
averaged conservation laws. A key benefit of this approach is the immediate
deduction of an additional conservation law compared to [5].

In light of (2.2), the definition (2.1) yields

$\mathbf{u}(Z)=\mathbf{k}\,\phi^{\prime}(Z)+\mathbf{v}\,,$ (2.1)

using primes to denote derivatives with respect to $Z$ for brevity. The fact
that $\phi(Z)$ is periodic implies

$\overline{{\mathbf{u}}}=\mathbf{v}\,,$ (2.2)

where throughout this work the overbar will denote the integral of a quantity
with respect to $Z$ over the unit period. Moreover, the definition (2.1)
implies the irrotationality condition

${\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\mathbf{u}=0\,.$ (2.3)

Hereafter, $\mathbf{v}\wedge\mathbf{w}$ is the $N$-dimensional wedge product,
which in two and three spatial dimensions can be replaced by the standard
cross product [17, 29]. We substitute the two-phase ansatz (2.2) into the
hydrodynamic equations (2.0a) and (2.0b) and collect the leading-order terms,
obtaining:

$\displaystyle-\omega\rho^{\prime}+2\mathbf{k}\cdot(\rho\mathbf{u})^{\prime}=0\,,$
(2.0a)
$\displaystyle-\omega\mathbf{u}^{\prime}+2(\mathbf{k}\cdot\mathbf{u})\,\mathbf{u}^{\prime}+2\mathbf{k}\rho^{\prime}-{\textstyle\frac{1}{4}}\Big{(}(\ln\rho)^{\prime\prime}+\frac{\rho^{\prime\prime}}{\rho}\Big{)}^{\prime}{\|\mathbf{k}\|^{2}}\mathbf{k}=\mathbf{0}\,.$
(2.0b)

Integrating (2.0a) and using (2.1) yields

$\phi^{\prime}(Z)=\frac{1}{\|\mathbf{k}\|}\bigg{(}U+\frac{J}{\rho}-\hat{\mathbf{k}}\cdot\bar{\mathbf{u}}\bigg{)}\,,$
(2.1)

where $U=\omega/(2\|\mathbf{k}\|)$ is the phase speed,
$\hat{\mathbf{k}}=\mathbf{k}/\|\mathbf{k}\|$, and the integration constant $J$
will be determined later. Using (2.1), we can rewrite (2.1) as:

$\mathbf{u}(Z)=\Big{(}\frac{J}{\rho}+U\Big{)}\,\hat{\mathbf{k}}+\bar{\mathbf{u}}_{\perp}\,,$
(2.2)

where
${\mathbf{u}}_{\perp}=\mathbf{u}-(\hat{\mathbf{k}}\cdot\mathbf{u})\,\hat{\mathbf{k}}$.
Importantly, the requirement that $\phi(Z)$ is periodic implies that
$\phi^{\prime}(Z)$ must have zero mean. Taking the inner product of (2.2) with
$\hat{\mathbf{k}}$ and averaging the result over the wave period yields a
relation between $\bar{\mathbf{u}}$ and $U$, and therefore determines
$\omega=2\|\mathbf{k}\|U$:

$U=\hat{\mathbf{k}}\cdot\bar{\mathbf{u}}-J\,\overline{\rho^{-1}}\,.$ (2.3)

Next, substituting (2.1) into (2.0b) and simplifying yields two ODEs for
$\rho$. Note that the two ODEs are consistent thanks to the constraint (2.3),
which becomes, to leading order,

$\mathbf{k}\wedge\mathbf{u}^{\prime}=0\,.$ (2.4)

Integrating the resulting ODE for $\rho$ one obtains [see Appendix A.2 for
details]

$\displaystyle(\rho^{\prime})^{2}=P_{3}(\rho)\,,$ (2.5) $\displaystyle
P_{3}(\rho)=\frac{4}{\|\mathbf{k}\|^{2}}(\rho-\lambda_{1})(\rho-\lambda_{2})(\rho-\lambda_{3})\,,$
(2.6)

whose solution is

$\rho(Z)=A+4m\|\mathbf{k}\|^{2}K_{m}^{2}\mathop{\rm
sn}\nolimits^{2}(2K_{m}z|m)\,,$ (2.7)

where $A$ is a free parameter, and with

$J^{2}=A\,\big{(}A+4\|\mathbf{k}\|^{2}K_{m}^{2}\big{)}\big{(}A+4m\|\mathbf{k}\|^{2}K_{m}^{2}\big{)}\,,$
(2.8)

The roots $\lambda_{1},\dots,\lambda_{3}$ are related to the coefficients in
the solution (2.7) as

$\displaystyle\lambda_{1}=A,\qquad\lambda_{2}=A+4mK_{m}^{2}\|\mathbf{k}\|^{2},\qquad\lambda_{3}=A+4K_{m}^{2}\|\mathbf{k}\|^{2}.$
(2.9)

Conversely, when $\lambda_{1},\lambda_{2},\lambda_{3}$ are known, $A$,
$\|\mathbf{k}\|$ and $m$ are given by

$A=\lambda_{1}\,\qquad\|\mathbf{k}\|^{2}=(\lambda_{3}-\lambda_{1})/4K_{m}^{2}\,,\qquad
m=(\lambda_{2}-\lambda_{1})/(\lambda_{3}-\lambda_{1})\,.$ (2.10)

The amplitude of the periodic oscillations of the density is
$\lambda_{2}-\lambda_{1}$. The requirements $\rho\geqslant 0$,
$\|\mathbf{k}\|\geqslant 0$ and $0\leqslant m\leqslant 1$ immediately yield
the constraints $A\geqslant 0$ as well as

$0\leqslant\lambda_{1}\leqslant\lambda_{2}\leqslant\lambda_{3}\,.$ (2.11)

The symmetric polynomials $e_{1},\dots,e_{3}$ defined by the roots
$\lambda_{1},\dots,\lambda_{3}$ will also be useful later:

$\displaystyle e_{1}=\lambda_{1}+\lambda_{2}+\lambda_{3}\,,\qquad
e_{2}=\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{3}\lambda_{1}\,,\qquad
e_{3}=\lambda_{1}\lambda_{2}\lambda_{3}=J^{2}\,.$ (2.12)

Note that (2.12) only determines $J$ up to a sign, i.e.,
$J=\sigma\sqrt{\lambda_{1}\lambda_{2}\lambda_{3}}$, with $\sigma=\pm 1$. Both
sign choices lead to valid solutions of the NLS equation (1.1). Some care is
deserved when determining the value of $\sigma$ in the presence of modulations
of the periodic solutions, as discussed in section 3.2.

The leading-order periodic solution of the hydrodynamic system (2.1) defined
by (2.1) and (2.7) contains the following independent parameters: $A$, $m$,
$\mathbf{k}$, $\bar{\mathbf{u}}$ and $\mu$. However, recall that, to derive
the hydrodynamic equation (2.0b) from the NLS equation (1.1), one
differentiates the real part with respect to the spatial variables. Imposing
that the solution of the dispersive hydrodynamic system (2.1) also solves the
NLS equation (by substituting into the undifferentiated imaginary part of the
NLS equation (1.1)) yields a constraint that determines $\mu$ in terms of the
other constants. Deriving this relation directly from the above expressions is
a bit cumbersome, but seeking a periodic solution of (1.1) without writing it
in hydrodynamic form [cf. Appendix A.1], one obtains

$\mu=4(1+m)\|\mathbf{k}\|^{2}K_{m}^{2}+3A+\|\mathbf{\bar{u}}\|-\big{(}J\overline{\rho^{-1}}\big{)}^{2}\,.$
(2.13)

One can now verify that adding this relation to the above solution of the
hydrodynamic system does indeed yield a solution of the NLS equation (1.1).
Alternatively, one can obtain (2.13) using the undifferentiated version of
(2.0b); see Appendix A.2. Thus, the periodic solutions of the NLS equation
(1.1) in $N$ spatial dimensions contain $2N+2$ scalar independent parameters:
$A$, $m$, $\mathbf{k}$ and $\mathbf{v}=\bar{\mathbf{u}}$, as one would expect
based on the invariances of the PDE (cf. [53]).

###  2.3 Harmonic and soliton limits of the periodic solutions

Recall that the harmonic ($m=0$) and soliton ($m=1$) limits of the Whitham
equations for the one-dimensional NLS (1DNLS) equation have special
significance [25]. The same will be true for the multi-dimensional NLS
equation. It is therefore useful to evaluate the corresponding limits of the
above periodic solutions.

In the limit $m\to 0$ (i.e., $\lambda_{2}\to\lambda_{1}^{+}$), the solution
(2.1) reduces to a plane wave. Indeed, in this limit, we have

$\displaystyle\rho(Z)=A,\,\qquad
B=0\,,\qquad\mu=2A+\|\bar{\mathbf{u}}\|^{2}\,,\qquad
J^{2}=A^{2}(\pi^{2}\|\mathbf{k}\|^{2}+A)\,,$ (2.14)

and

$\psi(\mathbf{x},t)=\sqrt{A}\,\,{\mathrm{e}}^{i(\bar{\mathbf{u}}\cdot\mathbf{x}-(\|\bar{\mathbf{u}}\|^{2}+2A)t)}\,.$
(2.15)

Therefore, the only independent parameters in this case are $A$ and
$\bar{\mathbf{u}}$.

In the opposite limit ($m\to 1$, i.e., $\lambda_{2}\to\lambda_{3}^{-}$), the
solution (2.1) reduces to the soliton solution of the NLS equation. Indeed, in
this limit, (2.7) and (2.10) yield

$\displaystyle\displaystyle\rho(Z)=\lambda_{1}+(\lambda_{3}-\lambda_{1})\tanh^{2}\big{[}\sqrt{\lambda_{3}-\lambda_{1}}\,\big{(}\hat{\mathbf{k}}\cdot\mathbf{x}-\omega
t/\|\mathbf{k}\|\big{)}\big{]},\,$ (2.0a) $\displaystyle
B=\lambda_{3}-\lambda_{1}\,,\quad J^{2}=\lambda_{1}\lambda_{3}^{2}\,,\quad
U=\hat{\mathbf{k}}\cdot\bar{\mathbf{u}}-\sigma\sqrt{\lambda_{1}}\,,\quad\mu=2\lambda_{3}+\|\bar{\mathbf{u}}\|^{2}\,.$
(2.0b)

Note that $\|\mathbf{k}\|\to 0$ as $m\to 1$, but $K_{m}\to\infty$ in such a
way that their product remains finite:
$\|\mathbf{k}\|K_{m}\to\sqrt{\lambda_{3}-\lambda_{1}}/2$. Using (2.1) we then
obtain

$\phi(Z)=\arctan\Big{[}\sqrt{\lambda_{3}-\lambda_{1}}\,\tanh\Big{(}\sqrt{\lambda_{3}-\lambda_{1}}\,\big{(}\hat{\mathbf{k}}\cdot\mathbf{x}-\omega
t/\|\mathbf{k}\|\big{)}\Big{)}/\sqrt{\lambda_{1}}\big{]}\,,$ (2.1)

implying

$e^{i\phi+iS}=e^{iS}\Big{[}\sqrt{\lambda_{1}}+i\sqrt{\lambda_{3}-\lambda_{1}}\tanh\Big{(}\sqrt{\lambda_{3}-\lambda_{1}}\,\big{(}\hat{\mathbf{k}}\cdot\mathbf{x}-\omega
t/\|\mathbf{k}\|\big{)}\Big{)}\Big{]}\Big{/}\sqrt{\rho(Z)}\,,$ (2.2)

with $S=\bar{\mathbf{u}}\cdot\mathbf{x}-\mu t$ as before. Putting everything
together, we obtain

$\displaystyle\psi(\mathbf{x},t)=A_{o}{\mathrm{e}}^{-2iA_{o}^{2}t}{\mathrm{e}}^{i(\bar{\mathbf{u}}\cdot\mathbf{x}-\|\bar{\mathbf{u}}\|^{2}t)}\big{\\{}\cos\theta+i\sin\theta\tanh[A_{o}\sin\theta\,[\hat{\mathbf{k}}\cdot\mathbf{x}-2(\hat{\mathbf{k}}\cdot\bar{\mathbf{u}}-A_{o}\cos\theta)t]]\big{\\}}\,,$

with $\bar{\mathbf{u}}$ as in (2.2), $A_{o}=\sqrt{\lambda_{3}}$ and
$\theta=\arctan\big{[}\sqrt{(\lambda_{3}-\lambda_{1})/\lambda_{1}}\big{]}$.
The independent parameters of the solution in this case are $\lambda_{1}$,
$\lambda_{3}$ (or equivalently $A_{o}$ or $\theta$), $\hat{\mathbf{k}}$ and
$\bar{\mathbf{u}}$. One can further reduce (2.3) to the more familiar form of
the dark soliton solutions of the defocusing NLS equation by choosing
$\bar{\mathbf{u}}=\mathbf{0}$.

##  3 Derivation of the NLS-Whitham equation by averaged conservation laws

We are now ready to study slow modulations of the periodic solutions described
above and derive the Whitham modulation equations that govern them.

###  3.1 Nonlinear modulations and averaged conservation laws

We begin by introducing the following multiple scales ansatz for the solution
of the NLS equation (1.1):

$\rho(\mathbf{x},t)=\rho(Z,\mathbf{X},T)\,,\qquad\Phi(\mathbf{x},t)=\phi(Z,\mathbf{X},T)+S\,,$
(3.1)

where $\mathbf{X}=\mathbf{x}$ and $T=t$, with $\rho$ and $\phi$ periodic in
$Z$ with period one and

$\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}Z=\frac{\mathbf{k}(\mathbf{X},T)}{\varepsilon}\,,\qquad
Z_{t}=-\frac{\omega(\mathbf{X},T)}{\varepsilon}\,,$ (3.0a)
$\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}S=\frac{\mathbf{v}(\mathbf{X},T)}{\varepsilon}\,,\qquad
S_{t}=-\frac{\mu(\mathbf{X},T)}{\varepsilon}\,,$ (3.0b)

where, as per the results of section 2.2, $\mathbf{v}=\bar{\mathbf{u}}$. The
above multiple scales ansatz implies

$\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\mapsto\frac{\mathbf{k}}{\varepsilon}\partial_{Z}+\frac{\mathbf{v}}{\varepsilon}\partial_{S}+\nabla_{\mathbf{X}}\,,\qquad\partial_{t}\mapsto-\frac{\omega}{\varepsilon}\partial_{Z}-\frac{\mu}{\varepsilon}\partial_{S}+\partial_{T}\,,$
(3.1)

Substituting (3.1) into (1.1), to leading order we recover the periodic
solution (2.1), but where all $2N+2$ parameters $A$, $m$, $\mathbf{k}$ and
$\bar{\mathbf{u}}$ are now slowly varying functions of $\mathbf{X}$ and $t$.
We then seek modulation equations to determine the space-time dependence of
these parameters. To avoid complicating the notation unncessarily, below we
will write derivatives in $\mathbf{X}$ and $T$ as derivatives in $\mathbf{x}$
and $t$. Equations (3.1) immediately yield the equations of conservation of
waves:

$\displaystyle\mathbf{k}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\omega=\mathbf{0}\,,$ (3.0a)
$\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\mathbf{k}=0\,,$ (3.0b)
$\displaystyle\bar{\mathbf{u}}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\mu=\mathbf{0}\,,$ (3.0c)
$\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\bar{\mathbf{u}}=0\,.$
(3.0d)

Of course only $N$ equations among (3.0a) and (3.0b) are independent, and
similarly for (3.0c) and (3.0d). Equations (3.0a) and (3.0c) form the first
two vectorial Whitham modulation equations, whereas (3.0b) and (3.0d) are
compatibility constraints.

Next, we obtain the remaining Whitham modulation equations by averaging the
conservation laws (2.1) over the fast variable $Z$. Using (3.1) to replace all
spatial and temporal derivatives in (2.1) and (2.1), expanding all terms in
powers of $\varepsilon$, and averaging, we obtain at order
$\mathcal{O}(\varepsilon^{0})$

$\displaystyle(\bar{\rho})_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\overline{\rho\mathbf{u}})=0\,,$
(3.0e) $\displaystyle(\overline{\rho\mathbf{u}})_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\overline{\rho\mathbf{u}\otimes\mathbf{u}})+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\overline{\rho^{2}})+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\bigg{(}\overline{\frac{(\rho^{\prime})^{2}}{2\rho}}\,\mathbf{k}\otimes\mathbf{k}\bigg{)}=\mathbf{0}\,,$
(3.0f) $\displaystyle\bar{h}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\Big{(}2\overline{h\mathbf{u}}+2\overline{\rho^{2}\mathbf{u}}+\Big{(}\mathbf{k}\cdot\\!\overline{\frac{\rho^{\prime}}{\rho}(\rho\mathbf{u})^{\prime}}\Big{)}\,\,\mathbf{k}-\|\mathbf{k}\|^{2}\overline{\rho^{\prime\prime}\mathbf{u}}\Big{)}=0\,,$
(3.0g)

where $\bar{h}$ denotes the averaged energy density:

$\bar{h}=\overline{\rho\|\mathbf{u}\|^{2}}+\overline{\rho^{2}}+\frac{1}{4}{\|\mathbf{k}\|^{2}}\overline{{(\rho^{\prime})^{2}}/{\rho}}\,.$
(3.1)

Together with (3.0a) and (3.0c), equations (3.0e)–(3.0g) are $3N+2$ scalar
PDEs for the $2N+2$ dependent variables $A$, $m$, $\mathbf{k}$ and
$\mathbf{v}=\bar{\mathbf{u}}$ subject to the $2N$ spatial constraints (3.0b),
(3.0d), and are the desired Whitham modulation equations in physical variables
in any number of spatial dimensions. Of course, not all of these equations are
independent. We will see later that choosing different subsets of equations
still leads to equivalent results, and in the end the number of independent
modulation equations is $2N+2$. At the same time, however, we emphasize the
simplicity and directness of this approach compared to [5] in deriving the
Whitham equations in multiple spatial dimensions.

###  3.2 Modified form of the modulation equations

In preparation for further simplification of the above system of Whitham
equations, it is convenient to express the periodic solutions in terms of the
roots $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, thereby replacing $A$, $m$
and $\|\mathbf{k}\|^{2}$ as dependent variables. Explicitly, (2.2) and (2.7)
become:

$\displaystyle\rho(Z)=\lambda_{1}+(\lambda_{2}-\lambda_{1})\,\mathop{\rm
sn}\nolimits^{2}(2K_{m}z|m)\,,$ (3.0a)
$\displaystyle\mathbf{u}(Z)=\mathbf{U}+\frac{J}{\rho(Z)}\,\hat{\mathbf{k}}\,,$
(3.0b)

with

$\displaystyle\mathbf{U}=\bar{\mathbf{u}}-J\overline{\rho^{-1}}\,\hat{\mathbf{k}}\,,\qquad$
(3.0c)
$\displaystyle\omega=2\mathbf{k}\cdot\mathbf{U}\,,\qquad\mu=\lambda_{1}+\lambda_{2}+\lambda_{3}+\|\mathbf{U}\|^{2}+2UJ\overline{\rho^{-1}}\,,$
(3.0d)

with $\hat{\mathbf{k}}=\mathbf{k}/\|\mathbf{k}\|$ as before and $J$, $A$,
$\|\mathbf{k}\|$ and $m$ given in terms of $\lambda_{1},\dots,\lambda_{3}$ by
(2.8) and (2.10). In turn, using (3.2), we can write the Whitham modulation
equations (3.1) as

$\displaystyle\mathbf{k}_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\mathbf{k}\cdot\mathbf{U})=\mathbf{0}\,,$
(3.0a) $\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\mathbf{k}=0\,,$ (3.0b)
$\displaystyle\big{(}\mathbf{U}+J\overline{\rho^{-1}}\,\hat{\mathbf{k}}\big{)}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\big{(}e_{1}+\|\mathbf{U}\|^{2}+2J\overline{\rho^{-1}}\,\mathbf{U}\cdot\hat{\mathbf{k}}\big{)}=\mathbf{0}\,,$
(3.0c) $\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\big{(}\mathbf{U}+J\overline{\rho^{-1}}\,\hat{\mathbf{k}}\big{)}=0\,,$
(3.0d)
$\displaystyle\bar{\rho}_{t}+2\nabla\cdot\big{(}J\hat{\mathbf{k}}+\bar{\rho}\mathbf{U}\big{)}=0\,,$
(3.0e)
$\displaystyle(J\,\hat{\mathbf{k}}+\overline{\rho}\mathbf{U})_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\overline{\rho^{2}})+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\Big{[}\Big{(}2\bar{\rho}\mathbf{U}+2J\,\hat{\mathbf{k}}\Big{)}\otimes\mathbf{U}+{2J}\mathbf{U}\otimes\hat{\mathbf{k}}+\frac{2}{3}\big{(}2e_{2}-e_{1}\bar{\rho}\big{)}\,\hat{\mathbf{k}}\otimes\hat{\mathbf{k}}\Big{]}=\mathbf{0}\,.$
(3.0f) $\displaystyle\bar{h}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\Bigg{[}2J(2\,\bar{\rho}+\overline{\|\mathbf{u}\|^{2}})\,\hat{\mathbf{k}}+2(\overline{\rho^{2}}+\bar{h})\,\mathbf{U}+\Bigg{(}\mathbf{U}\cdot\hat{\mathbf{k}}\,\Bigg{(}\overline{\frac{(\rho^{\prime})^{2}}{\rho}}\Bigg{)}-\frac{J}{2}\Bigg{(}\overline{\frac{(\rho^{\prime})^{2}}{\rho^{2}}}\Bigg{)}\Bigg{)}\,\mathbf{k}\Bigg{]}=0\,.$
(3.0g)

See Appendix A.2 for details on how to obtain (3.0f). The next step is the
evaluation of the elliptic integrals in (3.2). To this end, we have [46]

$\displaystyle\overline{\rho}=\int_{0}^{1}\rho(Z)\,{\mathrm{d}}z=\lambda_{3}-(\lambda_{3}-\lambda_{1})\frac{E_{m}}{K_{m}}\,,$
(3.0a)
$\displaystyle\overline{\rho^{-1}}=\int_{0}^{1}\rho^{-1}(Z)\,{\mathrm{d}}z=\frac{1}{\lambda_{1}K_{m}}\,\Pi\Big{(}1-\frac{\lambda_{2}}{\lambda_{1}}\Big{|}m\Big{)}\,,$
(3.0b)

where $K_{m}=K(m)$, $E_{m}=E(m)$ and $\Pi(\cdot|m)$ are the complete elliptic
integrals of the first, second and third kind respectively. We also note, for
convenience, that

$\displaystyle\bar{\mathbf{u}}=\int_{0}^{1}\mathbf{u}(Z)\,{\mathrm{d}}z=\mathbf{U}+\sigma\frac{\sqrt{\lambda_{2}\lambda_{3}}}{\sqrt{\lambda_{1}}K_{m}}\Pi\Big{(}1-\frac{\lambda_{2}}{\lambda_{1}}\Big{|}m\Big{)}\hat{\mathbf{k}}\,,$
(3.0a)
$\displaystyle\overline{\rho\mathbf{u}}=\int_{0}^{1}\rho(Z)\mathbf{u}(Z)\,{\mathrm{d}}z=\bar{\rho}\mathbf{U}+J\hat{\mathbf{k}}=\Big{(}\lambda_{3}-(\lambda_{3}-\lambda_{1})\frac{E_{m}}{K_{m}}\Big{)}\mathbf{U}+\sigma\sqrt{\lambda_{1}\lambda_{2}\lambda_{3}}\hat{\mathbf{k}}\,.$
(3.0b)

We reiterate that not all of the equations (3.2) are independent. For example,
one can obtain (3.0d) using (3.0c) and (3.0e). This is relevant because it
allows us to work with the most convenient subset of equations among all the
PDEs in (3.2), as long as the compatibility constraints (3.0b) and (3.0d) are
satisfied. To this end, recall that $\bar{h}$ is given by (3.1), and

$\displaystyle\overline{\|\mathbf{u}\|^{2}}=\|\mathbf{U}\|^{2}+2J\mathbf{U}\cdot\hat{\mathbf{k}}\,\overline{\rho^{-1}}+J^{2}\overline{\rho^{-2}}\,,$
(3.0a)
$\displaystyle\overline{\rho\|\mathbf{u}\|^{2}}=J^{2}\,\overline{\rho^{-1}}+2J\mathbf{U}\cdot\hat{\mathbf{k}}+\bar{\rho}\|\mathbf{U}\|^{2}\,.$
(3.0b)

Moreover, the terms $(\rho^{\prime})^{2}/\rho$ and
$(\rho^{\prime})^{2}/\rho^{2}$, which appear in (3.0g), can be computed using
(2.5). On the other hand, the averaged energy conservation law (3.0g) is the
most complicated among all of the equations (3.2). In section 7 we will show
that, thanks to the use of the two-phase ansatz and the resulting second
conservation of waves equations (3.0c) and (3.0d), one can avoid having to
deal with the averaged energy equation (3.0g), which greatly simplifies the
transformation to Riemann-type variables.

We also point out that the sign of $J$, as determined by the initial
conditions for the system through the value of $\sigma$—see the discussion
after (2.12)—affects $\bar{\mathbf{u}}$ via (3.0a) and
$\overline{\rho\mathbf{u}}$ via (3.0b). Therefore, when considering
modulations of the periodic solutions, the value of $\sigma$ depends on
$\mathbf{x}$ and $t$, and its value must be chosen in such a way to ensure
smoothness of $\overline{\rho\mathbf{u}}$. In particular, a sign change of $J$
occurs when the solution hits a vacuum point, i.e., $\lambda_{1}=0$. At such a
point, $\bar{\mathbf{u}}$ is singular but $\overline{\rho\mathbf{u}}$ is not.
See [33] for additional discussion.

##  4 Symmetries and reductions of the NLS-Whitham system in physical
variables

We now present several reductions of the Whitham modulation system (3.2) in
physical variables in arbitrary number of spatial dimensions. Further
symmetries and reductions in the two-dimensional case will be discusssed in
section 6 after we introduce Riemann-type variables in section 5.

###  4.1 Unidirectional reductions of the modulation equations

We begin by showing that the NLS-Whitham equations (3.2) reduce to the 1DNLS-
Whitham equations (i.e., the Whitham equations for the 1DNLS equation) when
$k_{2}=\cdots=k_{N}=v_{2}=\cdots=v_{N}=0$ and all quantities are independent
of $x_{2},\dots,x_{N}$. In this case, we have:

$\displaystyle\|\mathbf{k}\|^{2}=k_{1}^{2}\,,\qquad
u_{1}(Z)=\frac{J}{\rho}+U\,,\qquad\omega=2k\,U\,,\qquad
U=\bar{u}_{1}-J\overline{\rho^{-1}}\,,\qquad u_{2}(Z)=0\,.$ (4.0a)

The Whitham equations (3.0b) and (3.0d) and the second components of (3.0a),
(3.0c), and (3.0e) are satisfied trivially, while the rest simplify to:

$\displaystyle k_{t}+2(kU)_{x}=0\,,$ (4.0a)
$\displaystyle\big{(}U+J\overline{\rho^{-1}}\big{)}_{t}+\Big{(}e_{1}+2JU\,\overline{\rho^{-1}}+U^{2}\Big{)}_{x}=0\,,$
(4.0b) $\displaystyle(\bar{\rho})_{t}+2(U\,\bar{\rho}+J)_{x}=0\,,$ (4.0c)
$\displaystyle(U\,\bar{\rho}+J)_{t}+\Big{(}\overline{\rho^{2}}+2U^{2}\bar{\rho}+2J^{2}\overline{\rho^{-1}}+\frac{k^{2}}{2}\overline{\frac{(\rho^{\prime})^{2}}{\rho}}+4U\,J\Big{)}_{x}=0\,,$
(4.0d)
$\displaystyle\Big{(}\overline{\rho^{2}}+U^{2}\bar{\rho}+J^{2}\overline{\rho^{-1}}+\frac{k^{2}}{4}\overline{\frac{(\rho^{\prime})^{2}}{\rho}}+2JU\Big{)}_{t}+\Big{(}\frac{3U\,k^{2}}{2}\overline{\Big{(}\frac{(\rho^{\prime})^{2}}{\rho}\Big{)}}-\frac{k^{2}J}{2}\overline{\frac{(\rho^{\prime})^{2}}{\rho^{2}}}$
$\displaystyle\kern
100.00015pt+4U\overline{\rho^{2}}+(4J+2U^{3})\bar{\rho}+6JU^{2}+6J^{2}U\overline{\rho^{-1}}+2J^{3}\overline{\rho^{-2}}\Big{)}_{x}=0\,,$
(4.0e)

with $x=x_{1}$. The system (4.1) coincides with the modulation equations for
the 1DNLS equation [34] (cf. (4.41) and (4.42) in [34]) upon trivial
rescalings resulting from the different normalization of the NLS equation in
[34] compared to (1.1). Note that (4.1) comprise five PDEs for the four
solution parameters $A,m,k\&U$ (or equivalently
$\lambda_{1},\lambda_{2},\lambda_{3},U$). Once again, one can verify that the
modulation equation obtained from (4.1) is consistent with those obtained from
the first four PDEs above.

The above scenario is not the only one in which the Whitham modulation system
(3.2) reduces to that of the 1DNLS equation. Next we consider so-called
“rotated” one-dimensional reductions where the rotated coordinate frame is
determined by $\mathsf{R}$, an $N\times N$ orthogonal matrix. We introduce the
rotated vector ${\mathbf{w}}^{\sharp}=\mathsf{R}\,\mathbf{w}$ for any vector
$\mathbf{w}$. Then, the rotated one-dimensional reduction is obtained through
the requirement that $\mathbf{k}$ and $\bar{\mathbf{u}}$ (or equivalently
$\hat{\mathbf{k}}$ and $\mathbf{U}$) be parallel and that both depend only on
$t$ and the first component of ${\mathbf{x}}^{\sharp}$. We choose $\mathsf{R}$
so that $\hat{\mathbf{k}}^{\sharp}=(1,0,\dots,0)^{T}$, i.e.,
$k^{\sharp}_{2}=\dots=k^{\sharp}_{N}=0$, which also implies
$U^{\sharp}_{2}=\dots=U^{\sharp}_{N}=0$. Since the Whitham modulation
equations (3.2) are invariant under rotations of the coordinate axes (see
below), we recover the one-dimensional reduction (4.1) when all quantities are
independent of $x^{\sharp}_{2},\dots,x^{\sharp}_{N}$ in the rotated coordinate
frame, i.e., with $x$ and all modulation variables in (4.1) replaced by their
rotations $x^{\sharp}_{1}$, etc.

###  4.2 Invariances of the modulation equations

The Whitham modulation equations (3.2) are manifestly invariant under
translations of the spatial and temporal coordinates. Next we show that the
Whitham system (3.2) preserves the invariance of the NLS equation under
rotations of the Cartesian coordinates. Namely, if
$\mathbf{x}\mapsto\mathbf{x}^{\sharp}=\mathsf{R}\,\mathbf{x}$, where
$\mathsf{R}$ is an arbitrary constant rotation matrix, (3.2) remain unchanged
upon $\mathbf{U}\mapsto\mathbf{U}^{\sharp}=\mathsf{R}\,\mathbf{U}$ and
$\mathbf{k}\mapsto\mathbf{k}^{\sharp}=\mathsf{R}\,\mathbf{k}$. One can verify
that this is indeed the case using the following identities:

$\displaystyle\displaystyle\mathsf{R}\,{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\,,\quad\mathbf{U}\cdot\mathbf{k}=\mathbf{U}^{\sharp}\cdot\mathbf{k}^{\sharp}\,,\quad\|\mathbf{U}\|=\|\mathbf{U}^{\sharp}\|\,,$
(4.0a) $\displaystyle{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\cdot(\alpha\mathbf{k})={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\cdot(\alpha\mathbf{k}^{\sharp})\,,\quad{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\cdot(\alpha\mathbf{U})={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\cdot(\alpha\mathbf{U}^{\sharp})\,,$
(4.0b) $\displaystyle\mathsf{R}\,{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\cdot(\alpha\mathbf{U}\otimes\mathbf{U})={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\cdot(\alpha\mathbf{U}^{\sharp}\otimes\mathbf{U}^{\sharp})\,,\quad\mathsf{R}\,{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\cdot(\alpha\mathbf{k}\otimes\mathbf{k})={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\cdot(\alpha\mathbf{k}^{\sharp}\otimes\mathbf{k}^{\sharp})\,,$
(4.0c) $\displaystyle\mathsf{R}\,{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}}\cdot(\alpha\mathbf{k}\otimes\mathbf{U})={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{x}^{\sharp}}\cdot(\alpha\mathbf{k}^{\sharp}\otimes\mathbf{U}^{\sharp})\,.$
(4.0d)

where $\alpha$ is an arbitrary real number.

Next we show that the Whitham system (3.2) also preserves the invariance of
the NLS equation with respect to scaling and spatial reflections and Galilean
transformations. Recall that, if $q(\mathbf{x},t)$ is any solution of the NLS
equation, so are $q^{\sharp}(\mathbf{x},t)=\alpha
q(\alpha\,\mathbf{x},\alpha^{2}t)$,
$q^{\sharp}(\mathbf{x},t)=q(-\mathbf{x},t)$ and
$q^{\sharp}(\mathbf{x},t)=q(\mathbf{x}-2\mathbf{w}t,t)e^{i(\mathbf{w}\cdot\mathbf{x}-\|\mathbf{w}\|^{2}t)}$
where all transformation parameters are real-valued. We next show that the
modulation equations (3.2) are invariant under each of these transformations.
Specifically, letting
$q^{\sharp}(\mathbf{x},t)=[\rho^{\sharp}(\mathbf{x},t)]^{1/2}\,e^{i\phi^{\sharp}(\mathbf{x},t)}$,
we have, for the scaling symmetry,

$\rho^{\sharp}(\mathbf{x},t)=\alpha^{2}\rho(\alpha\mathbf{x},\alpha^{2}t)\,,\qquad\phi^{\sharp}(\mathbf{x},t)=\phi(\alpha\mathbf{x},\alpha^{2}t)\,,$
(4.1)

and the dependent variables of the Whitham equations become

$\displaystyle\lambda^{\sharp}_{j}(\mathbf{x},t)=\alpha^{2}\lambda_{j}(\alpha\mathbf{x},\alpha^{2}t)\,,\,j=1,2,3\,,$
(4.0a)
$\displaystyle\mathbf{k}^{\sharp}(\mathbf{x},t)=\alpha\mathbf{k}(\alpha\mathbf{x},\alpha^{2}t)\,,\quad\mathbf{U}^{\sharp}(\mathbf{x},t)=\alpha\mathbf{U}(\alpha\mathbf{x},\alpha^{2}t)\,,\quad
J^{\sharp}(\mathbf{x},t)=\alpha^{3}J(\alpha\mathbf{x},\alpha^{2}t)\,.$ (4.0b)

Using (4.1) and (4.1), one can show that the Whitham modulation equations
(3.2) remain unchanged. Similarly, it can be seen that spatial reflections
leave the modulation equations invariant upon the following transformation of
the dependent variables:

$\displaystyle\rho^{\sharp}(\mathbf{x},t)=\rho(-\mathbf{x},t)\,,\quad\lambda^{\sharp}_{j}(\mathbf{x},t)=\lambda_{j}(-\mathbf{x},t)\,,\,j=1,2,3\,,\quad\mathbf{k}^{\sharp}(\mathbf{x},t)=-\mathbf{k}(-\mathbf{x},t)\,,\quad$
(4.0a)
$\displaystyle\mathbf{U}^{\sharp}(\mathbf{x},t)=-\mathbf{U}(-\mathbf{x},t)\,,\quad
J^{\sharp}(\mathbf{x},t)=J(-\mathbf{x},t)\,.$ (4.0b)

Finally, with regards to Galilean transformations, writing
$q^{\sharp}(\mathbf{x},t)=\sqrt{\rho^{\sharp}(\mathbf{x},t)}\,e^{i\phi^{\sharp}(\mathbf{x},t)}$
implies

$\rho^{\sharp}(\mathbf{x},t)=\rho(\mathbf{x}-2\mathbf{w}t,t)\,,\qquad\phi^{\sharp}(\mathbf{x},t)=\phi(\mathbf{x}-2\mathbf{w}t,t)+\mathbf{w}\cdot\mathbf{x}-\|\mathbf{w}\|^{2}t\,.$
(4.1)

The dependent variables in the modulation equations (3.2) become

$\displaystyle\lambda^{\sharp}_{j}(\mathbf{x},t)=\lambda_{j}(\mathbf{x}-2\mathbf{w}t,t)\,,\,j=1,2,3\,,$
(4.0a)
$\displaystyle\mathbf{k}^{\sharp}(\mathbf{x},t)=\mathbf{k}(\mathbf{x}-2\mathbf{w}t,t)\,,\quad\mathbf{U}^{\sharp}(\mathbf{x},t)=\mathbf{U}(\mathbf{x}-2\mathbf{w}t,t)+\mathbf{w}\,,\quad
J^{\sharp}(\mathbf{x},t)=J(\mathbf{x}-2\mathbf{w}t,t)\,.$ (4.0b)

Using (4.1), one can verify that the modulation equations (3.2) remain
invariant under the above Galilean transformation. The Riemann-type variables,
which will be introduced in section 5, change as follows under the above
transformations:

$\displaystyle r_{j}^{\sharp}(\mathbf{x},t)=\alpha
r_{j}(\alpha\mathbf{x},\alpha^{2}t),\quad
r_{j}^{\sharp}(\mathbf{x},t)=r_{j}(-\mathbf{x},t),\quad
r^{\sharp}_{j}(\mathbf{x},t)=r_{j}(\mathbf{x}-2\mathbf{w}t,t)+\mathbf{w}\cdot\hat{\mathbf{k}}/2,\quad
j=1,2,3,4\,.$

###  4.3 Harmonic and soliton limits of the modulation equations in physical
variables

The harmonic and soliton limits of the Whitham equations for the KdV and 1DNLS
equations have proven to be quite useful to study various nonlinear dynamical
scenarios of practical interest [20, 42, 52]. The same is true for the
harmonic and soliton limits of the KP-Whitham equations [14, 48, 50, 49]. We
therefore expect that the same will also be true for the harmonic and soliton
limit of the NLS equation in multiple spatial dimensions.

Like with the periodic solution, the harmonic limit of the Whitham equations
is the limit $m\to 0$, corresponding to $\lambda_{2}\to\lambda_{1}^{+}$.
Recall that in this limit the solution becomes a plane wave. The integrals in
(3.2) simplify considerably:

$\overline{\rho}=\lambda_{1}\,,\quad\overline{\rho^{2}}=\lambda_{1}^{2}\,,\quad\Big{(}\overline{\frac{(\rho^{\prime})^{2}}{\rho}}\Big{)}=0\,,\quad\overline{\rho^{-1}}=1/\lambda_{1}\,,\quad\overline{\rho^{-2}}=1/\lambda_{1}^{2}\,,\qquad
J=\sigma\lambda_{1}\sqrt{\lambda_{3}}\,.$ (4.2)

Then, the linear dispersion relation is

$\omega=2\|\mathbf{k}\|\left(\hat{\mathbf{k}}\cdot\bar{\mathbf{u}}-\sigma\sqrt{\pi^{2}\|\mathbf{k}\|^{2}+\bar{\rho}^{2}}\,\right)\,,$
(4.3)

the averaged energy limits to
$\bar{h}=\bar{\rho}\|\bar{\mathbf{u}}\|^{2}+\bar{\rho}^{2}$ and the Whitham
equations (3.2) reduce to:

$\displaystyle\mathbf{k}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\omega=\mathbf{0}\,,$ (4.0a)
$\displaystyle\bar{\mathbf{u}}_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\big{(}2\bar{\rho}+\|\bar{\mathbf{u}}\|^{2}\big{)}=\mathbf{0}\,,\qquad{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\times\bar{\mathbf{u}}=0\,,$
(4.0b) $\displaystyle\bar{\rho}_{t}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot(\bar{\rho}\bar{\mathbf{u}})=0\,,$
(4.0c) $\displaystyle(\bar{\rho}\bar{\mathbf{u}})_{t}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}(\bar{\rho}^{2})+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\cdot\big{(}\bar{\rho}\bar{\mathbf{u}}\otimes\bar{\mathbf{u}}\big{)}=\mathbf{0}\,,$
(4.0d)
$\displaystyle\bar{h}_{t}+\nabla\cdot\big{(}2(\bar{h}+\bar{\rho}^{2})\bar{\mathbf{u}}\big{)}=0\,.$
(4.0e)

Again, not all of these equations are independent. For example, one can derive
(4.0d) using (4.0b) and (4.0c). Also, note that the variable $\mathbf{k}$ is
immaterial, since its value does not affect the solution, and (4.3) is
decoupled from the other PDEs. Thus, equations (4.0b) and (4.0c), which are
equivalent to the shallow water equations, are by themselves a closed system
of evolution PDEs for the parameters of the plane wave solution, $\bar{\rho}$
and $\bar{\mathbf{u}}$. Nonetheless, (4.0a) describes the evolution of a
harmonic wave propagating on top of the mean flow.

Finally, we discuss the soliton limit of the Whitham modulation system (3.2),
obtained for $m\to 1$ corresponding to $\lambda_{2}\to\lambda_{3}$. In this
limit, the integrals in (3.2) become:

$\overline{\rho}=\lambda_{3}\,,\quad\overline{\rho^{2}}=\lambda_{3}^{2}\,,\quad\|\mathbf{k}\|^{2}\,\overline{\Big{(}\frac{(\rho^{\prime})^{2}}{\rho}\Big{)}}=0\,,\quad\overline{\rho^{-1}}=1/\lambda_{3}\,,\quad\overline{\rho^{-2}}=1/\lambda_{3}^{2}\,.$
(4.1)

Then (3.0a) and (3.0b) are trivially satisfied, and the rest simplify to:

$\displaystyle\bar{\mathbf{u}}_{t}+\nabla\big{(}2\bar{\rho}+\|\bar{\mathbf{u}}\|^{2}\big{)}=\mathbf{0}\,,$
(4.0a)
$\displaystyle\bar{\rho}_{t}+2\nabla\cdot(\bar{\rho}\bar{\mathbf{u}})=0\,,$
(4.0b)
$\displaystyle(\bar{\rho}\bar{\mathbf{u}})_{t}+\nabla(\bar{\rho}^{2})+2\nabla\cdot\big{(}\bar{\mathbf{u}}\otimes\bar{\rho}\bar{\mathbf{u}}\big{)}=\mathbf{0}\,,$
(4.0c)
$\displaystyle\bar{h}_{t}+\nabla\cdot\big{(}2(\bar{h}+\bar{\rho}^{2})\bar{\mathbf{u}}\big{)}=0\,.$
(4.0d)

Note that, similar to before, we can derive equation (4.0c) and (4.0d) from
(LABEL:e:vw1_m=1) and (4.0b). Therefore, we have a system of $N+2$ PDEs for
the dependent variables $\bar{\mathbf{u}}$ and $\bar{\rho}=\lambda_{2}$. But
in this case, we are missing PDEs for $\lambda_{1}$ and $\hat{\mathbf{k}}$
that define the soliton amplitude and its propagation direction, which are
needed to completely determine the soliton solution. This deficiency is also
present in the one-dimensional case. The one-dimensional case is simpler,
however, because, in that case, $\mathbf{k}$ is a one-component vector, and
therefore $\hat{\mathbf{k}}=\pm 1$, constant. The soliton limit is singular so
care must be taken in its calculation. In any case, both in the one-
dimensional and higher-dimensional situation, the problem is eliminated by the
transformation to Riemann-type variables, as we will see later.

##  5 2DNLS-Whitham equations in Riemann-type variables

In this section and the next one we temporarily restrict our attention to the
two-dimensional case and perform suitable changes of dependent variables to
simplify the form of the 2DNLS-Whitham equations.

When $N=2$, the modulation system (3.2) consists of eight PDEs for six
dependent variables in the independent variables $\mathbf{x}=(x,y)^{T}$ and
$t$, plus the two scalar constraints (3.0b) and (3.0d). We will use the four
scalar conservation of waves equations (3.0a) and (3.0c) together with the
averaged conservation of mass (3.0e) and one of the components of the
conservation of momentum equations (3.0f), neglecting the compatibility
conditions (3.0b) & (3.0d) as well as the conservation of energy (3.0g).
Importantly, however, the resulting Whitham equations are equivalent to those
obtained by working with a different set of averaged equations [5].

As in the one-dimensional case, the transformation involves two steps. The
first step is the change of dependent variables from
$(A,k_{1},k_{2},m,\bar{u}_{1},\bar{u}_{2})$ to
$\mathbf{Y}=(\sqrt{\lambda_{1}},\sqrt{\lambda_{2}},\sqrt{\lambda_{3}},U_{1},U_{2},q)$,
with

$q=k_{2}/k_{1}=\tan\varphi\,,$ (5.1)

similar to [2], where $\varphi=\arctan(k_{2}/k_{1})$ [not to be confused with
the fast phase $\phi(Z)$ that was used in sections 2 and 3] identifies the
direction of the periodic wave’s fronts:

$\hat{\mathbf{k}}=(\cos\varphi,\sin\varphi)^{T}\,.$ (5.2)

The second step of the transformation is then defined by the map from
$\lambda_{1},\lambda_{2},\lambda_{3}$ and $U_{1}$ to the “Riemann-type”
variables $\dot{r}_{1},\dot{r}_{2},\dot{r}_{3},\dot{r}_{4}$ via the
transformation

$\displaystyle
U_{1}={\textstyle\frac{1}{2}}\cos\varphi\,(\dot{r}_{1}+\dot{r}_{2}+\dot{r}_{3}+\dot{r}_{4}),$
(5.0a)
$\displaystyle\lambda_{1}={\textstyle\frac{1}{4}}(\dot{r}_{1}-\dot{r}_{2}-\dot{r}_{3}+\dot{r}_{4})^{2},\quad\lambda_{2}={\textstyle\frac{1}{4}}(\dot{r}_{1}-\dot{r}_{2}+\dot{r}_{3}-\dot{r}_{4})^{2},\quad\lambda_{3}={\textstyle\frac{1}{4}}(\dot{r}_{1}+\dot{r}_{2}-\dot{r}_{3}-\dot{r}_{4})^{2}\,.$
(5.0b)

The variables $\dot{r}_{1},\dots,\dot{r}_{4}$ are one possible two-dimensional
generalization of the Riemann invariants of the Whitham equations for the
1DNLS equation. Note that in this work the overdot does not denote
differentiation with respect to time.

Recall that the existence of Riemann invariants for (1+1)-dimensional
hydrodynamic-type systems is intimately tied to the integrability properties
of the modulation equations. Using the one-dimensional Riemann invariants as
dependent variables in higher-dimensional systems diagonalizes their one-
dimensional reductions, and makes the equations more advantageous for analysis
(e.g., see [25]). We will show below that, for both the two-dimensional and
three-dimensional cases, a suitable generalization of the one-dimensional
Riemann invariants allows one to write the modulation equations in a concise
and convenient form.

In terms of $\dot{r}_{1},\dots,\dot{r}_{4}$, the periodic solution (2.7)
becomes

$\displaystyle\rho(Z)={\textstyle\frac{1}{4}}(\dot{r}_{1}-\dot{r}_{2}-\dot{r}_{3}+\dot{r}_{4})^{2}+(\dot{r}_{2}-\dot{r}_{1})(\dot{r}_{4}-\dot{r}_{3})\mathop{\rm
sn}\nolimits^{2}(2K_{m}\,Z|m)\,,$ (5.0a) $\displaystyle
m=\frac{(\dot{r}_{2}-\dot{r}_{1})(\dot{r}_{4}-\dot{r}_{3})}{(\dot{r}_{3}-\dot{r}_{1})(\dot{r}_{4}-\dot{r}_{2})}\,.$
(5.0b)

Moreover,
$\dot{\mathbf{R}}=(\dot{r}_{1},\dot{r}_{2},\dot{r}_{3},\dot{r}_{4},U_{\perp},q)^{T}$
satisfies the hydrodynamic system

$\dot{\mathbf{R}}_{t}+\mathsf{M}_{1}\dot{\mathbf{R}}_{x}+\mathsf{M}_{2}\dot{\mathbf{R}}_{y}=0\,.$
(5.1)

The matrices $\mathsf{M}_{1}$ and $\mathsf{M}_{2}$ are rather complicated, and
we therefore omit them for brevity. When $k_{2}=U_{\perp}=0$, however, the
last two equations in (5.1) are trivially satisfied, and the first four reduce
to the Whitham equations for the 1DNLS equation in Riemann invariant
(diagonal) form [28, 47]:

${\frac{\partial\dot{\mathbf{r}}}{\partial
t}}+\mathsf{V}\,{\frac{\partial\dot{\mathbf{r}}}{\partial x}}=0\,,$ (5.2)

with

$\dot{\mathbf{r}}=(\dot{r}_{1},\dots,\dot{r}_{4})^{T},\qquad\dot{\mathsf{V}}=\mathop{\rm
diag}\nolimits(\dot{\mathbf{V}}),\qquad\dot{\mathbf{V}}=(\dot{V}_{1},\dots,\dot{V}_{4})^{T},\qquad$
(5.3)
$\displaystyle\dot{V}_{1}=2V_{o}+\frac{2(\dot{r}_{2}-\dot{r}_{1})(\dot{r}_{4}-\dot{r}_{1})K_{m}}{(\dot{r}_{4}-\dot{r}_{2})E_{m}-(\dot{r}_{4}-\dot{r}_{1})K_{m}}\,,\quad\dot{V}_{2}=2V_{o}+\frac{2(\dot{r}_{2}-\dot{r}_{1})(\dot{r}_{3}-\dot{r}_{2})K_{m}}{(\dot{r}_{3}-\dot{r}_{2})K_{m}-(\dot{r}_{3}-\dot{r}_{1})E_{m}}\,,$
(5.0a)
$\displaystyle\dot{V}_{3}=2V_{o}+\frac{2(\dot{r}_{3}-\dot{r}_{2})(\dot{r}_{4}-\dot{r}_{3})K_{m}}{(\dot{r}_{4}-\dot{r}_{2})E_{m}-(\dot{r}_{3}-\dot{r}_{2})K_{m}}\,,\quad\dot{V}_{4}=2V_{o}+\frac{2(\dot{r}_{4}-\dot{r}_{1})(\dot{r}_{4}-\dot{r}_{3})K_{m}}{(\dot{r}_{4}-\dot{r}_{1})K_{m}-(\dot{r}_{3}-\dot{r}_{1})E_{m}}\,,$
(5.0b)

with $V_{o}=U_{1}$.

The Whitham modulation system (5.1) can be further simplified by introducing a
modified set of Riemann-type variables:

$r_{j}=\cos\varphi\,\,\dot{r}_{j}\,,\quad j=1,\dots,4\,,$ (5.0a)

with $q=\tan\varphi$ as before. Moreover, the curl-free constraint (3.0d)
yields [see section 7 for details]

$p=\sec\varphi\,\,U_{\perp}\,,$ (5.0b)

where the perpendicular component of $\mathbf{U}$ is defined by

$U_{\perp}=\mathbf{U}\cdot\hat{\mathbf{k}}_{\perp}\,,\qquad\hat{\mathbf{k}}_{\perp}=(-\sin\varphi,\cos\varphi)^{T}\,.$
(5.1)

The Whitham modulation equations (5.1) then reduce to the following form:

${\frac{\partial\mathbf{R}}{\partial
t}}+\mathsf{A}~{}{\frac{\partial\mathbf{R}}{\partial
x}}+\mathsf{B}~{}{\frac{\partial\mathbf{R}}{\partial y}}=0\,,$ (5.2)

where $\mathbf{R}=(r_{1},\dots,r_{4},q,p)^{T}$,

$\mathsf{A}=\left(\begin{array}[]{cc}A_{4\times 4}&A_{4\times 2}\\\ A_{2\times
4}&A_{2\times
2}\end{array}\right)\,,\qquad\mathsf{B}=\left(\begin{array}[]{cc}B_{4\times
4}&B_{4\times 2}\\\ B_{2\times 4}&B_{2\times 2}\end{array}\right)\,,$ (5.0a)

with $g=1+q^{2}$ as in [5] and

$\displaystyle A_{4\times
4}=\mathsf{V}-q^{2}U_{1}\,{\mathbb{1}}_{4}+q^{2}(\mathbf{1}\otimes\mathbf{r}+\mathbf{r}\otimes\mathbf{1})\,,\quad
A_{2\times 2}=2\left(\begin{array}[]{cc}(1-q^{2})U_{1}&-q^{2}\\\
q^{2}(2U_{1}^{2}-s_{2})&gU_{1}\end{array}\right)\,,$ (5.0d) $\displaystyle
A_{4\times
2}=-q\big{(}~{}2U_{1}\mathbf{r}-\mathsf{V}\mathbf{r}+\mathbf{a}/g~{},~{}U_{1}\mathbf{1}-(2(q^{2}-1)\mathbf{r}+\mathbf{V})/g~{}\big{)}\,,~{}~{}A_{2\times
4}=gq\,\big{(}\,-\mathbf{1}~{},~{}2(U_{1}\mathbf{1}-\mathbf{r})~{}\big{)}^{T},$
(5.0e) $\displaystyle B_{4\times
4}=q(\mathsf{V}+U_{1}\,{\mathbb{1}}_{4})+2p\,\mathsf{I}_{4}-q(\mathbf{1}\otimes\mathbf{r}+\mathbf{r}\otimes\mathbf{1})\,,\quad
B_{2\times 2}=2\left(\begin{array}[]{cc}p+2q\,U_{1}&q\\\
-q(2U_{1}^{2}-s_{2})&p\end{array}\right)\,,$ (5.0h) $\displaystyle B_{4\times
2}=1/(2g)\,\big{(}\,2\mathbf{a}~{},~{}(1-q^{2})(4\mathbf{r}-\mathbf{V})\,\big{)}\,,\qquad
B_{2\times 4}=-A_{2\times 4}/q\,,$ (5.0i)

with

$\displaystyle\mathbf{r}=(r_{1},\dots,r_{4})^{T},$ (5.0j)
$\displaystyle\mathsf{V}=\mathop{\rm
diag}\nolimits(\mathbf{V}),\quad\mathbf{V}=(V_{1},\dots,V_{4})^{T},\quad$
(5.0k)
$\displaystyle\mathbf{a}={\textstyle\frac{1}{3}}\big{[}4U_{1}(1-3q^{2})\mathbf{r}-2U_{1}\mathbf{V}-(1+3q^{2})((s_{2}-2U_{1}^{2})\mathbf{1}-\mathsf{V}\mathbf{r})\,\big{]}\,,$
(5.0l)

where $U_{1}=(r_{1}+r_{2}+r_{3}+r_{4})/2$, $V_{1},\dots,V_{4}$ are as in (5.3)
but with $(r_{1},\dots,r_{4})$ instead of $(\dot{r}_{1},\dots,\dot{r}_{4})$,
$\mathbf{1}=(1,\dots,1)^{T}$, $\mathsf{I}_{n}$ is the $n\times n$ identity
matrix, ${\mathbb{1}}_{n}$ denotes the $n\times n$ matrix with all entries
equal to one, and

$s_{n}=r_{1}^{n}+r_{2}^{n}+r_{3}^{n}+r_{4}^{n}\,.$ (5.1)

In component form, the Whitham modulation equations (5.2) are [5]

$\displaystyle{\frac{\partial r_{j}}{\partial t}}+V_{j}{\frac{\partial
r_{j}}{\partial x}}+(qV_{j}+2p){\frac{\partial r_{j}}{\partial
y}}+h_{j}=0,\quad j=1,2,3,4,$ (5.0a) $\displaystyle{\frac{\partial q}{\partial
t}}+2\big{(}g\,U_{1}+pq\big{)}{\frac{\partial q}{\partial
x}}+2\frac{D}{Dy}\big{[}gU_{1}+pq\,\big{]}=0\,,$ (5.0b)
$\displaystyle{\frac{\partial p}{\partial t}}+2g\,U_{1}{\frac{\partial
p}{\partial x}}+2p{\frac{\partial p}{\partial
y}}+\frac{D}{Dy}\big{[}g(s_{2}-2U_{1}^{2})\,\big{]}=0\,,$ (5.0c)

where

$\displaystyle
h_{j}=2q(U_{1}-r_{j})\frac{DU_{1}}{Dy}-\frac{1}{2}q\frac{Ds_{2}}{Dy}+q(V_{j}-2U_{1})\bigg{(}r_{j}{\frac{\partial
q}{\partial x}}+\frac{1}{2}{\frac{\partial p}{\partial
x}}\bigg{)}+\frac{a_{j}}{g}\frac{Dq}{Dy}-\frac{1-q^{2}}{2g}(V_{j}-4r_{j})\frac{Dp}{Dy}\,$
(5.0a)

and $D_{y}$ is the “convective” derivative as in [2]:

$\frac{D}{Dy}={\frac{\partial}{\partial y}}-q{\frac{\partial}{\partial x}}\,.$
(5.0b)

The steps to obtain (5.1) are just a special case of the ones needed to
simplify the Whitham equations for the three-dimensional NLS equation, which
will be discussed in section 7. All the calculations in section 7 can be
trivially reduced to the two-dimentional case by simply taking
$(q_{1},q_{2})=(q,0)$ and $(p_{1},p_{2})=(p,0)$ there. Therefore we omit the
details here for brevity.

Note the necessary compatibility condition for equations (5.1) in which the
initial data is subject to the curl-free constraints
$\nabla\times\bar{\mathbf{u}}=\nabla\times\mathbf{k}=0$, similarly to the KP
equation [2, 5]. In section 7 we will show how these constraints can be
written out explicitly in terms of the Riemann-type variables.

##  6 Further symmetries and reductions of the 2DNLS-Whitham equations

Both of the sets of Riemann-type variables $\dot{\mathbf{R}}$ and $\mathbf{R}$
introduced in section 5 are useful to study further symmetries of the 2DNLS-
Whitham system.

###  6.1 Reduction to the Whitham equations for the radial NLS equation

The Whitham equations for the 2DNLS equation admit a self-consistent reduction
to the Whitham equations for the radial NLS equation, which were recently
derived [4]. To show this, we first perform a change of independent variables
from the Cartesian coordinates $x$ and $y$ to the polar coordinates

$R=\sqrt{x^{2}+y^{2}}\,,\qquad\theta=\arctan(y/x)\,.$ (6.1)

Using the definition of the convective derivative $D_{y}$ in (5.0b), we find

$\frac{Df}{Dy}=\frac{(y-qx)}{R}\,{\frac{\partial f}{\partial
R}}+\frac{(x+qy)}{R^{2}}{\frac{\partial f}{\partial\theta}}\,.$ (6.2)

Equations (5.0b) and (5.0c) in polar coordinates then become, respectively,

$\displaystyle
q_{t}+g\sum_{i=1}^{4}\big{[}(\sin\theta-q\cos\theta)(r_{i})_{R}+\frac{(\cos\theta+q\sin\theta)}{R}(r_{i})_{\theta}\big{]}+2q(\sin\theta-q\cos\theta)p_{R}+\frac{2q}{R}(\cos\theta+q\sin\theta)p_{\theta}$
$\displaystyle\kern
40.00006pt+\big{[}2U_{1}(1-q^{2})\cos\theta+2(p+2qU_{1})\sin\theta\big{]}q_{R}+\big{[}2(p+2qU_{1})\cos\theta-2U_{1}(1-q^{2})\sin\theta\big{]}\frac{q_{\theta}}{R}=0\,,$
(6.0a) $\displaystyle
p_{t}+2g\sum_{i=1}^{4}(r_{i}-U_{1})\big{[}(\sin\theta-q\cos\theta)(r_{i})_{r}+(\cos\theta+q\sin\theta)\frac{(r_{i})_{\theta}}{r}\big{]}+2(gU_{1}\cos\theta+p\sin\theta)p_{R}$
$\displaystyle\kern 40.00006pt+2(p\cos\theta-
gU_{1}\sin\theta)\frac{p_{\theta}}{R}+2q(s_{2}-2U_{1}^{2})\big{[}(\sin\theta-q\cos\theta)q_{R}+(\cos\theta+q\sin\theta)\frac{q_{\theta}}{R}\big{]}=0\,.$
(6.0b)

We then look for a reduction of (6.2) and the remaining four Whitham equations
(5.0a) in which $q=\tan\theta=y/x$. With this assumption, (6.2) simplify
considerably. We also seek solutions in which the Riemann-type variables
$\dot{r}_{1},\dots,\dot{r}_{4}$ are independent of the angular coordinate
$\theta$. Recall that the variables $r_{1},\dots,r_{4}$ appearing in (6.2) are
related to $\dot{r}_{1},\dots,\dot{r}_{4}$ by (5.0a). Thus

${\frac{\partial r_{i}}{\partial
R}}=\frac{1}{\sqrt{g}}{\frac{\partial\dot{r}_{i}}{\partial
R}}-\frac{qr_{i}}{g^{3/2}}{\frac{\partial q}{\partial
R}}\,,\qquad{\frac{\partial
r_{i}}{\partial\theta}}=-\frac{q\dot{r}_{i}}{g^{3/2}}{\frac{\partial
q}{\partial\theta}}\,.$ (6.1)

Substituting the above expression into (6.0a) and (6.0b) yields, respectively,

$\displaystyle p_{\theta}+\cot\theta p=0\,,$ (6.0a) $\displaystyle
p_{t}+2(U_{1}\sec\theta+p\sin\theta)p_{R}+2(p\cos\theta-\tan\theta\sec\theta
U_{1})\,{p_{\theta}}/{R}=0\,.$ (6.0b)

Equation (6.0a) yields $p(R,\theta,t)=C(R,t)\,\csc\theta$, with $C(R,t)$ to be
determined. Then, substituting this expression into (6.0b) yields
$C_{t}+2(U_{1}\sec\theta+C)\,C_{R}-2(C\,\cot^{2}\theta-
U_{1}\sec\theta)\,{C}/{R}=0$, whose only self-consistent solution is $C=0$,
implying $p(R,\theta,t)=0$.

Now we turn our attention to the reduction of the first four Whitham
modulation equations, namely (5.0a). Tedious but straightforward calculations
show that, when written in the polar coordinates (6.1), and using
$q=\tan\theta$ and $p=0$ as well as (6.1), the four modulation equations
(5.0a) become exactly the Whitham equations for the radial NLS equation
derived in [4]:

${\frac{\partial\dot{\mathbf{r}}}{\partial
t}}+\dot{\mathsf{V}}\,{\frac{\partial\dot{\mathbf{r}}}{\partial
R}}+\frac{\dot{\mathbf{b}}}{R}=0\,,$ (6.1)

with $\dot{\mathbf{r}}=(\dot{r}_{1},\dots,\dot{r}_{4})^{T}$ and
$\dot{\mathsf{V}}=\mathop{\rm diag}\nolimits(\dot{\mathbf{V}})$ as in section
5, with $\dot{\mathbf{b}}=(\dot{b}_{1},\dots,\dot{b}_{4})^{T}$,

$\displaystyle\displaystyle\dot{b}_{1}$
$\displaystyle=2V_{o}^{2}-{\textstyle\frac{1}{3}}{(\dot{r}_{2}+\dot{r}_{3}+\dot{r}_{4})V_{1}}-{\textstyle\frac{1}{3}}{[(\dot{r}_{2}+\dot{r}_{3})^{2}+(\dot{r}_{3}+\dot{r}_{4})^{2}+(\dot{r}_{2}+\dot{r}_{4})^{2}]}\,,$
(6.0a) $\displaystyle\dot{b}_{2}$
$\displaystyle=2V_{o}^{2}-{\textstyle\frac{1}{3}}{(\dot{r}_{1}+\dot{r}_{3}+\dot{r}_{4})V_{2}}-{\textstyle\frac{1}{3}}{[(\dot{r}_{1}+\dot{r}_{3})^{2}+(\dot{r}_{3}+\dot{r}_{4})^{2}+(\dot{r}_{1}+\dot{r}_{4})^{2}]}\,,$
(6.0b) $\displaystyle\dot{b}_{3}$
$\displaystyle=2V_{o}^{2}-{\textstyle\frac{1}{3}}{(\dot{r}_{1}+\dot{r}_{2}+\dot{r}_{4})V_{3}}-{\textstyle\frac{1}{3}}{[(\dot{r}_{1}+\dot{r}_{2})^{2}+(\dot{r}_{2}+\dot{r}_{4})^{2}+(\dot{r}_{1}+\dot{r}_{4})^{2}]}\,,$
(6.0c) $\displaystyle\dot{b}_{4}$
$\displaystyle=2V_{o}^{2}-{\textstyle\frac{1}{3}}{(\dot{r}_{1}+\dot{r}_{2}+\dot{r}_{3})V_{4}}-{\textstyle\frac{1}{3}}{[(\dot{r}_{1}+\dot{r}_{2})^{2}+(\dot{r}_{2}+\dot{r}_{3})^{2}+(\dot{r}_{3}+\dot{r}_{1})^{2}]}\,,$
(6.0d)

and
$V_{o}={\textstyle\frac{1}{2}}(\dot{r}_{1}+\dot{r}_{2}+\dot{r}_{3}+\dot{r}_{4})$
as before. In terms of the physical variables, the assumption $q=\tan\theta$
implies that the wavefronts are oriented radially, and the requirement $p=0$
means that the mean flow has no transversal component either, which are both
conditions that are consistent with a radially symmetric reduction.

###  6.2 Harmonic limit and soliton limit of the 2DNLS-Whitham equations in
Riemann-type variables

In section 4 we studied the harmonic limit and the soliton limit of the
modulation equations in physical variables, and we saw that the singular
soliton limit yields fewer equations than are needed to describe the
parameters of the soliton solutions of the 2DNLS equation. We next study the
corresponding limits of the Whitham modulation equations in Riemann-type
variables, and we show how the transformation to Riemann-type variables
eliminates this problem and yields a closed system of equations.

The harmonic limit ($m\to 0$) corresponds to either $r_{2}\to r_{1}^{+}$ or
$r_{3}\to r_{4}^{-}$. In the former case, the PDE (5.0a) with $j=1$ and the
one with $j=2$ coincide, as needed for the limit to be a self-consistent
reduction, and the Whitham modulation system (5.2) then becomes

$\mathbf{R}_{t}+\mathsf{A}_{o.1}\mathbf{R}_{x}+\mathsf{B}_{o.1}\mathbf{R}_{y}=0\,,$
(6.0a)

with $\mathbf{R}=(r_{1},r_{3},r_{4},q,p)^{T}$. The matrices $\mathsf{A}_{o.1}$
and $\mathsf{B}_{o.1}$ are simply the matrices $\mathsf{A}$ and $\mathsf{B}$
from section 5 with $r_{2}=r_{1}$ and the second row and column omitted.
Moreover, the Riemann speeds reduce to

$V_{1}=V_{2}=4r_{1}-\frac{(r_{3}-r_{4})^{2}}{2r_{1}-r_{3}-r_{4}}\,,\qquad
V_{3}=3r_{3}+r_{4}\,,\qquad V_{4}=r_{3}+3r_{4}\,,$ (6.0b)

while $h_{1},\dots,h_{4}$ are still given by (5.0a) with $r_{2}=r_{1}$. In the
latter case (i.e., $r_{3}\to r_{4}^{-}$), the PDE (5.0a) with $j=3$ and the
one with $j=4$ coincide, and the Whitham modulation system (5.2) then becomes

$\mathbf{R}_{t}+\mathsf{A}_{o.2}\mathbf{R}_{x}+\mathsf{B}_{o.2}\mathbf{R}_{y}=0\,,$
(6.0a)

with $\mathbf{R}=(r_{1},r_{2},r_{3},q,p)^{T}$. The matrices $\mathsf{A}_{o.2}$
and $\mathsf{B}_{o.2}$ are just the matrices $\mathsf{A}$ and $\mathsf{B}$
from section 5 with $r_{4}=r_{3}$ and the fourth row and column omitted. The
Riemann speeds reduce to

$V_{1}=3r_{1}+r_{2}\,,\qquad V_{2}=r_{1}+3r_{2}\,,\qquad
V_{4}=V_{3}=4r_{3}+\frac{(r_{1}-r_{2})^{2}}{r_{1}+r_{2}-2r_{3}}\,,$ (6.0b)

with $h_{1},\dots,h_{4}$ now given by (5.0a) with $r_{4}=r_{3}$. In both
cases, it is straightforward to verify that, once the transformation to
Riemann-type variables is inverted and the modulation equations are written
back in terms of the physical variables, one recovers the system (4.3).

The soliton limit ($m\to 1$) corresponds to $r_{3}\to r_{2}^{+}$. In this
case, the PDEs (5.0a) with $j=3$ and the one with $j=2$ coincide, and the
remaining equations become

$\mathbf{R}_{t}+\mathsf{A}_{1}\mathbf{R}_{x}+\mathsf{B}_{1}\mathbf{R}_{y}=0\,,$
(6.0a)

with $\mathbf{R}=(r_{1},r_{2},r_{4},q,p)^{T}$. The matrices $\mathsf{A}_{1}$
and $\mathsf{B}_{1}$ are $\mathsf{A}$ and $\mathsf{B}$ from section 5 with
$r_{3}=r_{2}$ and the fourth row and column omitted. The Riemann speeds reduce
to

$V_{1}=3r_{1}+r_{4}\,,\qquad V_{2}=V_{3}=r_{1}+2r_{2}+r_{4}\,\qquad
V_{4}=r_{1}+3r_{4}\,,$ (6.0b)

where $h_{1},\dots,h_{4}$ are still given by (5.0a) with $r_{3}=r_{2}$. As in
the harmonic limit, it is straightforward to verify that, once the
transformation to Riemann-type variables is inverted and the modulation
equations are written back in terms of the physical variables, one recovers
the system (4.1). In this case, however, the equations in Riemann-type
variables also allow us to obtain the two previously missing modulation
equations, which determine the evolution of $\hat{\mathbf{k}}$ and the soliton
amplitude $a=\lambda_{3}-\lambda_{1}$. One of these equations is immediate,
since (5.0b) directly determines $q=\tan\varphi$ and therefore
$\hat{\mathbf{k}}$. As for the amplitude equation, note that (5.0b) yields
$\lambda_{3}-\lambda_{1}=\sec^{2}\varphi\,(r_{4}-r_{2})(r_{3}-r_{1})$.
Therefore, the modulation equation for $r_{1}$, $r_{2}=r_{3}$, $r_{4}$, and
$q$ determine the evolution of the soliton amplitude and direction.

##  7 Whitham modulation equations for the NLS equation in three spatial
dimensions

We now show how, thanks to the rotation-invariant form of all equations in
sections 2 and 3, the results of section 5 are easily generalized to the NLS
equation in three spatial dimensions.

###  7.1 Set-up and resulting 3DNLS-Whitham system

The Madelung transformation (2.1) yields the same hydrodynamic system of PDEs
(2.1) as well as the mass, momentum and energy conservation laws (2.1) in
differential form, now with $\mathbf{u}=(u_{1},u_{2},u_{3})^{T}$,
$\mathbf{x}=(x,y,z)^{T}$ and ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}=(\partial_{x},\partial_{y},\partial_{z})^{T}$.
The two-phase ansatz (2.2) is also the same, now with
$\mathbf{k}=(k_{1},k_{2},k_{3})^{T}$ and $\mathbf{v}=(v_{1},v_{2},v_{3})^{T}$,
and the curl-free condition (2.3) is now $\nabla\times\mathbf{u}=0$. The only
difference is the number of independent parameters in the periodic solutions:
eight in three spatial dimensions as opposed to six in two spatial dimensions.
The whole derivation in section 3 also remains the same, including the
averaged conservation laws (3.1) and the Whitham modulation equations (3.2),
again the only difference being the number of equations, which in three
dimensions is eleven evolutionary equations.

The first point at which the derivation for the three-dimensional case
diverges from the two-dimensional one is the transformation to Riemann-type
variables. Compared to [5], the process here is made much easier by the
availability of the second conservation of waves equation (3.0c), which allows
us to bypass the averaged conservation of energy, which, in turn, greatly
simplifies the calculations even in the presence of a third spatial dimension.
We begin with the natural generalization of the parametrization (5.2) for
$\hat{\mathbf{k}}$, namely:

$\displaystyle\hat{\mathbf{k}}=(\cos\varphi,\sin\varphi\cos\alpha,\sin\varphi\sin\alpha)^{T}\,.$
(7.0a) $\displaystyle q_{1}=k_{2}/k_{1}=\tan\varphi\,\cos\alpha\,,\quad
q_{2}=k_{3}/k_{1}=\tan\varphi\,\sin\alpha\,,\quad$ (7.0b) $\displaystyle
g=1+q_{1}^{2}+q_{2}^{2}=1/k_{1}^{2}=\sec^{2}\varphi\,.$ (7.0c)

The leading-order part (2.4) of the curl-free condition (2.3) now consists of
three equations. The first two of them are
$k_{1}u_{2}^{\prime}=k_{2}u_{1}^{\prime}$ and
$k_{1}u_{3}^{\prime}=k_{3}u_{1}^{\prime}$, which, when integrated, yield

$\mathbf{u}_{\flat}(Z)=(u_{2}(Z),u_{3}(Z))^{T}=u_{1}(Z)\,\mathbf{q}+\mathbf{p}\,,$
(7.1)

with $\mathbf{q}=(q_{1},q_{2})^{T}$, $\mathbf{p}=(p_{1},p_{2})^{T}$, and
$p_{1}$, $p_{2}$ are additional modulation variables depending on the slow
variables $\mathbf{x}$ and $t$ that appear due to integration in $Z$. For any
three-component vector $\mathbf{w}=(w_{1},w_{2},w_{3})^{T}$ we introduce the
“flat” notation $\mathbf{w}_{\flat}=(w_{2},w_{3})^{T}$, which we use
extensively, to denote the two-component vector comprised of the second and
third components of the vector $\mathbf{w}$. The third equation,
$k_{2}u_{3}^{\prime}=k_{3}u_{2}^{\prime}$ is automatically satisfied. Also,
averaging (7.1), we obtain the two additional relations:

$\displaystyle\bar{\mathbf{u}}_{\flat}=\bar{u}_{1}\mathbf{q}+\mathbf{p}\,,$
(7.0a) $\displaystyle\mathbf{U}_{\flat}=U_{1}\mathbf{q}+\mathbf{p}\,.$ (7.0b)

Similarly, the first component of (3.0b) yields

$\displaystyle u_{1}(Z)=U_{1}+Jk_{1}/(\|\mathbf{k}\|\rho(Z))\,,$ (7.0c)
$\displaystyle\omega=2k_{1}(gU_{1}+\mathbf{q}\cdot\mathbf{p})\,.$ (7.0d)

Finally, we define the Riemann-type variables $r_{1},\dots,r_{4}$ via the same
transformation as in section 5, namely:

$\displaystyle U_{1}={\textstyle\frac{1}{2}}(r_{1}+r_{2}+r_{3}+r_{4}),$ (7.0a)
$\displaystyle\lambda_{1}={\textstyle\frac{1}{4}}\,g(r_{1}-r_{2}-r_{3}+r_{4})^{2},\quad\lambda_{2}={\textstyle\frac{1}{4}}\,g(r_{1}-r_{2}+r_{3}-r_{4})^{2},\quad\lambda_{3}={\textstyle\frac{1}{4}}\,g(r_{1}+r_{2}-r_{3}-r_{4})^{2}.$

Then, in sections 7.2, 7.3 and 7.4 below, we show that the Whitham modulation
equations (3.2) yield the eight-component system of equations

$\displaystyle{\frac{\partial\mathbf{r}}{\partial
t}}+\mathsf{V}\,{\frac{\partial\mathbf{r}}{\partial
x}}+(\mathbf{q}\otimes\mathsf{V}+2\mathbf{p}\otimes\mathsf{I}_{4})\cdot{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\mathbf{r}+\mathbf{h}(\mathbf{r},\mathbf{q},\mathbf{p})=0,$
(7.0a) $\displaystyle{\frac{\partial\mathbf{q}}{\partial
t}}+2(U_{1}+\mathbf{q}\cdot\mathbf{U}_{\flat})\,{\frac{\partial\mathbf{q}}{\partial
x}}+2\mathbf{D}_{\flat}\,(U_{1}+\mathbf{q}\cdot\mathbf{U}_{\flat})=0\,,$
(7.0b) $\displaystyle{\frac{\partial\mathbf{p}}{\partial
t}}+2(U_{1}+\mathbf{q}\cdot\mathbf{U}_{\flat}){\frac{\partial\mathbf{p}}{\partial
x}}+\mathbf{D}_{\flat}\,(g(\tilde{e}_{1}-U_{1}^{2})+\|\mathbf{p}\|^{2})=0\,.$
(7.0c)

Here, as before, $\mathbf{r}=(r_{1},\dots,r_{4})^{T}$, $\mathsf{V}=\mathop{\rm
diag}\nolimits(\mathbf{V})$ with $\mathbf{V}=(V_{1},\dots,V_{4})^{T}$ as in
(5.0k), and the dot product in (7.0a) operates on the two-component vectors to
its left and its right. That is, in component form, for each $j=1,\dots,4$ the
third term in (7.0a) is the dot product between
$\mathbf{q}\,V_{j}+2\mathbf{p}$ and ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}r_{j}$. Additionally,
(7.0b) and (7.0c) contain the three-dimensional generalization of the
convective derivative of [2] and section 5, namely:

$\mathbf{D}_{\flat}=(D_{y},D_{z})^{T}={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}-\mathbf{q}\,\partial_{x}\,,$
(7.0a)

where ${\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}=(\partial_{y},\partial_{z})^{T}$
and

$\displaystyle\frac{D}{Dy}={\frac{\partial}{\partial
y}}-q_{1}{\frac{\partial}{\partial
x}}\,,\qquad\frac{D}{Dz}={\frac{\partial}{\partial
z}}-q_{2}{\frac{\partial}{\partial x}}\,.$ (7.0b)

The term
$\mathbf{h}(\mathbf{r},\mathbf{q},\mathbf{p})=(h_{1},\dots,h_{4})^{T}$ in
(7.0a) is given by

$\displaystyle
h_{j}=2(U_{1}-r_{j})\mathbf{q}\cdot\mathbf{D}_{\flat}U_{1}-{\textstyle\frac{1}{2}}\mathbf{q}\cdot\mathbf{D}_{\flat}s_{2}+(V_{j}-2U_{1})\,\mathbf{q}\cdot\Big{(}r_{j}{\frac{\partial\mathbf{q}}{\partial
x}}+\frac{1}{2}{\frac{\partial\mathbf{p}}{\partial
x}}\Big{)}-{\textstyle\frac{1}{4}}(V_{j}-4r_{j})\,\mathbf{D}_{\flat}\cdot\mathbf{p}$
$\displaystyle\kern
20.00003pt+a_{j}\mathbf{D}_{\flat}\cdot\mathbf{q}+(b_{j}/g)\,\tr[(\mathbf{q}\otimes\mathbf{q})(\mathbf{D}_{\flat}\otimes\mathbf{q})]+((V_{j}-4r_{j})/g)\,\tr[(\mathbf{q}\otimes\mathbf{q})(\mathbf{D}_{\flat}\otimes\mathbf{p})],$
(7.0a) $\displaystyle
a_{j}={\textstyle\frac{1}{3}}[2(2r_{j}-V_{j})U_{1}-s_{2}+2U_{1}^{2}+V_{j}r_{j}],\quad
b_{j}=r_{j}(V_{j}-4U_{1})-s_{2}+2U_{1}^{2}+a_{j},$ (7.0b)

for $j=1,\dots,4$. The $s_{n}$ are as in (5.1), and
$\tilde{e}_{1}=g(\lambda_{1}+\lambda_{2}+\lambda_{3})$ in (7.0c), is similar
to (2.12). Equations (7.0a), (7.0b), (7.0c) and (7.0a) should be compared to
(5.0a), (5.0b), (5.0c) and (5.0a) in the two-dimensional case. Note that,
while $\mathbf{h}_{4}(\mathbf{r},\mathbf{q},\mathbf{p})$ might give the
impression of a forcing term in (7.0a), that is not the case in reality, as
(7.0a) shows that $\mathbf{h}_{4}(\mathbf{r},\mathbf{q},\mathbf{p})$ is in
fact a homogenous first-order differential polynomial in $\mathbf{r}$,
$\mathbf{q}$ and $\mathbf{p}$, and therefore (7.1) is indeed a system of PDEs
of hydrodynamic type like its one-dimensional and two-dimensional
counterparts.

Similarly to the two-dimensional case, the 3DNLS-Whitham modulation equations
(7.1) are subject to the compatibility conditions
$\nabla\times\bar{\mathbf{u}}(\mathbf{x},0)=\mathbf{0}$ and
$\nabla\times\mathbf{k}(\mathbf{x},0)=\mathbf{0}$ at $t=0$. In Appendix A.3 we
show that, in terms of the dependent variables defined above, these
constraints become, respectively,

$k_{1}\mathbf{q}_{x}={\mathbf{D}_{\flat}k_{1}}\,,\qquad
k_{1}\mathbf{p}_{x}=2(({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}{k_{1}})^{T}\mathsf{R}_{4})\cdot\mathbf{D}_{\flat}\mathbf{r}-2U_{1}{\mathbf{D}_{\flat}k_{1}}\,,$
(7.1)

where ${\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}=(\partial_{r_{1}},\dots,\partial_{r_{4}})^{T}$,
$\mathsf{R}_{4}=\mathop{\rm diag}\nolimits(r_{1},\dots,r_{4})$, and the dot
product operates on the four-component vectors to its left and its right.
Equations (7.1) are conditions that must be satisfied by the initial
conditions for (7.1) in order for its solutions to represent modulations of
actual one-phase solutions of the NLS equation (1.1).

###  7.2 Derivation of the 3DNLS-Whitham system: Equations for the auxiliary
variables

To derive the evolution equation (7.0b) for $\mathbf{q}$, we split the first
conservation of waves equation (3.0a) and rewrite it using the convective
derivatives $D_{y}$ and $D_{z}$ defined in (7.0b), to obtain

$k_{1,t}+\omega_{x}=0\,,\qquad{\mathbf{q}}_{t}+(\mathbf{D}_{\flat}\omega)/k_{1}=0\,,\qquad{\mathbf{q}}_{x}=(\mathbf{D}_{\flat}k_{1})/k_{1}\,,$
(7.2)

with $\omega$ as in (7.0d). The first of equations (7.2) will be used later to
derive (7.0a). Substituting the third equation in (7.2) into the second one
and using (7.0d) yields the desired evolution equation (7.0b). Note also that
the third equation in (7.2) is precisely the first of the constraints (7.1).

Next, to derive the evolution equation (7.0c) for $\mathbf{p}$, we start with
the constraint (3.0d) for the second conservation of waves equation. Using
(7.1), (3.0d) yields

$(\bar{u}_{1}\mathbf{q}+\mathbf{p})_{x}={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\bar{u}_{1}\,.$ (7.3)

Averaging (7.0c) over the unit period, we can rewrite the above as

$\partial_{x}\big{[}(U_{1}+J\,\overline{\rho^{-1}}/g^{1/2})\mathbf{q}+\mathbf{p}\big{]}={\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\big{(}U_{1}+J\,\overline{\rho^{-1}}/g^{1/2}\big{)}\,,$
(7.4)

and simplifying further we obtain

${\mathbf{p}}_{x}=\mathbf{D}_{\flat}U_{1}+\mathbf{D}_{\flat}\big{(}J\,\overline{\rho^{-1}}/g^{1/2}\big{)}-\big{(}U_{1}+J\,\overline{\rho^{-1}}/g^{1/2}\big{)}(\mathbf{D}_{\flat}k_{1})/{k_{1}}\,.$
(7.5)

Now we use the second conservation of waves (3.0c), written in the form of
(3.0c). From the second and third components, together with the above
relations, we have

$\displaystyle\partial_{t}\big{(}\big{(}U_{1}+J\,\overline{\rho^{-1}}/g^{1/2}\big{)}\mathbf{q}+\mathbf{p}\big{)}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\big{(}g(U_{1}^{2}+\tilde{e}_{1})+\|\mathbf{p}\|^{2}+2U_{1}\mathbf{q}\cdot\mathbf{p}+2(gU_{1}+\mathbf{q}\cdot\mathbf{p})J\,\overline{\rho^{-1}}/g^{1/2}\big{)}=0\,.$

Simplifying yields

$\displaystyle{\mathbf{p}}_{t}+\mathbf{D}_{\flat}\big{(}g(U_{1}^{2}+\tilde{e}_{1})+\|\mathbf{p}\|^{2}+2U_{1}\mathbf{p}\cdot\mathbf{q}\big{)}-2U_{1}\,\mathbf{D}_{\flat}(gU_{1}+\mathbf{p}\cdot\mathbf{q})$
$\displaystyle\kern
100.00015pt+2(gU_{1}+\mathbf{p}\cdot\mathbf{q})\big{(}\mathbf{D}_{\flat}\big{(}J\,\overline{\rho^{-1}}/g^{1/2}\big{)}-\big{(}U_{1}+J\,\overline{\rho^{-1}}/g^{1/2}\big{)}(\mathbf{D}_{\flat}k_{1})/{k_{1}}\big{)}=0\,.$

Using (7.5) yields the desired equation (7.0c).

###  7.3 Derivation of the 3DNLS-Whitham system: Convective derivatives

It remains to derive the four equations in (7.0a) for the Riemann-type
variables $r_{1},\dots,r_{4}$. To this end, we use the two conservation of
waves equations (3.0a) and (3.0c) (as well as the compatibility conditions
(3.0b) and (3.0d)) along with the averaged conservation of mass and momentum
equations (3.0e), (3.0f). The process comprises three main steps.

The first step is the further simplification of the averaged conservation
laws. Note that

$\displaystyle k_{1}^{2}=g(\lambda_{3}-\lambda_{1})/4K_{m}^{2}\,.$ (7.8)

For convenience we also introduce the quantity
$\mathbf{M}=\overline{\rho\mathbf{u}}/g=(M_{1},M_{2},M_{3})^{T}$, with

$M_{1}:=\overline{\rho
u_{1}}/g=U_{1}\overline{\rho}/g+\tilde{J}\,,\qquad\mathbf{M}_{\flat}:=M_{1}\mathbf{q}+(\overline{\rho}/g)\,\mathbf{p}=(\overline{\rho}/g)\mathbf{U}_{\flat}+\tilde{J}\mathbf{q}\,,$
(7.9)

and $\tilde{J}=\hat{k}_{1}\,J/g$. Then, using (7.0d) one can rewrite the
modulation equations (3.0a) and (3.0e) as follows:

$\displaystyle k_{1,t}+2[k_{1}(U_{1}+\mathbf{q}\cdot\mathbf{U}_{\flat})]=0\,,$
(7.10)
$\displaystyle(\overline{\rho})_{t}+2(gM_{1})_{x}+2{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot(g~{}\mathbf{M}_{\flat})=0\,,$
(7.11)

while the first component of the second conservation of waves equation (3.0c)
becomes

$(U_{1}+J\,\overline{\rho^{-1}}/g^{1/2})_{t}+\big{[}g\big{(}\tilde{e}_{1}+U_{1}^{2}+2J\,\overline{\rho^{-1}}/g^{1/2}U_{1}\big{)}+\|\mathbf{p}\|^{2}+2(U_{1}+J\,\overline{\rho^{-1}}/g^{1/2})\mathbf{p}\cdot\mathbf{q}\big{]}_{x}=0\,.$

Moreover, using the equation (A.0b) we can write the averaged momentum
equation (3.0f) in component form as

$\displaystyle(gM_{1})_{t}+\big{[}g\big{(}2U_{1}(M_{1}+\tilde{J})+\tilde{e}_{2}-\overline{\rho}^{2}/g^{2}\big{)}+\overline{\rho}^{2}\big{]}_{x}$
$\displaystyle\kern 100.00015pt+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\big{[}2g(M_{1}+\tilde{J})\mathbf{U}_{\flat}-2g\tilde{J}\mathbf{p}+g(\tilde{e}_{2}-\overline{\rho}^{2}/g^{2})\mathbf{q}\big{]}=0\,,$
(7.0a)
$\displaystyle(g\mathbf{M}_{\flat})_{t}+\Big{[}g\big{(}2(M_{1}+\tilde{J})\mathbf{U}_{\flat}+(\tilde{e}_{2}-\overline{\rho}^{2}/g^{2})\mathbf{q}-2\tilde{J}\mathbf{p}\big{)}\Big{]}_{x}+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\overline{\rho}^{2})$
$\displaystyle\kern 100.00015pt+{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\Big{[}2g\mathbf{M}_{\flat}\otimes\mathbf{U}_{\flat}+2g\tilde{J}\mathbf{U}_{\flat}\otimes\mathbf{q}+g(\tilde{e}_{2}-\overline{\rho}^{2}/g^{2})\,\mathbf{q}\otimes\mathbf{q}\Big{]}=\mathbf{0}\,.$
(7.0b)

Next, we perform a second, intermediate step to write the Whitham modulation
equations in terms of convective derivatives. First, we derive some identities
that will be useful later. Equation (3.0b) and the definition of $\mathbf{q}$
in (7.0b) yield

$\mathbf{q}_{x}=\frac{1}{k_{1}}\mathbf{D}_{\flat}k_{1}\,,\quad\mathbf{q}_{y}=\frac{1}{k_{1}}\mathbf{D}_{\flat}(k_{1}q_{1})\,,\quad\mathbf{q}_{z}=\frac{1}{k_{1}}\mathbf{D}_{\flat}(k_{1}q_{2})\,.$
(7.1)

Moreover, in Appendix A.3, we show that these relations also yield the two
constraints

$\displaystyle D_{y}q_{2}=D_{z}q_{1}\,,$ (7.0a) $\displaystyle
D_{y}p_{2}=D_{z}p_{1}\,,$ (7.0b)

which will prove to be useful. We then define the additional convective
derivatives

$D_{x}={\frac{\partial}{\partial x}}+\mathbf{q}\cdot{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\,,\qquad
D_{t}={\frac{\partial}{\partial t}}+2U_{1}{\frac{\partial}{\partial
x}}+2\mathbf{U}_{\flat}\cdot{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\,.$ (7.1)

Now we rewrite the evolution equations for $\mathbf{q}$ using these convective
derivatives. Specifically, in Appendix A.3 we show that (7.10) and (7.3)
yield, respectively,

$\displaystyle\frac{D_{t}k_{1}}{k_{1}}+2D_{x}U_{1}+W_{1}=0\,,$ (7.2)
$\displaystyle
D_{t}(U_{1}+J\,\overline{\rho^{-1}}/g^{1/2})+(U_{1}-J\,\overline{\rho^{-1}}/g^{1/2})\frac{D_{t}k_{1}}{k_{1}}+D_{x}(\tilde{e}_{1}+U_{1}^{2})+2W_{2}=0\,,$
(7.3)

where

$\displaystyle
gW_{1}=\mathbf{q}\cdot[D_{t}\mathbf{q}+2U_{1}\,D_{x}\mathbf{q}+2D_{x}\mathbf{p}]\,,$
(7.0a) $\displaystyle
gW_{2}=\mathbf{q}\cdot\big{[}U_{1}\,D_{t}\mathbf{q}+s_{2}D_{x}\mathbf{q}+\frac{1}{2}D_{t}\mathbf{p}+U_{1}\,D_{x}\mathbf{p}\big{]}\,.$
(7.0b)

Moreover, in Appendix A.3 we also show that the conservation of mass equation
(7.11) and conservation of momentum equation (7.0a) yield, respectively,

$\displaystyle
g\Big{[}D_{t}(\overline{\rho}/g)-(\overline{\rho}/g)\frac{D_{t}k_{1}}{k_{1}}+2D_{x}\tilde{J}\Big{]}+(\overline{\rho}/g)(\mathbf{q}\cdot
D_{t}\mathbf{q})+6\tilde{J}\mathbf{q}\cdot D_{x}\mathbf{q}$
$\displaystyle\kern 80.00012pt+2M_{1}\big{(}g{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{q}-\mathbf{q}\cdot
D_{x}\mathbf{q}\big{)}+2(\overline{\rho}/g)\big{(}g{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{p}-\mathbf{q}\cdot
D_{x}\mathbf{p}\big{)}=0\,,$ (7.1) $\displaystyle
g\Big{[}(\overline{\rho}/g)D_{t}U_{1}+D_{t}\tilde{J}-2\tilde{J}\frac{D_{t}k_{1}}{k_{1}}+D_{x}\tilde{e}_{2}\Big{]}+\mathbf{q}\cdot[M_{1}D_{t}\mathbf{q}+(\overline{\rho}/g)D_{t}\mathbf{p}+4\tilde{e}_{2}D_{x}\mathbf{q}]$
$\displaystyle\kern
80.00012pt+(\tilde{e}_{2}-(\overline{\rho}^{2}/g^{2})+2U_{1}\tilde{J})\big{(}g{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{q}-\mathbf{q}\cdot
D_{x}\mathbf{q}\big{)}+2\tilde{J}\big{(}g{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{p}-\mathbf{q}\cdot
D_{x}\mathbf{p}\big{)}=0\,.$ (7.2)

Equations (7.2), (7.3), (7.1) and (7.2) comprise the four modified modulation
equations written in terms of the variables $U_{1}$, $\lambda_{1}$,
$\lambda_{2}$ and $\lambda_{3}$ and the convective derivatives $D_{x}$ and
$D_{t}$.

###  7.4 Derivation of the 3DNLS-Whitham system: Equations for Riemann-type
variables

The third and final step in the derivation of (7.0a) is to express the
modulation equations in terms of $r_{1},\dots,r_{4}$. Recall the
transformation (7.1) to the Riemann-type variables. Note that the arrangement
of indices in (7.1) is dictated by the requirement that the constraint (2.11)
be satisfied when $r_{1}\leqslant r_{2}\leqslant r_{3}\leqslant r_{4}$, since

$\displaystyle\lambda_{2}-\lambda_{1}=g(r_{4}-r_{3})(r_{2}-r_{1})\,,~{}~{}$
(7.0a)
$\displaystyle\lambda_{3}-\lambda_{1}=g(r_{4}-r_{2})(r_{3}-r_{1})\,,~{}~{}$
(7.0b) $\displaystyle\lambda_{3}-\lambda_{2}=g(r_{4}-r_{1})(r_{3}-r_{2})\,.$
(7.0c)

In Appendix A.3, using the above definitions, we show that (7.2), (7.3), (7.1)
and (7.2) yield, respectively,

$\displaystyle({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}[\,\log{k_{1}}])^{T}D_{t}\mathbf{r}+D_{x}s_{1}+W_{1}$
$\displaystyle=$ $\displaystyle 0\,,$ (7.0a) $\displaystyle
2({\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}[\,\log{k_{1}}])^{T}\mathsf{R}_{4}\,D_{t}\mathbf{r}+D_{x}s_{2}+2W_{2}$
$\displaystyle=$ $\displaystyle 0\,,$ (7.0b) $\displaystyle
3({\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}[\,\log{k_{1}}])^{T}\mathsf{R}_{4}^{2}D_{t}\mathbf{r}+D_{x}s_{3}+3W_{3}$
$\displaystyle=$ $\displaystyle 0\,,$ (7.0c) $\displaystyle
4({\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}[\,\log{k_{1}}])^{T}\mathsf{R}_{4}^{3}\,D_{t}\mathbf{r}+D_{x}s_{4}+4W_{4}$
$\displaystyle=$ $\displaystyle 0\,,$ (7.0d)

where $\mathbf{r}=(r_{1},\dots,r_{4})^{T}$, ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\mathbf{r}}=(\partial_{r_{1}},\dots,\partial_{r_{4}})^{T}$
and $\mathsf{R}_{4}=\mathop{\rm diag}\nolimits(r_{1},\dots,r_{4})$ as before,
with $W_{1}$ and $W_{2}$ as in (7.0a) and (7.0b), and

$\displaystyle
gW_{3}={\textstyle\frac{1}{4}}(s_{2}-2U_{1}^{2})gW_{1}+U_{1}gW_{2}+{\textstyle\frac{1}{2}}\mathbf{q}\cdot[(\overline{\rho}/g)D_{t}\mathbf{q}+6\tilde{J}D_{x}\mathbf{q}]$
$\displaystyle\kern
40.00006pt+(U_{1}\bar{\rho}/g+\tilde{J})\big{(}g({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{q})-\mathbf{q}\cdot
D_{x}\mathbf{q}\big{)}+(\overline{\rho}/g)\big{(}g({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{p})-\mathbf{q}\cdot
D_{x}\mathbf{p}\big{)}\,,$ (7.0a) $\displaystyle
gW_{4}={\textstyle\frac{1}{8}}(6\tilde{J}-U_{1}s_{2}+2U_{1}^{3})gW_{1}+{\textstyle\frac{1}{4}}(s_{2}-4U_{1}^{2})gW_{2}+{\textstyle\frac{3}{2}}U_{1}gW_{3}$
$\displaystyle\kern
40.00006pt+{\textstyle\frac{1}{4}}\big{[}\mathbf{q}\cdot(M_{1}D_{t}\mathbf{q}+(\overline{\rho}/g)D_{t}\mathbf{p}+4\tilde{e}_{2}D_{x}\mathbf{q})$
$\displaystyle\kern
40.00006pt+(\tilde{e}_{2}-(\overline{\rho}^{2}/g^{2})+2U_{1}\,\tilde{J})\big{(}g({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{q})-\mathbf{q}\cdot
D_{x}\mathbf{q}\big{)}+2\tilde{J}\big{(}g({\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}\cdot\mathbf{p})-\mathbf{q}\cdot
D_{x}\mathbf{p}\big{)}\big{]}\,.$ (7.0b)

Importantly, note that, even though the second conservation of waves equation
(7.3) contains the third complete elliptic integral $\Pi(\,\cdot\,,m)$ via
$\overline{\rho^{-1}}$ [cf. (3.0b)], the third elliptic integral does not
appear in the resulting modulation equation (7.0a). Note that $\Pi(\cdot,m)$
is also contained in the conservation of energy equation. Next, one can
collect the four equations (7.4) and rewrite them in matrix form as

$\mathsf{M}(\mathbf{r})\,\big{(}\nabla_{\mathbf{r}}[\log k_{1}]\cdot
D_{t}\mathbf{r}+D_{x}\mathbf{r}\big{)}+\mathbf{W}=\mathbf{0}\,,$ (7.1)

where $\mathbf{W}=(W_{1},\cdots,W_{4})^{T}$ and $\mathsf{M}(\mathbf{r})$ is
the Vandermonde matrix

$\mathsf{M}(\mathbf{r})=\left(\begin{array}[]{cccc}1&1&1&1\\\
r_{1}&r_{2}&r_{3}&r_{4}\vspace{1mm}\\\
r_{1}^{2}&r_{2}^{2}&r_{3}^{2}&r_{4}^{2}\vspace{1mm}\\\
r_{1}^{3}&r_{2}^{3}&r_{3}^{3}&r_{4}^{3}\end{array}\right)\,.$ (7.2)

Multiplying (7.1) by $\mathsf{M}^{-1}(\mathbf{r})$, we then finally obtain
(7.0a), with

$\displaystyle h_{j}=\frac{(-1)^{j+1}\Delta_{ilm}}{|\Delta|\,(\partial
k/\partial
r_{j})/k}[r_{i}r_{l}r_{m}W_{1}-(r_{i}r_{l}+r_{l}r_{m}+r_{m}r_{i})W_{2}+(r_{i}+r_{l}+r_{m})W_{3}-W_{4}],~{}~{}j=1,\dots,4\,,$
(7.3)

where $j\neq i,j\neq l,j\neq m$, $i<l<m$, summation of repeated indices is
implied, and

$|\Delta|=\prod_{j>l}^{4}(r_{j}-r_{l})\,,\quad\Delta_{ilm}=(r_{i}-r_{l})(r_{l}-r_{m})(r_{m}-r_{i})\,.$
(7.4)

Finally, using equations (A.0a) and (A.0b), one can simplify
$h_{1},\dots,h_{4}$ in (7.3) to obtain (7.0a).

##  8 Discussion and perspectives

In summary, we derived the Whitham modulation equations for the defocusing NLS
equation in two, three and higher spatial dimensions using a two-phase ansatz
and the averaged conservation laws of the NLS equation written in coordinate-
free vector form, and we elucidated various symmetries and reductions of the
resulting equations, including the reduction to the Whitham equations of the
radial NLS equation as well as the harmonic and soliton limits. We point out
that, long after this work was completed, we learned that modulation equations
for multi-dimensional equations of NLS type were written down in physical
variables using a general framework in [7], and the modulation equations were
used to study the stability of the plane wave solutions. On the other hand, no
transformation to Riemann-type variables was carried out in [7].

We reiterate that the use of a two-phase ansatz in this work (as opposed to a
one-phase ansatz as in [5]) greatly simplifies the derivation, since it
results in a second conservation of waves equation that allows us to avoid
using the conservation of energy equation, which is much more complicated in
comparison. Moreover, the advantage of using a two-phase ansatz increases with
the number of spatial dimensions. This is because the number of modulation
equations needed is $2N+2$. Therefore, if one tried to derive the modulation
equations in three spatial dimensions with a one-phase ansatz, one would need
to use additional conservation laws for the NLS equation. This would not only
lead to a much more complicated derivation, but one would quickly exhaust the
number of available conservation laws, since the NLS equation in more than one
spatial dimensions is not completely integrable, and therefore does not have
hidden symmetries resulting in an infinite number of conservation laws.

In contrast, the results of section 7 can be generalized in a straightforward
way to obtain the Whitham modulation equations in simplified form in an
arbitrary number of spatial dimensions. The system of modulation equations
(7.1) is already written in vectorial, dimension-independent form, with the
only caveat that, with $N$ spatial dimensions, $\mathbf{q}$ and $\mathbf{p}$
have $N-1$ components. Moreover, all the steps of the derivation in section 7
are written in a way that generalizes to any number of spatial dimensions.
Indeed, one can introduce spherical coordinates in $N$ spatial dimensions by
generalizing (7.0a) as $\hat{k}_{1}=\cos\varphi_{1}$,
$\hat{k}_{2}=\sin\varphi_{1}\cos\varphi_{2}$,
$\hat{k}_{3}=\sin\varphi_{1}\sin\varphi_{2}\cos\varphi_{3}$, etc., up to
$\hat{k}_{N-1}=\sin\varphi_{1}\cdots\sin\varphi_{N-2}\cos\varphi_{N-1}$ and
$\hat{k}_{N}=\sin\varphi_{1}\cdots\sin\varphi_{N-2}\sin\varphi_{N-1}$. Then,
one introduces $q_{1},\dots,q_{N}$ via the generalization of (7.0b), namely,
$q_{1}=k_{2}/k_{1}$, $q_{2}=k_{3}/k_{1}$, etc., up to $q_{N-1}=k_{N}/k_{1}$,
as well as $p_{1},\dots,p_{N}$ via the natural generalization of (7.1). In
this way, one obtains the generalization of (7.0c) as
$g=1+q_{1}^{2}+\cdots+q_{N}^{2}=\sec\varphi_{1}$, and all the calculations and
equations in section 7 remain valid as long as one also redefines the
operators ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}_{\flat}$ and
$\mathbf{D}_{\flat}$ accordingly.

We point out that, even though we have not done so explicitly in this work, it
would be straightforward to obtain the reduction to the Whitham equations for
the radial NLS equation from the 3DNLS-Whitham equations derived in section 7
using spherical coordinates. It would also be straightforward to write down
explicitly the harmonic and soliton limits in three spatial dimensions, as
well as all of these corresponding limits in higher dimensions.

We should comment on the importance of the constraints ${\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\mathbf{k}=0$ and
${\kern-0.83337pt\lower 0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower
0.20004pt\hbox{\kern 0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise
0.20004pt\hbox{\kern-0.50003pt\hbox{$\nabla$}}}\wedge\bar{\mathbf{u}}=0$. On
one hand, these constraints play a key role in the derivation. On the other
hand, the final Whitham equations [e.g., (7.0a), (7.0b) and (7.0c)] do not
automatically ensure that these constraints are satisfied, only that their
time derivative is (similar to [2]). Because of this, the compatibility
between solutions of the Whitham system and modulated one-phase solutions of
the NLS equation is not guaranteed a priori, and, similar to [2], one must
give initial conditions that are compatible with the one-phase assumption.
These constraints are also likely to be related to the integrability
properties of the system, as discussed below.

We emphasize that this work is foundational and that, similar to [2], the
results presented here open up a number of interesting problems, which are
expected to lead to several further advances in the near future. Specifically,
we next mention and briefly discuss some of these possible of avenues for
further research.

One direction for future work is the derivation of the Whitham equations for
the focusing NLS equation in three spatial dimensions. We expect that this
will be straightforward. Indeed, the Whitham equations in the two-dimensional
focusing case were already written in [5] (although not in rotation-invariant
form). Once the derivation of the one-phase solutions of the NLS equation is
done in dimension-invariant form, as was the case in section 2.2, the rest of
the machinery presented in this work will carry over to the case of the
focusing case in three and higher dimensions without significant changes. Of
course, as in the one-dimensional case, the resulting Whitham equations will
be elliptic (i.e., the characteristic velocities will be complex), and
therefore require suitable interpretation of initial value problems; see [15,
23, 32, 33, 34, 40] as well as [13, 25] and references therein. The Whitham
equations for the NLS equation in one spatial dimension have also proved to be
useful in some situations, even in the focusing case [11, 16, 24], so one can
expect that those in two and three spatial dimensions will be useful as well.

Another important direction for future work is a study to determine whether
the Whitham modulation system derived here, or any of its reductions, are
completely integrable. A notion of integrability for multidimensional systems
was put forth in [26, 27], based on the existence of infinitely many
$N$-component reductions. Of course, the NLS equation in more than one spatial
dimension is not integrable, and therefore one would have no reason to expect
that the corresponding NLS-Whitham systems are. Still, the reductions to one-
dimensional NLS-Whitham equations are indeed integrable, and therefore it is a
natural question whether there are other integrable reductions. In this
regard, we should point out that, even for the KP equation (which is
integrable), the original Whitham system derived in [2] appears not to be
integrable, but its harmonic and soliton limits are [14]. Moreover, so are
various less-trivial one-dimensional reductions beyond the obvious reduction
to the Whitham system for the KdV equation, once one properly takes into
account the analogue of the compatibility conditions (3.0a) and (3.0b) [12].

Yet another interesting problem for future work is the issue of whether one
can establish a precise relation between the 2DNLS-Whitham system and the KP-
Whitham system. It is well known that the 1DNLS-Whitham system admits a
reduction to the KdV-Whitham system [32]. It is also well known that the 2DNLS
equation admits a reduction to the KP equation [41]. A natural question is
therefore whether the 2DNLS-Whitham system admits a reduction to the KP-
Whitham system. It is straightforward to see that, if one considers the same
reduction as in [32], the PDEs for $r_{1},\dots,r_{4}$ in the 2DNLS-Whitham
system naturally reduce to those for $r_{1},\dots,r_{3}$ in the KP-Whitham
system. The PDE for $q$ also reduces to the corresponding equation in the KP-
Whitham system, since it just comes from the second component of the
conservation of waves equation in both systems. The open question, however, is
how one can obtain a PDE for $p$ that does not contain a time derivative, as
prescribed in the KP-Whitham system.

Finally, and most importantly from a practical point of view, an obvious
opportunity for future work will be the use of the modulation equations
derived here to characterize the dynamical behavior in physically significant
scenarios. One important application is to the description of dispersive shock
waves (DSWs) [32, 13]. Some of the earliest experiments on DSWs in nonlinear
optics and Bose-Einstein condensates (BECs)—where the defocusing NLS equation
is an excellent model—involved inherently multidimensional nonlinear wave
propagation [22, 51, 34, 55]. One intriguing feature, observed in both BEC and
optics [34, 55], is the coherent propagation of multidimensional DSWs with
stable ring/spherical and elliptical/ellipsoidal patterns. These observations
are at odds with the known transverse instability of planar cnoidal wave
solutions of (1.1) [54]. Further analysis of the 2D and 3DNLS-Whitham
modulation equations may provide some analytical insight in this. Moreover,
BECs are three-dimensional, so the (3+1)-dimensional modulation equations
derived here are needed to describe large amplitude matter waves. Three-
dimensional effects have been shown to be decisive in some BEC DSW experiments
[19, 43].

Various applications of the Whitham equations for the focusing and defocusing
NLS equations in one spatial dimension were already mentioned above. We should
also note that, while the full modulation system composed of equations (7.0a),
(7.0b) and (7.0c) might appear complicated, even its reductions can be useful
in this regard. For example, of particular interest from an applicative point
of view are the harmonic and soliton limits. In the one-dimensional case,
soliton modulation theory and its applications were studied for the KdV
equation in [42] and for the defocusing NLS equation in [52], while the
harmonic limit of the Whitham equations for the KdV equation was studied in
[20]. Similarly, the harmonic and soliton limits of the Whitham equations for
the KP equation, which were derived and analyzed in [2, 14] have found
concrete applications in [48, 50, 49]. These reductions analytically describe
the evolution of a soliton or linear waves in the presence of the slowly
varying mean field $\bar{\rho}$, $\bar{\mathbf{u}}$. Obtaining these
modulation equations using multiple scales and a soliton ansatz is quite
tedious and, to our knowledge, has apparently only been carried out for the
KdV equation in [31]. We believe that, like with Whitham equations for the KP
equation [2], the modulation equations derived in this work will prove to be
an effective tool to study several physically significant problems. The
soliton limit should prove to be particularly important in this respect,
similar to the KP equation [48, 49, 50].

We hope that the results of this work and the present discussion will provide
a stimulus for several further studies on these and related problems.

### Acknowledgments

We thank Alexandr Chernyavskiy and Dmitri Kireyev for many useful discussions
on topics related to this work. This research was partially supported by the
National Science Foundation under grant numbers DMS-1816934 and DMS-2009487.

## Appendix

###  A.1 Direct derivation of the periodic solutions of the NLS equation

Here we derive the periodic solutions of the NLS equation in an arbitrary
number of dimensions directly, without using the hydrodynamic system. We start
with the one-phase ansatz

$\psi(\mathbf{x},t)=\sqrt{\rho(z/\varepsilon)}\,{\mathrm{e}}^{i\Phi(z/\varepsilon)}\,,$
(A.1)

where, as before, the “fast variable” is $Z=\mathbf{k}\cdot\mathbf{x}-\omega
t$. Substituting (A.1) into (1.1) and separating into real and imaginary parts
yields respectively:

$\displaystyle(\sqrt{\rho})^{\prime\prime}-\sqrt{\rho}(\Phi^{\prime})^{2}+\frac{\omega}{\|\mathbf{k}\|^{2}}\sqrt{\rho}\,\Phi^{\prime}-\frac{2}{\|\mathbf{k}\|^{2}}\rho^{3/2}=0\,,$
(A.0a)
$\displaystyle\sqrt{\rho}\,\Phi^{\prime\prime}+\Big{(}2\Phi^{\prime}-\frac{\omega}{\|\mathbf{k}\|^{2}}\Big{)}(\sqrt{\rho})^{\prime}=0\,,$
(A.0b)

where, for brevity in this section, we denote $a=\|\mathbf{k}\|^{2}$.
Integrating (A.0b) yields $\Phi^{\prime}$ up to an integration constant $J$

$\Phi^{\prime}=\frac{J}{\|\mathbf{k}\|\,\rho}+\frac{\omega}{2a}\,.$ (A.1)

Substituting the phase relation (A.1), the real part (A.0a) reduces to:

$\displaystyle(\sqrt{\rho})^{\prime\prime}-\frac{J^{2}}{a\,\rho^{3/2}}+\Big{(}\frac{\omega}{2a}\Big{)}^{2}\sqrt{\rho}-\frac{2}{a}\rho^{3/2}=0.$
(A.2)

Multiplying by $2(\sqrt{\rho})^{\prime}$ and integrating with respect to $Z$
and letting $f=\rho$ yields:

$\displaystyle(f^{\prime})^{2}=\frac{4}{a}f^{3}-4\Big{(}\frac{\omega}{2a}\Big{)}^{2}f^{2}+4c_{1}f-\frac{4J^{2}}{a}\,.$
(A.3)

By substituting $f(Z)=A+By^{2}(Z)$, we get the following ODE for $y$:

$\displaystyle\displaystyle(y^{\prime})^{2}=\frac{1}{B^{2}}\Big{[}\frac{A^{3}}{a}-A^{2}\Big{(}\frac{\omega}{2a}\Big{)}^{2}+Ac_{1}-\frac{J^{2}}{a}\Big{]}\frac{1}{y^{2}}+\frac{1}{B}\Big{[}\frac{3A^{2}}{a}-2A\Big{(}\frac{\omega}{2a}\Big{)}^{2}+c_{1}\Big{]}+\Big{[}\frac{3A}{a}-\Big{(}\frac{\omega}{2a}\Big{)}^{2}\Big{]}y^{2}+\frac{B}{a}y^{4}.$
(A.4)

Now recall that the Jacobian elliptic sine $y(Z)=\mathop{\rm
sn}\nolimits(cZ|m)$ solves the ODE $(y^{\prime}/c)^{2}=(1-y^{2})(1-my^{2})$.
By requiring that (A.4) matches the ODE for the elliptic sine, one then
obtains (2.7), with $B=4m\|\mathbf{k}\|^{2}K_{m}^{2}$ as before, and with

$\displaystyle
J^{2}=4aK_{m}^{2}A\bigg{(}1+\frac{A}{4K_{m}^{2}a}\bigg{)}(A+4mK_{m}^{2}a)\,,$
(A.0a)
$\displaystyle\Big{(}\frac{\omega}{2a}\Big{)}^{2}=4K_{m}^{2}(1+m)+\frac{3A}{a}\,,\quad
c_{1}=\frac{1}{a}\,\big{[}\big{(}{4mK_{m}^{2}a}+A\big{)}(4K_{m}^{2}a+2A)+A\big{(}{4K_{m}^{2}a}+A\big{)}\big{]}\,.$
(A.0b)

Similar to section 2.2, we write the ODE (A.3) as $(f^{\prime})^{2}=P_{3}(f)$,
where

$P_{3}(f)=\frac{4}{a}\Big{[}f^{3}-a\Big{(}\frac{\omega}{2a}\Big{)}^{2}f^{2}+c_{1}\,a\,f-J^{2}\Big{]}=\frac{4}{a}(f-\lambda_{1})(f-\lambda_{2})(f-\lambda_{3})\,,$
(A.1)

with $\lambda_{1},\dots,\lambda_{3}$ given by (2.9). Note that the
requirements $a\geqslant 0$ and $0\leqslant m\leqslant 1$ again immediately
imply (2.11). The symmetric polynomials defined by $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$ are related to the above constants as

$e_{1}=\lambda_{1}+\lambda_{2}+\lambda_{3}={\omega^{2}}/{4a}\,,\quad
e_{2}=\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{3}\lambda_{1}=c_{1}\,a\,,\quad
e_{3}=\lambda_{1}\lambda_{2}\lambda_{3}=J^{2}.$ (A.2)

which also allow one to recover $A$, $a$, and $m$ when $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$ are known. The above solution contains $N+2$
independent parameters: $A$, $m$ and $\mathbf{k}$ [since $J$ and $\omega$ are
determined by (A.0a), (A.0b)]. Next we employ the Galilean invariance of the
NLS equation to apply a Galilean boost and thereby obtain the more general
family of solutions

$\tilde{\psi}(\mathbf{x},t)=\psi(\mathbf{x}-2\mathbf{v}t,t)\,{\mathrm{e}}^{i(\mathbf{v}\cdot\mathbf{x}-\|\mathbf{v}\|^{2}t)/\varepsilon}=\sqrt{\rho(\tilde{z}/\varepsilon)}\,{\mathrm{e}}^{i\tilde{\Phi}(\tilde{z}/\varepsilon,\mathbf{x},t)}\,,$
(A.0a)

where $\tilde{z}=\mathbf{k}\cdot\mathbf{x}-\tilde{\omega}t\,$, with
$\tilde{\omega}=\omega+2\mathbf{k}\cdot\mathbf{v}$, and where

$\tilde{\Phi}(\tilde{z}/\varepsilon,\mathbf{x},t)=\Phi(\tilde{z}/\varepsilon)+(\mathbf{v}\cdot\mathbf{x}-\|\mathbf{v}\|^{2}t)/\varepsilon\,.$
(A.0b)

The transformation adds the $N$ new independent parameters
$v_{1},\dots,v_{N}$. Therefore, the periodic solution of the NLS equation
(1.1) in $N$ spatial dimensions contains $2N+2$ independent real parameters:
$A$, $m$, $\mathbf{k}$ and $\mathbf{v}$, as expected.

###  A.2 Calculation of the solution amplitude and second frequency and
simplification of certain terms

Here we give a few additional details on the calculation of the fluid density.
Starting from (2.0b), using (2.2) and simplifying the resulting ODE, one has

$a\rho^{\prime\prime\prime}+a\rho^{\prime}\frac{(\rho^{\prime})^{2}-2\rho\rho^{\prime\prime}}{\rho^{2}}-4\rho\rho^{\prime}+\frac{4J^{2}\rho^{\prime}}{\rho^{2}}=0\,,$
(A.1)

where $a=\|\mathbf{k}\|^{2}$ for brevity. Integrating w.r.t. $Z$ yields

$a\rho^{\prime\prime}-a\frac{(\rho^{\prime})^{2}}{\rho}-2\rho^{2}+2c_{1}-\frac{4J^{2}}{\rho}=0\,,$
(A.2)

where $c_{1}$ is an arbitrary integration constant. Multiplying (A.2) by
$2\rho^{\prime}/\rho^{2}$ and integrating with respect to $Z$ again yields

$\displaystyle
a(\rho^{\prime})^{2}=4\rho^{3}-4c_{2}\rho^{2}+4c_{1}\rho-4J^{2}\,,$ (A.3)

with $c_{2}$ another arbitrary integration constant. Letting
$\rho(Z)=A+By^{2}(Z)$ yields the following ODE:

$\displaystyle(y^{\prime})^{2}=\frac{1}{B^{2}a}(A^{3}-A^{2}c_{2}+Ac_{1}-J^{2})\frac{1}{y^{2}}+\frac{1}{Ba}(3A^{2}-2Ac_{2}+c_{1})+\Big{(}\frac{3A}{a}-\frac{c_{2}}{a}\Big{)}y^{2}+\frac{B}{a}y^{4}.$
(A.4)

Now recall that the Jacobi elliptic sine $y(Z)=\mathop{\rm sn}\nolimits(cz|m)$
solves the ODE $(y^{\prime}/c)^{2}=(1-y^{2})(1-my^{2})$. By requiring that
(A.4) matches the ODE for the elliptic sine, one obtains (2.7), with the
coefficients as in (A.0a).

Next we obtain (2.13), which determines the frequency $\mu$ of the second
phase. As mentioned in section 2, to this end one can use the undifferentiated
version of (2.0b) [obtained from the real part of (1.1) using (2.1) and
(2.2)], which is

$-\omega\phi^{\prime}+2(\mathbf{k}\cdot\bar{\mathbf{u}})\,\phi^{\prime}+a(\phi^{\prime})^{2}+2\rho-\mu+\|\bar{\mathbf{u}}\|^{2}-\frac{a}{4}\Big{(}(\ln\rho)^{\prime\prime}+\frac{\rho^{\prime\prime}}{\rho}\Big{)}=0\,.$
(A.5)

Differentiating (A.5) w.r.t. $x$ and $y$ and collecting leading-order terms
yields (2.0b). However, (A.5) allows us to determine $\mu$ in a more
straightforward manner. Indeed, substituting (2.1) into equation (A.5) and
simplifying yields,

$2a\rho^{\prime\prime}-a\frac{(\rho^{\prime})^{2}}{\rho}-8\rho^{2}+C\rho-\frac{4J^{2}}{\rho}=0\,,$
(A.0a)

where

$C=4\mu-4\big{(}\|\bar{u}\|^{2}-(J\,\overline{\rho^{-1}}/g^{1/2})^{2}\big{)}$
(A.0b)

Multiplying (A.0a) by $\rho^{\prime}/\rho$ and integrating with respect to $Z$
yields

$\displaystyle a(\rho^{\prime})^{2}=4\rho^{3}-C\rho^{2}+4c_{3}\rho-4J^{2}\,,$
(A.1)

with an arbitrary integration constant $c_{3}$. Comparing the coefficients in
(A.3) and (A.1) we have $C=4c_{2}$ [as well as $c_{1}=c_{3}$], which, when
inserted in (A.0b), finally yields (2.13) for $\mu$.

Finally, we provide further details on how to simplify the modulation
equations (3.1) and in particular on how to obtain (3.0f). The averaged
conservation of momentum equation (3.0f), when written in terms of
$\lambda_{1},\dots,\lambda_{3}$ and $\mathbf{U}$, is

$\displaystyle(J\,\hat{\mathbf{k}}+\bar{\rho}\mathbf{U})_{t}+\nabla(\overline{\rho^{2}})+\nabla\cdot\bigg{[}2\bar{\rho}\mathbf{U}\otimes\mathbf{U}+2J\,(\hat{\mathbf{k}}\otimes\mathbf{U}+\mathbf{U}\otimes\hat{\mathbf{k}})+\bigg{(}\bigg{(}\overline{\frac{(\rho^{\prime})^{2}}{2\rho}}\bigg{)}+2\frac{J^{2}\overline{\rho^{-1}}}{\|\mathbf{k}\|^{2}}\bigg{)}\mathbf{k}\bigg{)}\otimes\mathbf{k}\bigg{]}=\mathbf{0}\,.$
(A.2)

Notice that averages containing $\rho_{z}$ can be evaluated by recalling that
$\rho(Z)$ satisfies the ODE (2.5). Differentiating (2.5) and using the
definition of the symmetric polynomials yields

$\frac{\|\mathbf{k}\|^{2}}{2}\rho^{\prime\prime}=3\rho^{2}-2e_{1}\rho+e_{2}\,,$
(A.0a)

and averaging over the fast variable $Z$ gives

$\overline{\rho^{2}}=\frac{2e_{1}\overline{\rho}-e_{2}}{3}\,.$ (A.0b)

Reordering the ODE (2.5) gives us

$\frac{\|\mathbf{k}\|^{2}\rho_{2}^{2}+4J^{2}}{\rho}=4\rho^{2}-4e_{1}\rho+4e_{2}\,.$
(A.0c)

Averaging again over the fast variable $Z$ and using (A.0b) yields

$\frac{\|\mathbf{k}\|^{2}}{4}\Big{(}\overline{\big{(}\frac{(\rho^{\prime})^{2}}{\rho}\big{)}}+\frac{4J^{2}}{\|\mathbf{k}\|^{2}}\overline{\rho^{-1}}\Big{)}=\frac{1}{3}{(2e_{2}-e_{1}\overline{\rho})}\,.$
(A.0d)

Finally, using (A.0d), equation (A.2) yields (3.0f).

###  A.3 Detailed steps in the derivation of the 3DNLS-Whitham system

We begin by expressing the modulation equations in terms of convective
derivatives. Using (7.1) one can see that

$D_{y}(q_{2})=q_{2,y}-q_{1}q_{2,x}=\frac{D_{z}(k_{1}q_{1})}{k_{1}}-q_{1}\frac{D_{z}k_{1}}{k_{1}}=D_{z}(q_{1})\,,$
(A.1)

which proves (7.0a). Moreover, using (3.0d) and the fact that
$\mathbf{p}=\bar{\mathbf{u}}_{\flat}-\bar{u}_{1}\mathbf{q}$, one has

$\displaystyle D_{y}(p_{2})$
$\displaystyle=(\bar{u_{3}}-q_{2}\bar{u_{1}})_{y}-q_{1}(\bar{u_{3}}-q_{2}\bar{u_{1}})_{x}=\bar{u}_{3,y}-q_{1}\bar{u}_{3,x}-q_{2}\bar{u}_{1,y}+q_{1}q_{2}\bar{u}_{1,x}-\bar{u}_{1}D_{y}(q_{2})\,,$
(A.0a)
$\displaystyle=\bar{u}_{2,z}-q_{1}\bar{u}_{1,z}-q_{2}\bar{u}_{2,x}+q_{1}q_{2}\bar{u}_{1,x}-\bar{u}_{1}D_{z}(q_{1})=D_{z}(p_{1})\,.$

which yields (7.0b).

Using the identity (7.0a) and straightforward algebra, we can rewrite (7.0b)
and (7.0c) as

$\displaystyle
D_{t}\mathbf{q}+2g\mathbf{D}_{\flat}U_{1}+2q_{1}(U_{1}\mathbf{D}_{\flat}q_{1}+\mathbf{D}_{\flat}p_{1})+2q_{2}(U_{1}\mathbf{D}_{\flat}q_{2}+\mathbf{D}_{\flat}p_{2})=0\,,$
(A.0a) $\displaystyle
D_{t}\mathbf{p}-2q_{1}U_{1}D_{y}\mathbf{p}-2q_{2}U_{1}D_{z}\mathbf{p}+\mathbf{D}_{\flat}(g(\tilde{e}_{1}-U_{1}^{2}))=0\,,$
(A.0b)

where $D_{t}$ is as in (7.1).

Next, we express the first conservation of waves equation in convective
derivative form. Recalling equations (7.10) and using (7.1) one can obtain the
following,

$\frac{D_{t}k_{1}}{k_{1}}+2(U_{1})_{x}+2\mathbf{q}\cdot(\mathbf{U}_{\flat})_{x}=0\,.$
(A.1)

Simplifying further we have (7.2), with $D_{x}$ as in (7.1) and with

$W_{1}=U_{1}(\|\mathbf{q}\|^{2})_{x}-2\mathbf{q}\cdot\mathbf{D}_{\flat}U_{1}+2\mathbf{q}\cdot\mathbf{p}_{x}\,.$
(A.2)

Moreover, using (A.0a) one can simplify $W_{1}$ further and obtain (7.0a).

Next, it can be easily seen that (7.3) becomes

$\displaystyle
D_{t}(U_{1}+J\,\overline{\rho^{-1}}/g^{1/2})+g(\tilde{e}_{1})_{x}+2gJ\,\overline{\rho^{-1}}/g^{1/2}(U_{1})_{x}+2U_{1}\|\mathbf{q}\|^{2}(J\,\overline{\rho^{-1}}/g^{1/2}+U_{1})_{x}$
$\displaystyle\kern
20.00003pt-2(U_{1}\mathbf{q}+\mathbf{p})\cdot{\kern-0.83337pt\lower
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\lower 0.20004pt\hbox{\kern
0.50003pt\hbox{$\nabla$}}\kern-8.33336pt\raise
0.20004pt\hbox{$\nabla$}\kern-8.33336pt\raise