2.1 Derivation of the MNLSE using Whitham’s method

The governing equation for irrotational inviscid two dimensional surface gravity waves is (see, for example, Phillips Reference Phillips1977)

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{\prime 2}\unicode[STIX]{x1D719}^{\prime }=0,\end{eqnarray}$$

together with boundary conditions

(2.2a,b ) $$\begin{eqnarray}\left.\unicode[STIX]{x1D719}_{t^{\prime }}^{\prime }+{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}^{\prime }\unicode[STIX]{x1D719}^{\prime })^{2}+gz^{\prime }=0\right|_{z^{\prime }=\unicode[STIX]{x1D702}^{\prime }};\quad \left.\unicode[STIX]{x1D702}_{t^{\prime }}^{\prime }+\unicode[STIX]{x1D719}_{x^{\prime }}^{\prime }\unicode[STIX]{x1D702}_{x^{\prime }}^{\prime }=\unicode[STIX]{x1D719}_{z^{\prime }}^{\prime }\right|_{z^{\prime }=\unicode[STIX]{x1D702}^{\prime }},\end{eqnarray}$$

along with the condition of no flow at the bottom $\unicode[STIX]{x1D719}_{z^{\prime }}^{\prime }\rightarrow 0$ as $z^{\prime }\rightarrow -\infty$ , where $\unicode[STIX]{x1D719}^{\prime }$ is the velocity potential, $\unicode[STIX]{x1D702}^{\prime }$ is the free surface displacement, $\unicode[STIX]{x1D735}^{\prime }=(\unicode[STIX]{x2202}_{x^{\prime }},\unicode[STIX]{x2202}_{z^{\prime }})$ and $g$ is the acceleration due to gravity. Although the governing equation is linear, the boundary conditions are nonlinear, and more severely, are evaluated at a dependent variable of the system, namely $\unicode[STIX]{x1D702}^{\prime }$ .

In order to make analytic progress, asymptotic approximations are employed. The wave slope $\unicode[STIX]{x1D716}\equiv a_{0}k_{0}$ is a measure of the nonlinearity of the system, and is assumed to be much less than 1. Furthermore, if $\unicode[STIX]{x1D6FF}k$ is a perturbation to the characteristic wavenumber $k_{0}$ , then we define the bandwidth as $\unicode[STIX]{x1D6E5}\equiv |\unicode[STIX]{x1D6FF}k|/k_{0}$ . Now, unlike the canonical approach where $\unicode[STIX]{x1D6E5}\sim \unicode[STIX]{x1D716}$ , we instead leave the exact relationship between the two parameters general (see also the discussion in Phillips (Reference Phillips1981)). Finally, the characteristic angular frequency, $\unicode[STIX]{x1D714}_{0}$ , is related to the wavenumber by the linear dispersion relation for deep water surface gravity waves, i.e.  $\unicode[STIX]{x1D714}_{0}^{2}=gk_{0}$ .

Following (Lo & Mei Reference Lo and Mei1985), we non-dimensionalize our system as follows:

(2.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D719}=(a_{0}^{2}\unicode[STIX]{x1D714}_{0})^{-1}\unicode[STIX]{x1D719}^{\prime },\quad \unicode[STIX]{x1D702}=a_{0}^{-1}\unicode[STIX]{x1D702}^{\prime },\\ x=k_{0}x^{\prime },\quad t=\unicode[STIX]{x1D714}_{0}t^{\prime },\quad z=k_{0}z^{\prime }.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Expanding the velocity potential and free surface displacement as a series, we have (Chu & Mei Reference Chu and Mei1970)

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=\overline{\unicode[STIX]{x1D719}}(X,T;\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})+\mathop{\sum }_{l=1}^{M}\mathop{\sum }_{m=1}^{M}\mathop{\sum }_{n\geqslant m}^{N}\unicode[STIX]{x1D716}^{n-1}\unicode[STIX]{x1D6E5}^{l-1}\unicode[STIX]{x1D719}_{lmn}(X,T,Z)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}\text{e}^{mz}+\text{c.c.}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\overline{\unicode[STIX]{x1D702}}(X,T;\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})+\mathop{\sum }_{l=1}^{M}\mathop{\sum }_{m=1}^{M}\mathop{\sum }_{n\geqslant m}^{N}\unicode[STIX]{x1D716}^{n-1}\unicode[STIX]{x1D6E5}^{l-1}\unicode[STIX]{x1D702}_{lmn}(X,T)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}+\text{c.c.}, & \displaystyle\end{eqnarray}$$

where $X=\unicode[STIX]{x0394}x,Z=\unicode[STIX]{x0394}z,T=\unicode[STIX]{x0394}t,\unicode[STIX]{x1D703}=x-t$ and $\text{c.c.}$ means complex conjugate; $\overline{\unicode[STIX]{x1D719}}$ and $\overline{\unicode[STIX]{x1D702}}$ are the mean (phase-averaged) velocity potential and free surface displacement, respectively, with their scale to be determined. Finally, we take the complex valued coefficient $\unicode[STIX]{x1D719}_{111}\equiv B$ as our dependent variable. This has important implications for the structure of the resulting governing equation (see the discussion in Kit & Shemer (Reference Kit and Shemer2002) and Trulsen (Reference Trulsen2006)).

The expansions of the dependent variables to third order in $\unicode[STIX]{x1D716}$ and second order in $\unicode[STIX]{x1D6E5}$ , appropriate for the MNLSE, are (Trulsen & Dysthe Reference Trulsen and Dysthe1997; Trulsen Reference Trulsen2006; Pizzo & Melville Reference Pizzo and Melville2016)

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D716}\overline{\unicode[STIX]{x1D719}}+\frac{1}{2}\left(\left[B-\text{i}\unicode[STIX]{x0394}zB_{X}-\frac{\unicode[STIX]{x1D6E5}^{2}}{2}z^{2}B_{XX}\right]\text{e}^{z}\text{e}^{\text{i}\unicode[STIX]{x1D703}}+\text{c.c.}\right),\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702} & = & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x1D6E5}\overline{\unicode[STIX]{x1D702}}+\frac{1}{2}\left(\left[\text{i}B+\frac{\unicode[STIX]{x1D6E5}}{2}B_{X}+\frac{\text{i}\unicode[STIX]{x1D6E5}^{2}}{8}B_{XX}+\frac{\text{i}\unicode[STIX]{x1D716}^{2}}{8}|B|^{2}B\right]\text{e}^{\text{i}\unicode[STIX]{x1D703}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\left[-\frac{\unicode[STIX]{x1D716}}{2}B^{2}+\text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x0394}BB_{X}\right]\text{e}^{2\text{i}\unicode[STIX]{x1D703}}+\frac{3\text{i}\unicode[STIX]{x1D716}^{2}}{8}B^{3}\text{e}^{3\text{i}\unicode[STIX]{x1D703}}+\text{c.c.}\right).\end{eqnarray}$$

Recall, equations (2.1), (2.2) can be reformulated in terms of a variational principle through the action (Luke Reference Luke1967; Zakharov Reference Zakharov1968; Miles Reference Miles1977)

(2.8) $$\begin{eqnarray}\int L\,\text{d}x\,\text{d}t=\int \left\{\unicode[STIX]{x1D716}\unicode[STIX]{x1D713}\unicode[STIX]{x1D702}_{t}-\unicode[STIX]{x1D716}\left(\int _{-\infty }^{\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}}\frac{1}{2}(\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D719})^{2}\,\text{d}z+\frac{1}{2}\unicode[STIX]{x1D702}^{2}\right)\right\}\,\text{d}x\,\text{d}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D719}(x,z=\unicode[STIX]{x1D716}\unicode[STIX]{x1D702},t)$ is the velocity potential evaluated at the free surface.

Next, we substitute the form of the velocity potential and free surface displacement into (2.8). Furthermore, following Lo & Mei (Reference Lo and Mei1985), we map into a spatial reference frame by introducing the following transformation

(2.9a,b ) $$\begin{eqnarray}(2X-T)=\unicode[STIX]{x1D70F};\quad \unicode[STIX]{x0394}X=\unicode[STIX]{x1D712}.\end{eqnarray}$$

This reference frame is travelling at the linear group velocity, and for a fixed position $X$ , $\unicode[STIX]{x1D70F}$ is proportional to the negative of the elapsed time, while $\unicode[STIX]{x1D712}$ is related to the fetch (or stretched distance). Finally, for clarity of presentation we transform from the dependent variable $B$ to the dependent variable in the mapped frame  $A$ .

Employing the mapping given by (2.9), the averaged Lagrangian (Whitham Reference Whitham1965, Reference Whitham1974), then becomes

(2.10) $$\begin{eqnarray}{\mathcal{L}}\equiv \frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}L\,\text{d}\unicode[STIX]{x1D703}=\frac{\text{i}}{2}\left(A\frac{\unicode[STIX]{x2202}A^{\ast }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}-A^{\ast }\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}\right)-{\mathcal{H}},\end{eqnarray}$$

where ${\mathcal{H}}$ is defined as (cf. Gramstad & Trulsen Reference Gramstad and Trulsen2011, appendix C)

(2.11) $$\begin{eqnarray}{\mathcal{H}}=|A_{\unicode[STIX]{x1D70F}}|^{2}-\frac{\unicode[STIX]{x1D6FE}_{0}}{2}|A|^{4}-\frac{\unicode[STIX]{x1D6FC}_{0}}{2}|A|^{2}\mathbb{H}(|A|_{\unicode[STIX]{x1D70F}}^{2})+\text{i}\frac{\unicode[STIX]{x1D6FD}_{0}}{4}|A|^{2}(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast }),\end{eqnarray}$$

and $\unicode[STIX]{x1D6FE}_{0}=\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}^{2}$ , $\unicode[STIX]{x1D6FC}_{0}=2\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ , $\unicode[STIX]{x1D6FD}_{0}=8\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ , while $^{\ast }$ denotes the complex conjugate. Note, $\unicode[STIX]{x1D6FE}_{0}$ is a measure of the ratio of the $O(\unicode[STIX]{x1D716}^{2})$ nonlinearity to dispersion (this term has been shown to be important for statistical analogues of the NLSE (Alber Reference Alber1978; Onorato et al. Reference Onorato, Osborne, Serio and Bertone2001; Janssen Reference Janssen2003) as well as the deterministic dynamics of the NLSE (Anderson Reference Anderson1983)) and the Hilbert transform $\mathbb{H}$ is defined as

(2.12) $$\begin{eqnarray}\mathbb{H}(|A|_{\unicode[STIX]{x1D70F}}^{2})=\frac{1}{\unicode[STIX]{x03C0}}P.V.\int _{-\infty }^{\infty }\frac{\unicode[STIX]{x2202}|A(\unicode[STIX]{x1D712},\unicode[STIX]{x1D70F}^{\prime })|^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}^{\prime }}\frac{\text{d}\unicode[STIX]{x1D70F}^{\prime }}{\unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}^{\prime }},\end{eqnarray}$$

with $P.V.$ meaning we are to take the principal value of the integral (Titchmarsh Reference Titchmarsh1948). Following Janssen (Reference Janssen1983, see also Akylas Reference Akylas1989) we have rewritten our action in terms of one variable by recognizing that $\overline{\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D70F}}}$ and $\overline{\unicode[STIX]{x1D711}_{Z}}$ are harmonic conjugates, where $\unicode[STIX]{x1D711}$ is the velocity potential in this spatial reference frame. That is, these variables are real and imaginary parts of an analytic function in the lower half plane and hence we can connect the two variables, evaluated at $Z=0$ , via the Hilbert transform $\mathbb{H}$ .

As the linear dispersion is exactly represented in this spatial reference frame (Trulsen et al. Reference Trulsen, Kliakhandler, Dysthe and Velarde2000; Kit & Shemer Reference Kit and Shemer2002), only one term is needed to completely describe it. Besides making this equation more physically accurate, this also offers numerical advantages (Lo & Mei Reference Lo and Mei1985), and simplifies the equations algebraically.

The action associated with the Lagrangian density in the spatial reference frame, i.e. equation (2.10), is defined as

(2.13) $$\begin{eqnarray}{\mathcal{S}}(A,A^{\ast })=\int _{\unicode[STIX]{x1D712}_{0}}^{\unicode[STIX]{x1D712}_{f}}\int _{\mathbb{R}}{\mathcal{L}}\,\text{d}\unicode[STIX]{x1D70F}\,\text{d}\unicode[STIX]{x1D712},\end{eqnarray}$$

which is a functional over all admissible functions satisfying the prescribed conditions at $A(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D712}_{0})$ and $A(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D712}_{f})$ . Hamilton’s principle states that the governing equations of the system are found by requiring that ${\mathcal{S}}$ be stationary. That is, we seek solutions $(A,A^{\ast })$ such that

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}{\mathcal{S}}={\mathcal{S}}(A+\unicode[STIX]{x1D6FF}A,A^{\ast }+\unicode[STIX]{x1D6FF}A^{\ast })-{\mathcal{S}}(A,A^{\ast })=0,\end{eqnarray}$$

for infinitesimal $(\unicode[STIX]{x1D6FF}A,\unicode[STIX]{x1D6FF}A^{\ast })$ . Substituting the Lagrangian density into the action, then applying Hamilton’s principle, and recalling the anti-self-adjoint nature of the Hilbert transform (Titchmarsh Reference Titchmarsh1948), gives us (Lo & Mei Reference Lo and Mei1985; Kit & Shemer Reference Kit and Shemer2002)

(2.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}+\text{i}\frac{\unicode[STIX]{x2202}^{2}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}^{2}}+\text{i}\unicode[STIX]{x1D6FE}_{0}|A|^{2}A+\text{i}\unicode[STIX]{x1D6FC}_{0}A\mathbb{H}\left(\frac{\unicode[STIX]{x2202}|A|^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\right)+\unicode[STIX]{x1D6FD}_{0}|A|^{2}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=0.\end{eqnarray}$$

The constants $(\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FD}_{0})$ also make it easy to track the induced mean flow (which is related to the term with coefficient $\unicode[STIX]{x1D6FC}_{0}$ ), and the asymmetric self-steepening term (which is related to the term with coefficient $\unicode[STIX]{x1D6FD}_{0}$ ) in the ensuing calculations. Note, for $\unicode[STIX]{x1D6E5}=\unicode[STIX]{x1D716}$ , we return the standard form of the MNLSE.

This equation is usually written as a set of coupled partial differential equations (PDEs) in two dependent variables, $(A,\overline{\unicode[STIX]{x1D711}}|_{Z=0})$ . The two variables are connected by the relation

(2.16) $$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}Z}\right|_{Z=0}=\frac{\unicode[STIX]{x2202}|A|^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}},\end{eqnarray}$$

which implies that gradients in the radiation stress lead to the generation of a mean flow (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1964; Dysthe Reference Dysthe1979; McIntyre Reference McIntyre1981).

2.2 Conservation laws

A benefit of the above use of Whitham’s method is that symmetries of the action (2.13) may be connected with conservation laws using Noether’s theorem. In particular, Pizzo & Melville (Reference Pizzo and Melville2016) were able to identify several symmetries of the action of the Lagrangian associated with (2.15). They found that (2.15) conserves wave action,

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}|A|^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\left(\text{i}(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast })+\frac{\unicode[STIX]{x1D6FD}_{0}}{2}|A|^{4}\right)=0,\end{eqnarray}$$

as well as the wave action weighted frequency (Trulsen & Dysthe Reference Trulsen and Dysthe1997)

(2.18) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}(\text{i}(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast }))+\frac{\unicode[STIX]{x2202}{\mathcal{G}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+2\unicode[STIX]{x1D6FC}_{0}|A|_{\unicode[STIX]{x1D70F}}^{2}\mathbb{H}(|A|_{\unicode[STIX]{x1D70F}}^{2})=0,\end{eqnarray}$$

where

(2.19) $$\begin{eqnarray}{\mathcal{G}}=2|A_{\unicode[STIX]{x1D70F}}|^{2}-(AA_{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70F}}^{\ast }-A^{\ast }A_{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70F}})-\unicode[STIX]{x1D6FE}_{0}|A|^{4}+\text{i}\unicode[STIX]{x1D6FD}_{0}|A|^{2}(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast })-2\unicode[STIX]{x1D6FC}_{0}|A|^{2}\mathbb{H}(|A|_{\unicode[STIX]{x1D70F}}^{2}).\end{eqnarray}$$

As was noted in Pizzo & Melville (Reference Pizzo and Melville2016), equation (2.18) does not take the form of a local conservation law due to the last term in (2.18), which cannot be written as a perfect derivative. However, the integral quantity $\int \text{i}(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast })\,\text{d}\unicode[STIX]{x1D70F}$ is still conserved, as the final term on the left-hand side of (2.18) integrates to zero due to the anti-self-adjoint nature of the Hilbert transform (Titchmarsh Reference Titchmarsh1948).

3.1 Variational equations

Following Pizzo & Melville (Reference Pizzo and Melville2016), we use the variational method (Anderson Reference Anderson1983; Agrawal Reference Agrawal2007; Chandrasekhar Reference Chandrasekhar2013), to derive approximate analytic forms for the evolution of physically important packet parameters. This qualifies and quantifies the packet evolution, and forms a useful tool to understand the bulk scale properties of the wave envelope during focusing. The success of this approach relies on a trial function which captures the essence of the phenomenon in question, without being too algebraically complicated.

To this end, we recall that in the absence of linear dispersion, as was discussed by Pizzo & Melville (Reference Pizzo and Melville2016), equations (2.17) and (2.18) decouple, and may be solved using the method of characteristics. For small $\unicode[STIX]{x1D6FD}_{0}\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D712}$ , solutions to the wave action equation (i.e. (2.17)) takes the form

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D70C}\approx \unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D6FD}_{0}\unicode[STIX]{x1D70C}_{0}^{2}\unicode[STIX]{x1D70C}_{0}^{\prime }\unicode[STIX]{x1D712},\end{eqnarray}$$

where the subscript $0$ means we are evaluating these functions at $\unicode[STIX]{x1D712}=0$ .

This implies that that for an initially chirped Gaussian wave packet, we should try a trial function of the form (Anderson Reference Anderson1983)

(3.3) $$\begin{eqnarray}{\mathcal{A}}=(a(\unicode[STIX]{x1D712})\text{e}^{-\unicode[STIX]{x1D70F}^{2}/2w(\unicode[STIX]{x1D712})^{2}}+\unicode[STIX]{x1D6FD}_{0}s(\unicode[STIX]{x1D712})\unicode[STIX]{x1D70F}\text{e}^{-3\unicode[STIX]{x1D70F}^{2}/2w(\unicode[STIX]{x1D712})^{2}})\text{e}^{\text{i}(C(\unicode[STIX]{x1D712})\unicode[STIX]{x1D70F}^{2}/2+G(\unicode[STIX]{x1D712})+\unicode[STIX]{x1D6FD}_{0}F(\unicode[STIX]{x1D712})\unicode[STIX]{x1D70F})},\end{eqnarray}$$

where we are solving for the six real functions $(a,s,w,C,F,G)$ , governing the amplitude of the symmetric component of the packet $a$ , the asymmetric amplitude $s$ , the packet width $w$ , the chirp $C$ , the mean frequency $F$ and the phase shift $G$ .

To employ the variational technique we then substitute our trial wave function into the action, i.e. (2.13), so that to the same order of approximation as (2.15), we find

(3.4) $$\begin{eqnarray}\displaystyle \mathfrak{L} & = & \displaystyle \int _{-\infty }^{\infty }L\,\text{d}\unicode[STIX]{x1D70F}=\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\unicode[STIX]{x1D6FC}_{0}\unicode[STIX]{x1D718}_{0}a^{4}-\sqrt{2}\frac{a^{2}}{w}+\unicode[STIX]{x1D6FE}_{0}a^{4}w\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & & \displaystyle -\,\sqrt{2}a^{2}C^{2}w^{3}+\frac{a^{2}w^{3}C^{\prime }}{\sqrt{2}}+2\sqrt{2}a^{2}wG^{\prime },\end{eqnarray}$$

where $\unicode[STIX]{x1D718}_{0}$ is a constant given by

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D718}_{0}=\int _{-\infty }^{\infty }\text{e}^{-(\unicode[STIX]{x1D70F}/w)^{2}}\mathbb{H}((\text{e}^{-(\unicode[STIX]{x1D70F}/w)^{2}})_{\unicode[STIX]{x1D70F}})\,\text{d}\unicode[STIX]{x1D70F}=1,\end{eqnarray}$$

as the Hilbert transform of the Gaussian function may be written in terms of the Dawson function (Temme Reference Temme2010), so that this integral may be computed analytically.

We begin by noting that the integral form of conservation of wave action, i.e. equation (2.17), becomes

(3.7) $$\begin{eqnarray}E_{0}=\int _{-\infty }^{\infty }|A|^{2}\,\text{d}\unicode[STIX]{x1D70F}=\sqrt{\unicode[STIX]{x03C0}}a^{2}w+O(\unicode[STIX]{x1D6FD}_{0}^{2}),\end{eqnarray}$$

where $E_{0}$ is a specified constant (independent of the scales $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6E5}$ ) based on the initial wave form. Furthermore, conservation of mean frequency, i.e. the integral of (2.18), gives

(3.8) $$\begin{eqnarray}P_{0}=\frac{\text{i}}{E_{0}}\int _{-\infty }^{\infty }(A^{\ast }A_{\unicode[STIX]{x1D70F}}-AA_{\unicode[STIX]{x1D70F}}^{\ast })\,\text{d}\unicode[STIX]{x1D70F}=\frac{aw\unicode[STIX]{x1D6FD}_{0}}{8}\left(aF+\frac{Csw^{2}}{2\sqrt{2}}\right),\end{eqnarray}$$

which for the initially chirped Gaussian function (3.1) must be taken to be zero, so that

(3.9) $$\begin{eqnarray}F=-\frac{1}{2\sqrt{2}}\frac{sw^{5/2}C}{E_{0}^{1/2}}.\end{eqnarray}$$

We can now apply Hamilton’s principle to each of the dependent variables, yielding a set of coupled ordinary differential equations. Taking variations with respect to $w$ implies

(3.10) $$\begin{eqnarray}G^{\prime }+\frac{3}{4}w^{2}C^{\prime }-\frac{3}{2}C^{2}w^{2}+\frac{1}{2w^{2}}+\frac{\unicode[STIX]{x1D6FE}_{0}a^{2}}{2\sqrt{2}}=0,\end{eqnarray}$$

while variations with respect to $G$ gives

(3.11) $$\begin{eqnarray}(a^{2}w)^{\prime }=0,\end{eqnarray}$$

which is consistent with equation (3.7).

Next, varying $\mathfrak{L}$ with respect to $a$ we find

(3.12) $$\begin{eqnarray}G^{\prime }+\frac{1}{4}w^{2}C^{\prime }-\frac{1}{2}C^{2}w^{2}+\frac{\unicode[STIX]{x1D6FC}_{0}a^{2}}{\sqrt{\unicode[STIX]{x03C0}}w}-\frac{1}{2w^{2}}+\unicode[STIX]{x1D6FE}_{0}\frac{a^{2}}{\sqrt{2}}=0,\end{eqnarray}$$

while variations with respect to $C$ yields

(3.13) $$\begin{eqnarray}4a^{2}w^{3}C+(a^{2}w^{3})^{\prime }=0.\end{eqnarray}$$

Next, recall Pizzo & Melville (Reference Pizzo and Melville2016) found that the centroid of packets governed by the MNLSE, in the frame of reference moving at the linear group velocity, have a finite amplitude correction analogous to the Stokes correction for permanent progressive waves. There, it was found that (for $P_{0}=0$ )

(3.14) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D712}}\frac{1}{E_{0}}\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}|A|^{2}\,\text{d}\unicode[STIX]{x1D70F}=\frac{\unicode[STIX]{x1D6FD}_{0}}{2E_{0}}\int _{-\infty }^{\infty }|A|^{4}\,\text{d}\unicode[STIX]{x1D70F}.\end{eqnarray}$$

Substituting our trial function into this equation, we find a relationship for $s$ ,

(3.15) $$\begin{eqnarray}s^{\prime }+\frac{5}{2}\frac{w^{\prime }}{w}s=\frac{8E_{0}^{3/2}}{w^{7/2}}.\end{eqnarray}$$

Next, from equation (3.7), equation (3.13) can be rewritten to find

(3.16) $$\begin{eqnarray}C=-\frac{1}{2}\frac{w^{\prime }}{w}=-\frac{1}{2}\frac{\text{d}\ln w}{\text{d}\unicode[STIX]{x1D712}}.\end{eqnarray}$$

Then, subtracting (3.10) from (3.12) and substituting in (3.16), we have

(3.17) $$\begin{eqnarray}w^{\prime \prime }+\unicode[STIX]{x1D6FE}_{0}E_{0}\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{1}{w^{2}}+4\left(-1+\frac{\unicode[STIX]{x1D6FC}_{0}E_{0}}{\unicode[STIX]{x03C0}}\right)\frac{1}{w^{3}}=0.\end{eqnarray}$$

Multiplying by $w^{\prime }$ and integrating once, we find an equation describing a particle trapped in a potential well,

(3.18) $$\begin{eqnarray}\frac{1}{2}\left(\frac{\text{d}w}{\text{d}\unicode[STIX]{x1D712}}\right)^{2}+\unicode[STIX]{x1D6F1}(w)=0,\end{eqnarray}$$

where the potential $\unicode[STIX]{x1D6F1}(w)$ is given by

(3.19) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}(w)=2\left(1-\frac{\unicode[STIX]{x1D6FC}_{0}E_{0}}{\unicode[STIX]{x03C0}}\right)\frac{1}{w^{2}}-\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D6FE}_{0}E_{0}}{w}+\unicode[STIX]{x1D6F1}_{0},\end{eqnarray}$$

and $\unicode[STIX]{x1D6F1}_{0}$ is an integration constant chosen to satisfy the initial conditions, that is

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{0}=-1/2w_{0}^{^{\prime }2}-2\left(1-\frac{\unicode[STIX]{x1D6FC}_{0}E_{0}}{\unicode[STIX]{x03C0}}\right)\frac{1}{w_{0}^{2}}+\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D6FE}_{0}E_{0}}{w_{0}}.\end{eqnarray}$$

The initial velocity of the particle in the well is set by the chirp of the packet and the initial width, as is evident by (3.16). Equation (3.19) makes it clear that the dynamics of $w$ is set by the size of $(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6E5})^{2}$ and $\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ . Note, throughout the rest of this discussion we take $E_{0}$ to be $\sqrt{\unicode[STIX]{x03C0}}$ , so that $a_{0}^{2}w_{0}=1$ .

Examples of potentials for two values of $\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ are shown in figure 2. There, the black dot represents the initial packet width (taken to be $w_{0}=1$ ). We take $C_{0}<0$ so that the packet width $w$ decreases, leading to focusing. For increasing values of $\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ , this focusing event becomes more pronounced as linear dispersion is retarded by the induced mean flow.

The packet width may be solved for implicitly, as we may rewrite (3.18) as

(3.21) $$\begin{eqnarray}\sqrt{2}\unicode[STIX]{x1D712}=\int _{1}^{w}\frac{\tilde{w}\text{d}\tilde{w}}{\sqrt{-\tilde{w}^{2}\unicode[STIX]{x1D6F1}_{0}+\sqrt{2}\unicode[STIX]{x1D6FE}_{0}\tilde{w}-2(1-\unicode[STIX]{x1D6FC}_{0}/\sqrt{\unicode[STIX]{x03C0}})}}.\end{eqnarray}$$

Figure 2.

The potential well given by (3.19) for $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D6E5}=0.1$ (black line) and $\unicode[STIX]{x1D716}=4\unicode[STIX]{x1D6E5}/3=0.1$ (red line). The dashed lines show the corresponding potentials for the nonlinear Schrödinger equation. Note, the red line has $\unicode[STIX]{x1D6F1}^{\prime }(w_{0}=1)<0$ , implying the nonlinearity is stronger than the dispersion initially, aiding the chirp in compressing the packet, i.e. leading to smaller values of $w$ .

Next, equation (3.15) may be integrated to find

(3.22) $$\begin{eqnarray}s=\frac{8E_{0}^{3/2}}{w^{5/2}}\int \frac{\text{d}\tilde{\unicode[STIX]{x1D712}}}{w(\tilde{\unicode[STIX]{x1D712}})}.\end{eqnarray}$$

The solutions to these equations are discussed in more detail below.

3.3 Numerical simulations of the MNLSE

To better understand these predictions, we now look at several numerical examples. The numerical scheme is a split-step pseudo-spectral method, where the step size is chosen so that the energy is conserved to at least 1 part in $10^{6}$ and the mean frequency to 1 part in $10^{4}$ . The number of modes was increased until convergence was found. In particular, we use $2^{12}$ modes and step sizes of the order of $10^{-5}$ .

Note, one of the difficulties in applying ideas from weakly nonlinear narrow banded theory to areas outside of their strict asymptotic region of validity is the choice in defining the (small) parameters $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6E5}$ (see the discussion in § 2 of Perlin et al. (Reference Perlin, Choi and Tian2013)). For example, does one define the parameter $\unicode[STIX]{x1D716}$ as the initial slope, or the largest slope, as predicted by some theory? A technique used in laboratory experiments is to characterize the packet slope $S$ based on the linear prediction of the maximum slope at focusing. This has been shown to be effective in describing the strength of breaking waves for laboratory wave packets (Drazen et al. Reference Drazen, Melville and Lenain2008), which, due to wave channel length limitations, must be broadband in order to start as linear waves and progress to breaking over a reasonable fetch (usually on the order of $10$  m). Based on this, we connect the parameters $(\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})$ with the predictions of the maximum modulations to amplitude and spatial/temporal scales according to linear theory.

Figure 3 shows a comparison of the predictions of the trial functions with the MNLSE. Here, $\unicode[STIX]{x1D716}=\unicode[STIX]{x1D6E5}=0.1$ . We can see good agreement between the predictions of the amplitude $|A|$ and the trial function amplitude ${\mathcal{A}}=|{\mathcal{A}}|\text{e}^{\text{i}\unicode[STIX]{x1D717}}$ , while the frequency evolution $\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}$ (where we recall $A=|A|\text{e}^{\text{i}\unicode[STIX]{x1D703}}$ ) agrees with trial function predictions until focusing, when the $\unicode[STIX]{x1D70F}$ structure is more complicated than the low-order polynomials we assumed for $\unicode[STIX]{x1D717}$ . The lowest row shows the predictions of the magnitude of the envelope slope $|A_{\unicode[STIX]{x1D70F}}|$ of the envelope during focusing, and there is relatively good agreement for small values of $\unicode[STIX]{x1D712}$ . However, this simple model is unable to predict the slope at focusing, as the asymmetric term is not part of the energy balance to this order and grows secularly.

3.4 Maximum slope at focusing

In this section we seek to model the maximum slope at focusing, based on the initial packet parameters $(\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})$ . The variational method hinted at the importance of $\unicode[STIX]{x1D6FE}_{0},\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FD}_{0}$ in describing the focusing event. In particular, the values of $w_{\ast }$ , which we recall is the minimum packet width predicted by the variational method, sets the scale of the envelope at focusing. These values are not uniquely defined as a function of $\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5}$ , and to first order in $\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}$ , it may be seen that when $\unicode[STIX]{x1D6E5}=d_{0}\unicode[STIX]{x1D716}+d_{1}\unicode[STIX]{x1D716}^{2}$ , then lines of constant $w_{\ast }$ may be found for $d_{0}^{2}/2=w_{\ast }(1-w_{\ast })/(4w_{\ast }^{2}-2)$ and $d_{1}$ a more complicated function of $w_{\ast }$ which is not presented here as it is not relevant for the ensuing discussion. This hints at a quasi self-similarity in the focusing event.

Next, the free surface $\unicode[STIX]{x1D702}$ is given by (2.7). To better understand the dependence of max $(|\unicode[STIX]{x1D702}_{x}|)$ on $(\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})$ , the MNLSE was numerically integrated with the initial conditions given by (3.1) for 400 values of $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6E5}$ , and the maximum of the modulus of the slope is shown in figure 4. Note, the three decibel width of the initial spectrum is chosen as a measure of the packet bandwidth, $\unicode[STIX]{x1D6E5}$ . For comparison, the same simulations are conducted for the NLSE and are shown. The slope predictions are larger for the MNLSE, and this becomes more pronounced with increasing $\unicode[STIX]{x1D716}$ and decreasing $\unicode[STIX]{x1D6E5}$ (which corresponds to increasing $\unicode[STIX]{x1D6FE}_{0},\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FD}_{0}$ ).

Figure 4.

The maximum of the modulus of the slope, max $(|\unicode[STIX]{x1D702}_{x}|)$ , according to the modified nonlinear Schrodinger equation (upper surface) and nonlinear Schrodinger equation (lower surface) as a function of $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6E5}$ . Here, 400 runs were conducted, and the displayed surfaces are interpolated based on these data points. Note, the MNLSE predicts noticeably larger slopes than its lower-order counterpart, in accordance with the theoretical considerations of § 3. Slopes are largest for large $\unicode[STIX]{x1D716}$ and small $\unicode[STIX]{x1D6E5}$ , corresponding to larger values of $\unicode[STIX]{x1D6FE}_{0},\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FD}_{0}$ .

We define the maximum of the modulus of the free surface slope as ${\mathcal{T}}$ , that is

(3.27) $$\begin{eqnarray}{\mathcal{T}}\equiv \text{max}(|\unicode[STIX]{x1D702}_{x}|).\end{eqnarray}$$

For linear waves, we have ${\mathcal{T}}=\unicode[STIX]{x1D716}$ , by definition. From (3.26), the definition of $\unicode[STIX]{x1D702}$ (i.e. (2.7)), and the asymptotic self-similarity of $w_{\ast }$ described above, we propose that for the MNLSE

(3.28) $$\begin{eqnarray}{\mathcal{T}}=\unicode[STIX]{x1D716}(1+f(\unicode[STIX]{x1D705})),\end{eqnarray}$$

with $\unicode[STIX]{x1D705}$ as defined in (3.25) and $f$ a function, to be constrained by the numerical simulations.

Next, by inspection of figure 4 we expect the growth of the slope to be sensitive to values of $\unicode[STIX]{x1D6FE}_{0}$ and $\unicode[STIX]{x1D6FC}_{0}$ . To this end, we assume that the function $f$ is exponential, and find the following fit

(3.29) $$\begin{eqnarray}{\mathcal{T}}=\unicode[STIX]{x1D716}\left(h_{0}+h_{1}\text{e}^{h_{2}\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D6E5}^{2}(1+h_{3}\unicode[STIX]{x1D6E5})}\right),\end{eqnarray}$$

with $(h_{0},h_{1},h_{2},h_{3})=(-0.37,1.18,0.20,10.05)$ given by a least squares fit with the data. A comparison of this model with the data is shown in figure 5 in the form of a contour plot. There, bulk scale agreement between this simple model and the data is shown. Note, these coefficients also have a dependence on the chirp of the packet  $C$ .

Finally, note that (3.29) implies

(3.30) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}=\frac{h_{2}h_{3}\unicode[STIX]{x1D716}^{2}+\sqrt{h_{2}}\unicode[STIX]{x1D716}\sqrt{4\ln ({\mathcal{T}}/\unicode[STIX]{x1D716}-h_{0})+h_{2}h_{3}^{2}\unicode[STIX]{x1D716}^{2}-4\ln (h_{1})}}{2\ln ({\mathcal{T}}/\unicode[STIX]{x1D716}-h_{0})-2\ln (h_{1})}.\end{eqnarray}$$

Figure 5.

A comparison of the contours of max $(|\unicode[STIX]{x1D702}_{x}|)$ as a function of $(\unicode[STIX]{x1D716},\unicode[STIX]{x1D6E5})$ , with the model ${\mathcal{T}}$ (dashed lines) given by (3.29). There is a bulk scale agreement between this relatively simple model and the numerical simulations.