#  $\bm{SU(5)}$ grand unification and $W$-boson mass

Goran Senjanović1,2 and Michael Zantedeschi1,3 1 Arnold Sommerfeld Center,
Ludwig-Maximilians University, Munich, Germany  2 International Centre for
Theoretical Physics, Trieste, Italy  3 Max-Planck-Institute for Physics,
Munich, Germany

###### Abstract

A realistic extension of the minimal $SU(5)$ theory consisting of the addition
of an adjoint fermion is known to predict light real fermion and scalar weak
triplets, potentially accessible at the LHC. These particles, in addition to
playing a key role in gauge coupling unification, have profound
phenomenological implications. The fermion triplet, that through the seesaw
mechanism offers a testable origin of neutrino mass, has been already
extensively discussed. The scalar triplet develops a vacuum expectation value
that modifies the $W$-boson mass. We show that its low-energy effective theory
is remarkably predictive: in the leading approximation, all the relevant
physical processes involving the scalar triplet depend only on its mass and
the deviation from the Standard Model $W$-mass value.

Introduction Some years ago it was shown that a simple realistic extension
Bajc and Senjanović (2007) of the minimal $SU(5)$ theory Georgi and Glashow
(1974), based on small representations, predicts a fermion and a scalar weak
triplets at today’s energies. These multiplets are necessary remnants from the
high-energy theory. This goes against the long-held prejudice of a desert at
energies between the electroweak ($M_{W}$) and grand unifying ($M_{G}$)
scales. Simply, the structural consistency of the model requires the existence
of an oasis close to $M_{W}$, as we recapitulate below.

Notice that we are not just saying that there could be such an oasis - after
all, in baroque models it is trivially possible - rather we wish to call
attention to the fact that such light states are a must in this strongly
motivated simple model.

The fermion triplet plays a central role in providing neutrino masses through
the seesaw mechanism, and thus allows to test the origin of neutrino mass at
the LHC or the next hadron collider. Its physical consequences have been amply
discussed Bajc  _et al._ (2007); Arhrib  _et al._ (2010). At first glance, the
physics of the scalar triplet seems less exciting. However, due to the mixing
with the Standard Model (SM) Higgs boson, it develops a small vacuum
expectation value (VEV) which modifies the $W$-boson mass. This fits nicely
with the recent CDF experiment claim Aaltonen  _et al._ (2022) of $W$-mass
deviation from the SM value.

Before we proceed, a comment is in order. The scalar triplet was also
introduced by hand at low energies Ross and Veltman (1975), and its
phenomenology has been studied over the years, see e.g., Gunion  _et al._
(1990); Blank and Hollik (1998); Chen  _et al._ (2006); Fileviez Perez  _et
al._ (2009). Recently, it was pointed out that it could obviously explain the
CDF result see e.g. Strumia (2022); Perez  _et al._ (2022). However, grand
unification, besides providing the rationale for its existence and its
accessibility at today’s energies, also makes the low-energy effective theory
of the triplet completely predictive. It turns out that, in the leading
approximation, all the relevant physical processes of new particles depend
only on a single dimensionless parameter, the ratio of the possible $W$-mass
deviation over its mass. The CDF measurement, if confirmed, would fix this
parameter, making the theory testable at hadron colliders. This situation is
analogous to the SM Higgs case, where the knowledge of its mass suffices to
predict all of its physics. It is noteworthy that without the grand unified
framework, the same is not true. Some of the triplet couplings happen to be
arbitrary, to the extent that its neutral component could even be artificially
decoupled Perez  _et al._ (2022) from SM fermions. In our case, these
couplings are uniquely determined and related to the couplings with gauge
bosons.

We show here how this remarkable result comes about - but first we review
briefly the unification conditions of the model.

$\mathbf{SU(5)}$ and the oasis Let us start with a discussion of gauge
coupling unification within the SM. It is well-known that the three gauge
couplings do not unify at a single point: $\alpha_{1}$ meets $\alpha_{2}$ too
early, at about $10^{13}\rm GeV$, while $\alpha_{2}$ and $\alpha_{3}$ meet at
around $10^{17}\rm GeV$.

In order to keep proton stable enough, the only way out of this impasse is to
increase the $\alpha_{1}-\alpha_{2}$ meeting point. Since there are no new
gauge bosons that could slow down the rise of $\alpha_{1}$, the only option is
to slow down $\alpha_{2}$ \- through new fermions or/and scalars with small
(ideally zero) hypercharge. Clearly, the real $SU(2)_{L}$ scalar triple
$T_{S}$ does the needed job, since it carries zero hypecharge. In this sense,
it is an ideal candidate.

However, even when light, not far from the electro-weak scale, this particle
does not suffice to ensure gauge coupling unification. On top of that, the
minimal model Georgi and Glashow (1974) predicts a massless neutrino. More is
needed.

A simple addition is the adjoint fermion representation $24_{F}$ Bajc and
Senjanović (2007) which contains a real weak triplet fermion $F$, leading to
the following particle content

$24_{H};\quad 24_{F};\quad 5_{H};\quad\overline{5}_{F}^{a};\quad 10_{F}^{a},$
(1)

where $a=1,2,3$ denotes the generation index. As well known, compatibility
with the SM fermion spectrum requires the presence of higher dimensional
operators.

The scalar adjoint representation $24_{H}$ provides the GUT symmetry breaking
to the SM and that’s where the scalar triplet $T_{S}$ resides, while the
triplet $T_{F}$ belongs to the fermion adjoint representation $24_{F}$. As
already mentioned above, these triplet states allow for the successful
unification of gauge couplings as can be seen from the unification condition
Bajc and Senjanović (2007)

$\frac{M_{G}}{M_{W}}=\exp{\left[\frac{15\pi}{42}(\alpha_{1}^{-1}-\alpha_{2}^{-1})_{M_{W}}\right]\left(\frac{M_{W}^{5}}{m_{T_{S}}m_{T_{F}}^{4}}\right)^{\frac{5}{84}}}.$
(2)

Proton stability $\tau_{p}\gtrsim 10^{34}\rm yrs$ Takenaka  _et al._ (2020)
requires $M_{G}\gtrsim 10^{15}\rm GeV$, and thus both $T_{S}$ and $T_{F}$ need
to lie close to present-day energies, below $10\rm TeV$.

In the rest of this work, we focus on the scalar triplet $T_{S}$, whose
physics has not yet been fully explored. Before plunging into its interesting
phenomenology, we need check the consistency of its mass lying orders of
magnitude below the GUT scale. For this, we turn our attention to the Higgs
potential. It is easy to see that the relevant terms are given by

$V_{\rm relevant}=m^{2}5_{H}^{\dagger}5_{H}+M5_{H}^{\dagger}24_{H}5_{H}+\beta
5_{H}^{\dagger}24_{H}^{2}5_{H}.$ (3)

This potential serves to ensure the following physical conditions: the masses
of the SM Higgs doublet $\Phi$ and the scalar triplet $T_{S}$ must be made to
vanish at this stage, while at the same time, the colored triplet $T_{C}$ from
$5_{H}$ must lie close to the GUT scale, since it mediates proton decay. Even
without the new light triplet, this leads to the infamous doublet-colored
triplet splitting problem, which requires huge fine-tuning.

To be precise, the standard argument of fine-tuning, which we follow here,
assumes no judicious cancellations of Yukawa coupling matrices, which could
lower arbitrarily color triplet mass Dvali (1992). In this case, $T_{C}$ could
lie even at TeV energies, with a potentially interesting collider
phenomenology Dvali (1992); Dvali  _et al._ .

We have nothing new to say on this subject, only that the lightness of $T_{S}$
requires additional fine tuning. This can be easily achieved in the case of
the general potential for $24_{H}$ Guth and Weinberg (1981). There is more to
it, though, since $T_{S}$ has cubic coupling $\mu\Phi^{\dagger}T_{S}\Phi$,
that must be also kept small in order to preserve the stability of the SM
symmetry breaking. We show now how this works. In general there are higher-
dimensional terms in the potential, but for the sake of simplicity and
illustration, we ignore them (including them only facilitate things).

Upon breaking of $SU(5)$ gauge symmetry with $\langle
24_{H}\rangle=v_{G}\,{\rm diag}\left(1,1,1,-3/2,-3/2\right)$, the relevant
mass terms become

$\begin{split}&m_{\Phi}^{2}=m^{2}-\frac{3}{2}Mv_{G}+\frac{9}{4}\beta
v_{g}^{2}\\\ &m_{T_{C}}^{2}=m^{2}+Mv_{G}+\beta v_{G}^{2}\\\ &\mu=M+3\beta
v_{G}.\end{split}$ (4)

At this level, before turning the SM symmetry breaking on, the weak doublet
and triplet masses must effectively vanish: $\mu\simeq 0$ and $m_{\Phi}\simeq
0$. This leads to the mass for the heavy colored triplet
$m_{T_{C}}^{2}=25\beta v_{G}^{2}/4$, which tells us that $\beta$ must be
positive, with $\beta\gtrsim 10^{-3}$ on order to keep $T_{C}$ heavy enough
for the sake of proton’s stability.

In short, as in the case of the minimal $SU(5)$ theory, one can achieve the
desired spectra by fine tuning the couplings in the potential. While this may
not be appealing, it is perfectly legitimate and simply a must for the sake of
the consistency of the theory.

Low energy effective theory It has been known for decades Buras  _et al._
(1978) that the scalar triplet $T_{S}=T$ attains a non-vanishing VEV $v_{T}$
due to the tadpole interaction with the SM Higgs doublet $\Phi$. To see this
explicitly, let us write the Lagrangian

$\begin{split}&\mathcal{L}=|D_{\mu}\Phi|^{2}+{\rm
Tr}(D_{\mu}T)^{2}+m_{\Phi}^{2}\Phi^{\dagger}\Phi-m_{T}^{2}{\rm Tr}\,T^{2}\\\
&-\lambda_{T}\left(\rm{Tr}\,T^{2}\right)^{2}-\lambda_{\Phi}(\Phi^{\dagger}\Phi)^{2}-\rho\Phi^{\dagger}\Phi{\rm
Tr}\,T^{2}-\mu\Phi^{\dagger}T\Phi,\end{split}$ (5)

where the hermitian matrix $T$ transforms as $T\rightarrow UTU^{\dagger}$ and
its covariant derivative is given by
$D_{\mu}T=\partial_{\mu}T+ig\left[A_{\mu},T\right]$. These fields are
decomposed as

$\Phi=\frac{1}{\sqrt{2}}\begin{pmatrix}\sqrt{2}\phi^{+}\\\
v+\phi_{R}^{0}+i\phi_{I}^{0}\end{pmatrix},\,T=\frac{1}{2}\begin{bmatrix}v_{T}+t^{0}&\sqrt{2}t^{+}\\\
\sqrt{2}t^{-}&-v_{T}-t^{0}\end{bmatrix}.$ (6)

Once the $SU(2)$ breaking is turned on, the $\mu$-term provides a tadpole for
$T$ which therefore develops a non-vanishing expectation value. From
(LABEL:eq:lagrangian), $v_{T}$ is readily computed

$v_{T}\simeq\frac{\mu}{g^{2}}\left(\frac{M_{W}}{m_{T}}\right)^{2},$ (7)

where we used the tree-level SM mass relation $M_{W}^{SM}\simeq gv/2$. In the
above expression, $m_{T}$ actually stands for the triplet mass, corrected by
the $\rho$ term. For simplicity, we use the same expression as in
(LABEL:eq:lagrangian).

Eq. (7) holds true if $v_{T}\ll M_{W}$. This fits nicely with the high-
precision success of the SM, which, as it is well-known, leads to
$v_{T}\lesssim\mathcal{O}({\rm GeV})\ll M_{W}$. In this case, the mass
spectrum for the triplet becomes

$\begin{split}m_{T}^{2}\simeq
m_{t^{\pm}}^{2}=m_{t^{0}}^{2}+\mathcal{O}(v_{T}^{2}).\end{split}$ (8)

Note that the masses in (8) are independent of the sign in front of
$m_{T}^{2}$ in (LABEL:eq:lagrangian). The phenomenological bound
$m_{t^{\pm}}\gtrsim 250\rm GeV$ Chiang  _et al._ (2021), forces the same bound
on $m_{T}$ since the last contribution in (8) is negligible. It follows then,
independently of whether or not the triplet is Higgsed111We are grateful to G.
Dvali for helping us clarify this essential point., that the mass spectrum is
dominated by the tadpole-coupling $\mu$, and is therefore degenerate.

Radiative corrections further split the masses, but are negligible Cirelli
_et al._ (2006) for the present discussion.

This is an example of the hierarchy of weak isospin breaking. In theories with
a large mass scale $M\gg M_{W}$, higher isospin fields have VEVs proportional
to $M_{W}^{n+2}/M^{n}$, where $n$ is an integer. This can be understood as the
decoupling of heavy SM singlet fields: when $M\rightarrow\infty$, all physical
effects must disappear. In the case of a triplet, $n=2$ and (7) is recovered
in the limit $\mu\simeq m_{T}=M$.

Another celebrated example of the hierarchy of weak isospin breaking is the
type II seesaw mechanism Magg and Wetterich (1980); Cheng and Li (1980);
Mohapatra and Senjanović (1981) in which case $T$ is traded for a
$SU(2)$-triplet with two units of hypercharge. Here also the VEV of the
charged triplet gets suppressed by its large mass which can be used nicely to
explain the smallness of neutrino mass. This, just as the usual type I seesaw
mechanism, is again an example of decoupling and it is represented by the
effective field theory argument Weinberg (1979). Interestingly, it is claimed
Heeck (2022) that this could also account for the CDF-result. Notice, however,
that there is no independent physical reason for this charged triplet to be
relatively light, contrary to what is being discussed here.

Once $v_{T}$ is turned on, the $\mu$-term leads to the physically important
mixing $\vartheta$ between $\Phi$ and $T$. By defining

$\vartheta=\frac{g\,v_{T}}{M_{W}},$ (9)

the physical states in the neutral sector become

$\begin{split}&h^{0}\simeq\phi^{0}+\vartheta\left(\frac{m_{H}^{2}}{m_{H}^{2}-m_{h}^{2}}\right)t^{0}\\\
&H^{0}\simeq
t^{0}-\vartheta\left(\frac{m_{H}^{2}}{m_{H}^{2}-m_{h}^{2}}\right)\phi^{0},\end{split}$
(10)

where $h^{0}$ and $H^{0}$ stand for the SM Higgs boson and the neutral
component of the triplet, respectively. In (10), $m_{h}=m_{\Phi}$ and
$m_{H}=m_{T}$ are the masses of the physical states. Similarly, for the
charged components, the physical states are

$\begin{split}&G^{+}\simeq\phi^{+}-\vartheta t^{+},\\\ &H^{+}\simeq
t^{+}+\vartheta\phi^{+}.\end{split}$ (11)

where $G^{+}$ is the would-be Goldstone boson eaten by $W^{+}$, and $H^{+}$ is
a heavy charged particle from the triplet. Recall that phenomenological
consistency requires $m_{H^{0}}\simeq m_{H^{+}}\simeq m_{H}$.

In (9) $\mu$ was traded for $v_{T}$ using (7). Eqs. (10) and (11) are given in
the physical limit of small $\vartheta$, justified by $v_{T}\ll M_{W}$. Apart
from a small phenomenological region for $m_{H}\gtrsim 250\rm GeV$,
corresponding to its lower bound, the mixings in (10) and (11) are, up to a
sign, identical. In what follows, we will give analytical formulas relevant
for most of the parameter space $500{\rm GeV}\lesssim m_{H}\lesssim 10\rm
TeV$. The small corrections in the lower end of the spectrum will be shown in
the Figs. 2-3 below. In what follows, we will continue working in the
excellent approximation of $\vartheta\ll 1$.

Strictly speaking, the mixing gets corrected by the $\rho$ term, through small
$v_{T}$. However, this contribution turns out to be small (by an order of
magnitude or more), since $\rho\ll 1$, dictated by the requirement that the
theory stays perturbative all the way to the unification scale. This can be
easily seen from the 1-loop renormalization group flow Cheng  _et al._ (1974)
of the quartic couplings in (LABEL:eq:lagrangian)

$\begin{split}&\frac{d\lambda_{\Phi}}{dt}\simeq\left(\frac{d\lambda_{\Phi}}{dt}\right)_{SM}+\frac{2}{3\pi^{2}}\rho^{2},\\\
&\frac{d\lambda_{T}}{dt}\simeq\frac{1}{16\pi^{2}}\left({24{g_{2}}^{4}-24{g_{2}}^{2}\lambda_{T}+11\lambda_{T}^{2}+16\rho^{2}}\right),\\\
&\frac{d\rho}{dt}\simeq\frac{1}{32\pi^{2}}\left(12\rho
y_{t}^{2}-3{g_{1}}^{2}\rho+6{g_{2}}^{4}-33{g_{2}}^{2}\rho\right.\\\
&\,\,\quad\qquad\qquad\qquad\qquad\left.+10\lambda_{T}\rho+24{\lambda_{\Phi}}\rho+16\rho^{2}\right).\end{split}$
(12)

where $t=\log\mu$, the first term on the right-hand side of the first equation
denotes the SM running of the quartc Higgs coupling and Yukawa couplings other
than the top one have been omitted, being small. This includes the Yukawa
couplings of $T_{F}$, which must be tiny in order to guarantee the smallness
of neutrino masses through the seesaw mechanism. Eqs. 12, together with the
gauge couplings and top Yukawa coupling renormalization group equations were
then solved numerically between $M_{W}$ and $M_{G}$.

First of all, $\rho$ increases the value of $\lambda_{\Phi}$ at $M_{G}$,
which, in the SM, is known to be negative at high energies Degrassi  _et al._
(2012). While the inclusion of uncertainties in the SM input parameters
combined with a higher-loop analysis seem to somehow alleviate the issue
Degrassi  _et al._ (2012) \- bringing the coupling close to zero - it is
remarkable that the triplet can easily ease this tension, as can be seen from
the first equation in (12).

The message is clear: the central value of the quartic Higgs coupling turns
negative at around $10^{10}\rm GeV$ if no other new particles contribute to
the running. This poses a big challenge to the desert picture, providing yet
another rationale for the oasis near $M_{W}$.

Figure 1: Allowed region from the 1-loop renormalization group flow.

Next, there is a condition $\lambda_{T}(M_{G})\lesssim 1$, in order to
preserve perturbativity at $M_{G}$. Moreover, vacuum stability at the
electroweak scale requires $\lambda_{T}(M_{W})>0$ and, for negative $\rho$,
$\rho(M_{W})^{2}<4\lambda_{\Phi}(M_{W})\lambda_{T}(M_{W})$, implying
$-0.2\lesssim\rho(M_{G})\lesssim 0.6$, or $|\rho(M_{W})|\lesssim 0.3$ as shown
in Fig. 1. This translates into a correction to $\vartheta$ in (10) which is
negligible for the triplet mass above $500$ GeV, as shown in Fig. 2 below. The
rotation to the physical basis in the charged sector (11) is unaffected by
$\rho$-value.

Finally, the $\lambda_{T}$ term is negligible for tiny $v_{T}$. This shows
that, at leading order, the only relevant coupling in (LABEL:eq:lagrangian) is
the $\mu$-term. In what follows, our results will be presented in this limit.

The end result is striking. The apparently complicated structure in
(LABEL:eq:lagrangian), gets tremendously simplified. The only physical
parameters are $\vartheta$ and $m_{H}$, which tells us that $\vartheta$ will
govern all the physical processes involving the triplet. This will be
confirmed by explicit calculations below.

The first and foremost impact of $v_{T}$, or better $\vartheta$, is the
modification of the $W$-mass

$M_{W}^{2}\simeq(M_{W}^{SM})^{2}(1+\vartheta^{2}),\quad
M_{Z}^{2}=(M_{Z}^{SM})^{2},$ (13)

where tiny loop effects from $T$ are neglected. Note that the $Z$-mass is
unchanged since $T$ carries no hypercharge. If one accepts the CDF result
$M_{W}=80.433\pm 0.009$ Aaltonen  _et al._ (2022), by taking
$M_{W}^{SM}=80.357\pm 0.006\rm GeV$ Zyla  _et al._ (2020), for
$g(M_{W})=0.6517$ Zyla  _et al._ (2020), it follows from (13) that
$\vartheta\simeq 0.04$. We should stress that the degeneracy of charged and
neutral triplet components masses plays a crucial role in this - otherwise
there would be also the usual radiative correction to the $W$ and $Z$-mass
ratio, which would complicate things.

This is the key point we want to convey: minimal $SU(5)$ makes a strong case
for the $W$-mass deviation from its SM value. Even more excitingly, the
resulting low-energy effective theory is completely predictive.

Before turning to the phenomenological implications of the theory, we provide
the relevant couplings of $H^{0}$ and $H^{+}$ with the SM particles. For the
sake of clarity, and illustration, we write these coupling for $500{\rm
GeV}\lesssim m_{H}\lesssim 10\rm TeV$, where as we argued above, the mixings
in the neutral and charged sectors are basically equal

$\begin{split}&g\vartheta H^{0}:\\\
&\quad\,\,\,M_{W}WW,\,\frac{M_{Z}Z^{2}}{2c_{W}},\,\frac{(m_{H}h^{0})^{2}}{4M_{W}},\frac{m_{f}\overline{f}f}{2M_{W}},\\\
&g\vartheta H^{+}:\\\
&\quad\,\,\,\frac{M_{W}ZW^{-}}{c_{W}},\,\frac{i}{2}\overset{\leftrightarrow}{\partial}h^{0}W^{-},\frac{m_{d}\overline{u}_{L}d_{R}-m_{u}\overline{u}_{R}d_{L}}{\sqrt{2}M_{W}}.\end{split}$
(14)

Phenomenological implications The triplet $T$, as given in (6), contains a
real CP-even field $H^{0}$ and the charged field $H^{+}$. They can be produced
through the $W$ and $Z$ fusion, both singly (suppressed by the small $v_{T}$)
and pairwise. $H^{+}$ can obviously also be produced pairwise through the
photon and the $Z$-boson exchange. Since this is equally applicable to the
heavy triplet situation as well as the Higgsed version, one can profit from
the detailed studies of the latter case Gunion  _et al._ (1990); Blank and
Hollik (1998); Chen  _et al._ (2006); Fileviez Perez  _et al._ (2009).

Since the production rates are dominated by gauge couplings and do not involve
the $\vartheta$ parameter, we do not describe them here. Rather, we turn our
attention to their decays to SM particles, which in our case, arise through
$v_{T}$, or equivalently through the mixing $\vartheta$. The knowledge of
$\vartheta$ determines and correlates all the relevant decay rates as a
function of a single parameter $m_{H}$. We exemplify this on some selected
processes, working in the physical limit $m_{H}\gg M_{W}$.

## Neutral Component

Once produced, it decays into pairs of $W$, $Z$ and Higgs bosons, with the
rates computed from (14)

$\begin{split}{\frac{1}{2}}\Gamma(H^{0}\rightarrow&W^{+}W^{-})\simeq\Gamma(H^{0}\rightarrow
ZZ)\\\ &\simeq\Gamma(H^{0}\rightarrow
h^{0}h^{0})\simeq\vartheta^{2}\frac{g^{2}}{128\pi}\frac{m_{H}^{3}}{M_{W}^{2}}.\qquad\end{split}$
(15)

The monochromatic final states, especially $W$ and $Z$-bosons, provide clean
signature for the $H^{0}$ detection. It is worth noting that we are dealing
with a rather narrow width particle. For example, for $m_{H}\simeq{\rm TeV}$
and the maximum allowed value $\vartheta=0.04$, $\Gamma_{H}/M_{H}\simeq
10^{-3}$, ideal from the experimental point of view. If the triplet mass is
pushed to the maximally allowed value $10\rm TeV$, this ratio becomes a few
percent, similar to the $W$ and $Z$-bosons, still quite narrow. Obviously, for
smaller $\vartheta$, this is even more striking.

It is interesting that the gauge boson and Higgs final state rates are equal,
as expected in the high-energy limit $E\gg M_{W}$, when the SM symmetry
breaking becomes negligible. In this regime, the symmetry tells us that the
decay rate into the charged scalar (eaten by $W$) and the neutral ones ought
to be the same up to a factor of 2 (counting of degrees of freedom).

Next, using (14), the decay rate of $H^{0}$ into fermion-antifermion pair is
given by

$\Gamma(H^{0}\rightarrow
f\overline{f})\simeq{N_{c}}\vartheta^{2}\frac{g^{2}}{{32}\pi}\frac{m_{f}^{2}m_{H}}{M_{W}^{2}},$
(16)

where $N_{c}$ denotes the number of colours of the decay product.

Figure 2: Decay rates of neutral component for $H^{0}$, $\vartheta\simeq
0.04$.

In Fig. 2, the full decay rates of $\Gamma(H^{0}\rightarrow ZZ)$,
$\Gamma(H^{0}\rightarrow h^{0}h^{0})$ and $\Gamma(H^{0}\rightarrow
t\overline{t})$ are shown as a function of $m_{H}$. The filled area indicates
the uncertainty in $\vartheta$ due to the presence of $\rho$-term in the
potential of (LABEL:eq:lagrangian). As it can be seen, already for
$m_{H}\gtrsim 500\rm GeV$, this area collapses into a line.

Clearly, the fermion modes are negligible, except for the top quark, in which
case (16) is comparable to (15) for a small region around $m_{H}\simeq 500\rm
GeV$. As noted above, the $hh$ and $ZZ$ decay rates become practically
identical for $m_{H}\gtrsim 1\rm TeV$.

Since the SM Higgs boson decays into two photons, so does $H^{0}$, with the
rate suppressed by $\vartheta$

$\begin{split}\Gamma(H^{0}\rightarrow\gamma\gamma)\sim\vartheta^{2}\frac{m_{H}}{M_{W}}\,\Gamma(h^{0}\rightarrow\gamma\gamma).\end{split}$
(17)

This expression is just a rough dimensional estimate. In the SM, as we know,
$\gamma\gamma$ channel was the path towards the Higgs boson discovery: the
small rate was compensated by the clear signature of two monochromatic
photons, and by the fact that the SM Higgs cannot decay into two on shell
$W$-bosons. In this case, the effect is less important since the on-shell $W$
and $Z$ final states decay channels are now allowed and enhanced by large
$m_{H}$. Still, it is a clean channel, worth keeping in mind.

## Charged Component

The charged component $H^{+}$ decays mainly into a $WZ$ and $Wh$, with the
following respective decay rates computed from (14)

$\begin{split}\Gamma(H^{+}\rightarrow W^{+}Z)&\simeq\Gamma(H^{+}\rightarrow
W^{+}h^{0})\\\ &\simeq\Gamma(H^{0}\rightarrow W^{+}W^{-}).\end{split}$ (18)

In analogy with (16), the fermionic decay rates of $H^{+}$ are negligible,
except for $\Gamma(H^{+}\rightarrow t\overline{b})$ which can be computed from
(14) and gives

$\Gamma(H^{+}\rightarrow t\overline{b})\simeq\Gamma(H^{0}\rightarrow
t\overline{t}).$ (19)

Notice that since $H_{0}$ and $H^{+}$ are degenerate, all the decay rates are
correlated. This can be seen in Fig. 3. The situation is analogous to the
neutral scalar case, except that, for $m_{H}\lesssim 500\rm GeV$, the
fermionic decay dominates, since the kinematics helps (bottom quark is almost
massless at these energies). For higher values, once again, the decays into
bosons take over.

Figure 3: Decay rates of charged component $H^{+}$, $\vartheta\simeq 0.04$.

## Standard Model Higgs

Due to the mixing between the doublet and the triplet, the usual Higgs boson
couplings are modified by $\vartheta^{2}$, implying a small deviation from its
couplings to gauge bosons and fermions, on the order of a permil. While
difficult to measure, especially at the hadron colliders, this would
eventually serve as an important independent test of the theory.

Summary In this work we have revisited a well-motivated extension of the
minimal $SU(5)$ grand unified theory. The theory predicts a light fermion and
a light scalar real weak triplet Bajc and Senjanović (2007). The fermion
triplet plays a key-role in generating a small neutrino mass through the
seesaw mechanism and its physics has been studied extensively in the past.

Here we focus on its scalar counterpart. It naturally develops a small
expectation value, implying generically a deviation of $W$-mass from its SM
value. The emerging low-energy effective theory is then truly predictive:
there are a plethora of calculable physical processes - which depend only on
the triplet mass and a single parameter $\vartheta$ that measures the $W$-mass
deviation.

In view of the recent CDF announcement Aaltonen  _et al._ (2022) of the
$W$-mass disagreement with its SM value, this theory potentially becomes of
great phenomenological importance. If true, the CDF result would fix
$\vartheta$, making the theory completely predictive, and easily testable at
the LHC or next hadron collider. If the triplet mass was to be pushed into the
TeV region, this would become a high-precision theory, almost at the level of
the SM.

It is noteworthy that the minimal $SO(10)$ grand unified theory based on small
Higgs representations Preda  _et al._ (2022), also predicts, among other
states, the same light scalar triplet. If it was not for other light
multiplets - required by the consistency of the theory - the same low-energy
predictivity would emerge. However, the complexity of the resulting
phenomenology calls for a separate detailed study Preda  _et al._ .

Defying the decades-long prejudice of the existence of a desert at energies
above $M_{W}$, this simple and well-motivated grand unified theory provides a
nearby oasis. The end result is profound: the theory makes the case for the
$W$-mass deviation and ties it to a number of calculable physical processes.
The far reaching testable consequences may have just started to unravel
themselves.

Note added As we were about to post this work, a new paper Evans  _et al._
(2022) appeared in which a particular $SU(5)$ model - very different from the
theory presented here - is built in order to accommodate the CDF-anomaly.

Acknowledgments We are grateful to Damir Lelas and Gia Dvali for illuminating
discussions, and to Damir Lelas and Alejandra Melfo for important comments and
indispensable help in the preparation of this manuscript.

## References

  * Bajc and Senjanović (2007) B. Bajc and G. Senjanović, JHEP 08, 014 (2007), arXiv:hep-ph/0612029 .
  * Georgi and Glashow (1974) H. Georgi and S. L. Glashow, Phys. Rev. Lett. 32, 438 (1974).
  * Bajc  _et al._ (2007) B. Bajc, M. Nemevšek, and G. Senjanović, Phys. Rev. D 76, 055011 (2007), arXiv:hep-ph/0703080 .
  * Arhrib  _et al._ (2010) A. Arhrib, B. Bajc, D. K. Ghosh, T. Han, G.-Y. Huang, I. Puljak, and G. Senjanovic, Phys. Rev. D 82, 053004 (2010), arXiv:0904.2390 [hep-ph] .
  * Aaltonen  _et al._ (2022) T. Aaltonen _et al._ (CDF), Science 376, 170 (2022).
  * Ross and Veltman (1975) D. A. Ross and M. J. G. Veltman, Nucl. Phys. B 95, 135 (1975).
  * Gunion  _et al._ (1990) J. F. Gunion, R. Vega, and J. Wudka, Phys. Rev. D 42, 1673 (1990).
  * Blank and Hollik (1998) T. Blank and W. Hollik, Nucl. Phys. B 514, 113 (1998), arXiv:hep-ph/9703392 .
  * Chen  _et al._ (2006) M.-C. Chen, S. Dawson, and T. Krupovnickas, Phys. Rev. D 74, 035001 (2006), arXiv:hep-ph/0604102 .
  * Fileviez Perez  _et al._ (2009) P. Fileviez Perez, H. H. Patel, M. J. Ramsey-Musolf, and K. Wang, Phys. Rev. D 79, 055024 (2009), arXiv:0811.3957 [hep-ph] .
  * Strumia (2022) A. Strumia, (2022), arXiv:2204.04191 [hep-ph] .
  * Perez  _et al._ (2022) P. F. Perez, H. H. Patel, and A. D. Plascencia, (2022), arXiv:2204.07144 [hep-ph] .
  * Takenaka  _et al._ (2020) A. Takenaka _et al._ (Super-Kamiokande), Phys. Rev. D 102, 112011 (2020), arXiv:2010.16098 [hep-ex] .
  * Dvali (1992) G. R. Dvali, Phys. Lett. B 287, 101 (1992).
  * (15) G. Dvali, O. Sakhelashvili, and A. Stuhlfauth, in preparation .
  * Guth and Weinberg (1981) A. H. Guth and E. J. Weinberg, Phys. Rev. D 23, 876 (1981).
  * Buras  _et al._ (1978) A. J. Buras, J. R. Ellis, M. K. Gaillard, and D. V. Nanopoulos, Nucl. Phys. B 135, 66 (1978).
  * Chiang  _et al._ (2021) C.-W. Chiang, G. Cottin, Y. Du, K. Fuyuto, and M. J. Ramsey-Musolf, JHEP 01, 198 (2021), arXiv:2003.07867 [hep-ph] .
  * Cirelli  _et al._ (2006) M. Cirelli, N. Fornengo, and A. Strumia, Nucl. Phys. B 753, 178 (2006), arXiv:hep-ph/0512090 .
  * Magg and Wetterich (1980) M. Magg and C. Wetterich, Phys. Lett. B 94, 61 (1980).
  * Cheng and Li (1980) T. P. Cheng and L.-F. Li, Phys. Rev. D 22, 2860 (1980).
  * Mohapatra and Senjanović (1981) R. N. Mohapatra and G. Senjanović, Phys. Rev. D 23, 165 (1981).
  * Weinberg (1979) S. Weinberg, Phys. Rev. Lett. 43, 1566 (1979).
  * Heeck (2022) J. Heeck, (2022), arXiv:2204.10274 [hep-ph] .
  * Cheng  _et al._ (1974) T. P. Cheng, E. Eichten, and L.-F. Li, Phys. Rev. D 9, 2259 (1974).
  * Degrassi  _et al._ (2012) G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori, and A. Strumia, JHEP 08, 098 (2012), arXiv:1205.6497 [hep-ph] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020), and 2021 update.
  * Preda  _et al._ (2022) A. Preda, G. Senjanović, and M. Zantedeschi, (2022), arXiv:2201.02785 [hep-ph] .
  * (29) A. Preda, G. Senjanović, and M. Zantedeschi, in preparation .
  * Evans  _et al._ (2022) J. L. Evans, T. T. Yanagida, and N. Yokozaki, (2022), arXiv:2205.03877 [hep-ph] .


# Stokes drift and its discontents

Jacques Vanneste1, William R. Young2

(1School of Mathematics and Maxwell Institute for Mathematical Sciences,
University of Edinburgh, Edinburgh EH9 3FD, UK
2Scripps Institution of Oceanography, University of California at San Diego,
La Jolla CA 92093-0213, USA)

###### Abstract

The Stokes velocity $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$, defined
approximately by Stokes (1847, Trans. Camb. Philos. Soc., 8, 441–455), and
exactly via the Generalized Lagrangian Mean, is divergent even in an
incompressible fluid. We show that the Stokes velocity can be naturally
decomposed into a solenoidal component,
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$, and a remainder that
is small for waves with slowly varying amplitudes. We further show that
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ arises as the sole
Stokes velocity when the Lagrangian mean flow is suitably redefined to ensure
its exact incompressibility. The construction is an application of Soward &
Roberts’s glm theory (2010, J. Fluid Mech., 661, 45–72) which we specialise to
surface gravity waves and implement effectively using a Lie series expansion.
We further show that the corresponding Lagrangian-mean momentum equation is
formally identical to the Craik–Leibovich equation with
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ replacing
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$, and we discuss the form of the
Stokes pumping associated with both $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$
and $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$.

##  1 Introduction

Surface gravity waves induce a rectified motion of fluid particles and thus a
wave-averaged difference between the mean Eulerian velocity,
$\bm{u}^{\mathrm{\scriptscriptstyle{E}}}$, and the mean Lagrangian velocity
$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ [25, 26]

$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\bm{u}^{\mathrm{\scriptscriptstyle{E}}}+\bm{u}^{\mathrm{\scriptscriptstyle{S}}}\,.$
(1)

Above, $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ is the Stokes velocity, also
known as Stokes drift. The Stokes velocity is fundamental to understanding
wave-averaged effects, such as the Craik–Leibovich vortex force and the
Stokes–Coriolis force, in the wave-averaged momentum and vorticity equations
[4, 14, 8, 13, 17].

The Stokes velocity can be defined exactly at finite wave amplitude using
Generalized Lagrangian Mean (GLM) theory [1, 2]. This exact GLM
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ is rotational and compressible, even
if the underlying fluid motion is irrotational and incompressible [16].
Expansion in powers of a wave-amplitude parameter $\epsilon$ produces the
standard approximation [22, 26] to the Stokes velocity

$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\,\bm{u}_{1}}\,.$
(2)

The overbar in (2) denotes a running time mean, or phase average. In (2),
$\bm{u}_{1}$ is the linear (first order in $\epsilon$) velocity of the wave
and the associated displacement $\bm{\xi}_{1}$ is defined by
$\partial_{t}\bm{\xi}_{1}=\bm{u}_{1}$ and $\overline{\bm{\xi}_{1}}=0$. (The
subscript $1$ indicates first-order fields throughout.) The small-amplitude
approximation to $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in (2) is also
rotational and compressible: assuming only that
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{\xi}_{1}=0$, McIntyre [16] shows from (2) that

$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\partial_{t}\left(\,\tfrac{1}{2}\overline{\xi_{1i}\xi_{1j}}\,\right)_{,ij}\,.$
(3)

The time derivative of an averaged quadratic quantity in (3) entails the same
slow-modulation assumption that underlies the concept of group velocity and so
introduces a second small parameter, $\mu$. The Eulerian mean velocity,
$\bm{u}^{\mathrm{\scriptscriptstyle{E}}}$ in (1), is incompressible and thus
the divergent $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in (3) implies a
divergent Lagrangian mean velocity.

Figure 1: The sea-surface displacement, $s(x,t)$, of a packet of surface
gravity waves at $t=0$. The envelope is Gaussian,
$a\exp[-(x-x_{0}-\tfrac{1}{2}ct)^{2}/2\ell^{2}]$, and the carrier wavenumber
is $k=2\pi/100$m. The $100$ m wave length corresponds to an 8 second period
and a group velocity $c/2$ of $6.24$ m s-1. The maximum surface displacement
$a=1.27$ m corresponds to maximum wave orbital speed $1$ m s-1. The modulation
parameter is $\mu=1/(k\ell)=0.05$ and the wave slope is $\epsilon=ka=0.08$.

In figure 1 we illustrate the role of the two small parameters $\epsilon$ and
$\mu$, by considering a weakly nonlinear, slowly modulated two-dimensional
packet of deep-water surface gravity waves. The Stokes expansion [1,10] is
justified by the weakly-nonlinear assumption that the wave slope is small:
$\epsilon=ak\ll 1$, where $k$ is the wavenumber and $a$ is the amplitude of
the surface displacement. In this example the slow-modulation parameter is
$\mu=(k\ell)^{-1}\ll 1$ where $\ell$ is length scale of the packet envelope.
Figure 2 shows the motion of a fluid particle in the velocity field of this
wave.

Figure 2: Trajectory of a fluid particle in the linear velocity field of the
wave packet in figure 1. The top panels show the $x$\- and $z$\- displacements
as functions of time; the bottom panel shows the trajectory. In this
computation we assume that the depth $d$ is much greater than the packet
length scale $\ell$ so that the second-order Eulerian mean flow is negligible
in the wave-active zone.

Despite (3), and the exact results provided by GLM, some authors are reluctant
to accept the reality of non-zero
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$. Moreover
discontent with
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}\neq 0$ is
sometimes confounded with unease over the vertical component of the Stokes
velocity,
$w^{\mathrm{\scriptscriptstyle{S}}}=\bm{\hat{z}}\bm{\cdot}\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$.
For example, rather than taking the vertical component of (2), McWilliams,
Restrepo & Lane [18] defines a “vertical Stokes pseudo-velocity” which,
together with the horizontal components of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$, makes an incompressible three-
dimensional “Stokes pseudo-velocity”.

Mellor [19], while emphasizing that $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$
is divergent, is unwilling to accept a non-zero vertical component
$w^{\mathrm{\scriptscriptstyle{S}}}$: “a mean vertical drift is not
acceptable”. In view of concerns with the mean vertical drift
$w^{\mathrm{\scriptscriptstyle{S}}}$ it is reassuring that particle tracking
velocimetry can be used to observe vertical Lagrangian displacements,
$w^{\mathrm{\scriptscriptstyle{S}}}\neq 0$, beneath groups of deep-water waves
[27, 20]: the mean vertical drift is upward as a wave packet arrives and
downwards as the packet departs. (For the wave packet in figure 1, the maximum
vertical displacement resulting from $w^{\mathrm{\scriptscriptstyle{S}}}$ is
about $4.8$cm – this is not visible in figure 2.) After the passage of the
packet a fluid particle returns to its initial depth. It is a net vertical
displacement that is unacceptable: transient vertical motion, on time scales
longer than a ten-second wave period, and shorter than the hundred-second
packet transit time, is not a concern.

Mellor critiques the Craik & Leibovich ([3,4] CL hereafter) vortex-force
formulation of wave-mean interaction by arguing that CL and subsequent authors
incorrectly assume that the divergence of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ is zero. A different interpretation
is that CL and many followers assume that the wave field has no temporal
modulation so that the right of (3) is conveniently zero. For example,
McWilliams & Restrepo [17] claims to prove
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=0$. But
examination of this argument shows that [17] assumes that there is no temporal
modulation of the wave field. This raises the issue of whether the CL
formulation is incomplete or misleading in situations with temporal modulation
of the wave field e.g. in ocean observations[23] and in modeling the growth of
swell [28].

In this paper, we revisit the concept of Stokes velocity. For surface gravity
waves we exhibit a natural Helmholtz decomposition of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in (2) and argue that the solenoidal
component, $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$, can
advantageously replace $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in most
situations. We emphasize that the familiar form (2) of the Stokes velocity is
not unique but depends on a specific definition of the Lagrangian-mean flow,
namely the GLM definition of Andrews & McIntyre [1]. An alternative
definition, proposed by Soward & Roberts [24] and closely related to classical
averaging and its Lie series implementation (e.g. [15, 21]), leads to a
solenoidal Lagrangian-mean velocity with
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ as the corresponding
Stokes velocity. This alternative definition, known as ‘glm’ but better
characterized as ‘solenoidal Lagrangian mean’, has the added benefit of
coordinate independence in any geometry, unlike standard GLM (see [6] for
other coordinate-independent definitions of the Lagrangian mean). We show
that, for surface gravity waves, the associated Lagrangian-mean momentum
equation governing the dynamics of the Eulerian mean flow is the CL equation
with $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ replacing
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$.

The difference between GLM and glm is vividly illustrated with an example
proposed by O. Bühler (personal communication). Consider a bucket of initially
motionless water. If the water is agitated, for example by pressure forcing at
the surface, then the potential energy of the water is increased, or
equivalently the center of mass of the water is elevated above its initial
height. The fluid at the bottom of Bühler’s bucket, however, cannot move in
the vertical and so any definition of “Lagrangian mean” that tracks the
position of the center of mass of the fluid – such as GLM – will be divergent
in this situation, even though the velocity of the water is entirely
incompressible. Conversely, any definition of “Lagrangian mean” resulting in a
strictly solenoidal Lagrangian mean velocity – such as glm – cannot track the
center of mass of the fluid.

The plan of the paper is as follows. In §2, we sketch the derivation of the
standard form (2) of the Stokes velocity, give its Helmholtz decomposition,
and show how a simple modification of this derivation, implementing an
alternative Lagrangian-mean flow definition, naturally brings about the
velocity $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$. In §3, we
examine the respective role of $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ in Stokes pumping,
which is the mechanism whereby the horizontal divergence of the Stokes
transport drives an Eulerian mean flow. In §4, we show how the glm approach
enables the systematic construction of solenoidal Lagrangian-mean and Stokes
velocities up to arbitrary algebraic accuracy in $\epsilon$. We explain how
Lie series provide both an interpretation and an efficient implementation of
this construction, and we derive the glm version of the CL equations. Section
5 gives conclusions.

##  2 The Stokes velocity $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and its
solenoidal part $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$

###  2.1 Derivation of the Stokes velocity

We start by recalling the traditional derivation of the Stokes velocity in
(2). The position $\bm{x}(t)$ of a fluid particle is determined by solving

$\frac{\mathrm{d}\bm{x}(t,\alpha,\epsilon)}{\mathrm{d}t}=\bm{u}\big{(}\bm{x}(t,\alpha,\epsilon),t,\alpha,\epsilon\big{)}\,,$
(4)

where $\alpha$ is a wave phase, regarded as an ensemble parameter. The fluid
velocity $\bm{u}(\bm{x},t,\alpha,\epsilon)$ is incompressible,
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}=0$, and has the form

$\bm{u}(\bm{x},t,\alpha,\epsilon)=\epsilon\bm{u}_{1}(\bm{x},t,\alpha)+\epsilon^{2}\bm{u}_{2}(\bm{x},t,\alpha)+\cdots,$
(5)

where $\epsilon$ is the wave amplitude parameter. The leading order term
$\bm{u}_{1}(x,t,\alpha)$ is a fast wavy flow, so $\overline{\bm{u}_{1}}=0$,
where the mean, denoted by the overbar, is an average over the phase $\alpha$.

To average the fast wave oscillations in (4) we consider the ansatz

$\displaystyle\bm{x}(t,\alpha,\epsilon)$
$\displaystyle=\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon)+\bm{\xi}\big{(}\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon),t,\alpha,\epsilon\big{)}\,,$
(6a)
$\displaystyle=\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon)+\epsilon\bm{\xi}_{1}\big{(}\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon),t,\alpha\big{)}+\epsilon^{2}\bm{\xi}_{2}\big{(}\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon),t,\alpha\big{)}+\cdots,$
(6b)

where $\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon)$ is the slow motion
of a Lagrangian mean position. Think of
$\bm{x}^{\mathrm{\scriptscriptstyle{L}}}$ as a ‘guiding center’ such that
rapid wavy oscillations are confined to the displacements $\bm{\xi}_{n}$;
these displacements from $\bm{x}^{\mathrm{\scriptscriptstyle{L}}}$ do not grow
with time, i.e. all members of the ensemble remain close to the guiding center
$\bm{x}^{\mathrm{\scriptscriptstyle{L}}}$. The motion of
$\bm{x}^{\mathrm{\scriptscriptstyle{L}}}$ is written as

$\displaystyle\frac{\mathrm{d}\bm{x}^{\mathrm{\scriptscriptstyle{L}}}}{\mathrm{d}t}$
$\displaystyle=\epsilon^{2}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ (7)
$\displaystyle=\epsilon^{2}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}_{2}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t)+\epsilon^{3}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}_{3}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t)+\cdots\,,$
(8)

where $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}(\bm{x},t,\epsilon)$ is the
Lagrangian mean velocity, yet to be defined and determined. The ansatz (6b) is
ambiguous because requiring only that the $\bm{\xi}_{n}$’s do not grow with
time does not uniquely determine $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ and
$\bm{\xi}_{n}$. We return to this point below.

Substituting (6b) and (8) into (4) and (5) and matching powers of $\epsilon$
at the first two orders results in

$\displaystyle\partial_{t}\bm{\xi}_{1}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)$
$\displaystyle=\bm{u}_{1}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)\,,$
(9a)
$\displaystyle\bm{u}^{\mathrm{\scriptscriptstyle{L}}}_{2}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t)+\partial_{t}\bm{\xi}_{2}(t,\alpha)$
$\displaystyle=\bm{u}_{2}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)+\bm{\xi}_{1}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)\\!\bm{\cdot}\\!\bm{\nabla}\,\bm{u}_{1}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)\,.$
(9b)

Choosing to follow Stokes [25] and Andrews & McIntyre [1], one disambiguates
(6b) by requiring that
$\overline{\bm{\xi}_{2}(\bm{x}^{\mathrm{\scriptscriptstyle{L}}},t,\alpha)}=0$.
In this case averaging (9b) produces the familiar result

$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\underbrace{\overline{\bm{u}_{2}}}_{\bm{u}^{\mathrm{\scriptscriptstyle{E}}}}+\underbrace{\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\,\bm{u}_{1}}}_{\bm{u}^{\mathrm{\scriptscriptstyle{S}}}}\qquad\text{(with
$\,\overline{\bm{\xi}_{2}}=0$)}\,.$ (10)

In (10) we have now omitted the subscript 2 on
$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$. With (10) we recover $\eqref{eq2}$
and the small wave-amplitude version of (1).

###  2.2 The solenoidal Stokes velocity
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$

Let us examine the divergence of $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in
more detail. The ‘un-averaged Stokes velocity’ can be written exactly as

$\displaystyle\bm{\xi}_{1}\bm{\cdot}\bm{\nabla}\bm{u}_{1}$
$\displaystyle=\partial_{t}\tfrac{1}{2}(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}+\tfrac{1}{2}(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{u}_{1}-\tfrac{1}{2}(\bm{u}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}\,,$
(11a)
$\displaystyle=\partial_{t}\tfrac{1}{2}(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}+\bm{\nabla}\\!\times\\!\tfrac{1}{2}(\bm{u}_{1}\times\bm{\xi}_{1})\,.$
(11b)

In passing from (11a) to (11b) we have used wave incompressibility

$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}_{1}=\bm{\nabla}\\!\bm{\cdot}\\!\bm{\xi}_{1}=0\,,$
(12)

to simplify the standard vector identity for the curl of the cross product
$\bm{u}_{1}\times\bm{\xi}_{1}$. The average of (11b),

$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\underbrace{\partial_{t}\tfrac{1}{2}\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}}}_{O(\epsilon^{2}\mu)}+\underbrace{\bm{\nabla}\\!\times\\!\tfrac{1}{2}\left(\overline{\bm{u}_{1}\times\bm{\xi}_{1}}\right)}_{O(\epsilon^{2})}\,,$
(13)

identifies a solenoidal part of the Stokes velocity as

$\displaystyle\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$
$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\bm{\nabla}\\!\times\\!\left(\tfrac{1}{2}\,\overline{\bm{u}_{1}\times\bm{\xi}_{1}}\right)\,,$
(14a)
$\displaystyle=\tfrac{1}{2}\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\,\bm{u}_{1}}-\tfrac{1}{2}\overline{(\bm{u}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\,\bm{\xi}_{1}}\,.$
(14b)

The solenoidal vector $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$
is the incompressible part of the Stokes velocity for all types of weakly
nonlinear waves in an incompressible fluid. We propose that
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ can advantageously
replace the traditional form of the Stokes velocity (2) in many circumstances.

The solenoidal Stokes velocity arises naturally if a small change is made in
the derivation in §22.1: suppose we decide from the outset to change the
definition of ‘Lagrangian mean’ so that
$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ in (8) is solenoidal. To implement
this choice we disambiguate (6b) with
$\overline{\bm{\xi}_{2}}=\tfrac{1}{2}\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}}$
(rather than $\overline{\bm{\xi}_{2}}=0$). Averaging (9b) then results in

$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\underbrace{\overline{\bm{u}_{2}}}_{\bm{u}^{\mathrm{\scriptscriptstyle{E}}}}+\underbrace{\bm{\nabla}\\!\times\\!\left(\tfrac{1}{2}\,\overline{\bm{u}_{1}\times\bm{\xi}_{1}}\right)}_{\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}}\,,\qquad\text{(with
$\overline{\bm{\xi}_{2}}=\tfrac{1}{2}\overline{(\bm{\xi}_{1}\\!\bm{\cdot}\\!\bm{\nabla})\bm{\xi}_{1}}$).}$
(15)

The incompressible Lagrangian mean velocity in (15) is an alternative to the
traditional compressible $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ in (10).
The $O(\epsilon^{2})$ shift $\overline{\bm{\xi}_{2}}$ in the position of
guiding centre $\bm{x}^{\mathrm{\scriptscriptstyle{L}}}(t,\epsilon)$ produces
Lagrangian mean and Stokes velocities that are divergence free at order
$\epsilon^{2}$. In §4 we discuss a systematic framework – Soward &
Robert’s[24] glm alternative to GLM – that generalizes this property to
arbitrarily high order in $\epsilon$.

Note that the difference between $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ is a time derivative
so has no impact on particle dispersion for waves that are represented by
stationary random processes as discussed by Holmes-Cerfon and Bühler [3].

###  2.3 Irrotational linear waves

In addition to incompressibility in (12), surface gravity waves are also
irrotational,

$\bm{\nabla}\\!\times\\!\bm{u}_{1}=\bm{\nabla}\\!\times\\!\bm{\xi}_{1}=0\,.$
(16)

Using (16), the final term in the vector identity
$\bm{\xi}_{1}\bm{\cdot}\bm{\nabla}\bm{\xi}_{1}=\bm{\nabla}\tfrac{1}{2}|\bm{\xi}_{1}|^{2}+(\bm{\nabla}\\!\times\\!\bm{\xi}_{1})\times\bm{\xi}_{1}$
is zero and the average of (13) simplifies further to a Helmholtz
decomposition of the surface-wave Stokes velocity

$\displaystyle\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$
$\displaystyle=\partial_{t}\bm{\nabla}\tfrac{1}{4}\overline{|\bm{\xi}_{1}|^{2}}+\underbrace{\bm{\nabla}\\!\times\\!\tfrac{1}{2}\overline{\bm{u}_{1}\times\bm{\xi}_{1}}}_{\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}}\,.$
(17)

The divergence of the surface-wave Stokes velocity in (17) is

$\displaystyle\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$
$\displaystyle=\underbrace{\partial_{t}\triangle\,\tfrac{1}{4}\overline{|\bm{\xi}_{1}|^{2}}}_{O(\epsilon^{2}\mu)}\,,$
(18)

where $\triangle=\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2}$ is the
Laplacian. The expression for
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ in (18)
is simpler than that in (3), but restricted to waves with
$\bm{\nabla}\\!\times\\!\bm{\xi}_{1}=0$.

The curl of the Stokes velocities does not vanish:

$\displaystyle\bm{\nabla}\\!\times\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\bm{\nabla}\\!\times\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$
$\displaystyle=\bm{\nabla}\\!\times\\!\bm{\nabla}\\!\times\\!\tfrac{1}{2}\overline{\bm{u}_{1}\times\bm{\xi}_{1}}\,,$
(19)
$\displaystyle=-\triangle\tfrac{1}{2}\overline{\bm{u}_{1}\times\bm{\xi}_{1}}\,,$
(20)

where we have used
$\bm{\nabla}\\!\times\\!\bm{\nabla}\\!\times\\!=\bm{\nabla}\bm{\nabla}\\!\bm{\cdot}\\!-\triangle$
and
$\bm{\nabla}\\!\bm{\cdot}\\!\tfrac{1}{2}\overline{\bm{u}_{1}\times\bm{\xi}_{1}}=0$.
With an error of order $\mu^{2}\epsilon^{2}$ we can replace $\triangle$ with
$\partial_{z}^{2}$ in (20).

##  3 Stokes transport and Stokes pumping

Discussions of the deep return flow associated with a surface-gravity wave
packet argue [26] that:

> “ the return flow $\cdots$ can be explained as the irrotational response to
> balance the Stokes transport $\cdots$ that acts to ‘pump’ fluid from the
> trailing edge of the packet to the leading edge”.

Similar sentiments are expressed in [9, 7]. With this physical picture in
mind, it is instructive to compare the Stokes pumping associated with
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ with that of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$. The interpretation
of the return flow in the quotation above is most consistent with
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$.

###  3.1 Stokes pumping and
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}=\bm{\nabla}\\!\times\\!\tfrac{1}{2}\overline{\bm{u}_{1}\times\bm{\xi}_{1}}$

We start with the easy solenoidal case, with Stokes transport

$\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}(x,y,t)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int_{-\infty}^{0}\\!\\!\\!\\!\\!\\!u^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}\,\bm{\hat{x}}+v^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}\,\bm{\hat{y}}\,\mathrm{d}z\,,$
(21)

where the vertical integration above is from the bottom of the wave-active
zone (denoted $-\infty$) to the mean sea surface at $z=0$. (In this section we
confine attention to deep-water waves so that the lower limit, $-\infty$, is
well above the distant bottom.) Vertical integration of
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}=0$
over the wave-active region produces the unsurprising result

solenoidal Stokes pumping
$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\bm{\nabla}\\!\bm{\cdot}\\!\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}\,,$
(22a)
$\displaystyle=w^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}\big{|}_{0}\,,$
(22b)

where we use $\big{|}_{0}$ to denote evaluation at $z=0$ e.g.
$w^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}\big{|}_{0}=w^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}(x,y,0,t)$.
The result (22b) is consistent with the idea that horizontal convergence
within the wave active zone pumps fluid downwards, out of the wave-active
zone, with the vertical velocity
$w^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}(x,y,0,t)$.

Vertical integration of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}=\tfrac{1}{2}\overline{(\bm{\xi}_{1}\bm{\cdot}\bm{\nabla})\bm{u}_{1}}-\tfrac{1}{2}\overline{(\bm{u}_{1}\bm{\cdot}\bm{\nabla})\bm{\xi}_{1}}$
over the wave active zone results in

$\displaystyle\bm{\hat{x}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}$
$\displaystyle=\tfrac{1}{2}\big{(}\overline{\zeta_{1}u_{1}}-\overline{w_{1}\xi_{1}}\big{)}\big{|}_{0}+\partial_{y}\chi\,,$
(23a)
$\displaystyle\bm{\hat{y}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}$
$\displaystyle=\tfrac{1}{2}\big{(}\overline{\zeta_{1}v_{1}}-\overline{w_{1}\eta_{1}}\big{)}\big{|}_{0}-\partial_{x}\chi\,,$
(23b)

where

$\chi\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\tfrac{1}{2}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\left(\overline{\eta_{1}u_{1}}-\overline{v_{1}\xi_{1}}\right)\,\mathrm{d}z\,,$
(24)

is a streamfunction for horizontally circulating Stokes transport. This
horizontal Stokes circulation is previously unremarked, perhaps because the
terms involving $\chi$ in (23a) and (23b) are order $\mu$ smaller than the
other terms. Using the coordinate expressions in (23a) and (23b) the
solenoidal pumping can be expressed entirely in terms of surface quantities:

$\text{solenoidal Stokes
pumping}=\partial_{x}\tfrac{1}{2}\big{(}\overline{\zeta_{1}u_{1}}-\overline{w_{1}\xi_{1}}\big{)}\big{|}_{0}+\partial_{y}\tfrac{1}{2}\big{(}\overline{\zeta_{1}v_{1}}-\overline{w_{1}\eta_{1}}\big{)}\big{|}_{0}\,.$
(25)

###  3.2 Stokes pumping and
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\overline{(\bm{\xi}_{1}\bm{\cdot}\bm{\nabla})\bm{u}_{1}}$

We turn now to the traditional definition of the Stokes velocity. Define the
Stokes transport via

$\bm{T}^{\mathrm{\scriptscriptstyle{S}}}(x,y,t)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!u^{\mathrm{\scriptscriptstyle{S}}}\,\bm{\hat{x}}+v^{\mathrm{\scriptscriptstyle{S}}}\,\bm{\hat{y}}\,\mathrm{d}z\,,$
(26)

and the Stokes pumping by

$\text{Stokes
pumping}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\bm{\nabla}\\!\bm{\cdot}\\!\bm{T}^{\mathrm{\scriptscriptstyle{S}}}\,.$
(27)

Because
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}\neq 0$
there is no analogy to (22b). Instead, integrating (18) over the wave-active
zone

$\bm{\nabla}\\!\bm{\cdot}\\!\bm{T}^{\mathrm{\scriptscriptstyle{S}}}+w^{\mathrm{\scriptscriptstyle{S}}}\big{|}_{0}=\partial_{t}\partial_{z}\tfrac{1}{4}\overline{|\bm{\xi}_{1}|^{2}}\big{|}_{0}+\partial_{t}(\partial_{x}^{2}+\partial_{y}^{2})\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\tfrac{1}{4}|\bm{\xi}_{1}|^{2}\,\mathrm{d}z\,.$
(28)

With vertical integration of the horizontal components of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ over the wave-active zone

$\displaystyle\bm{\hat{x}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}$
$\displaystyle=\big{(}\overline{\zeta_{1}u_{1}}\big{)}\big{|}_{0}+\partial_{t}\partial_{x}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\tfrac{1}{2}\overline{\xi_{1}^{2}}\,\mathrm{d}z+\partial_{y}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\overline{\eta_{1}u_{1}}\,\mathrm{d}z\,,$
(29a)
$\displaystyle\bm{\hat{y}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}$
$\displaystyle=\big{(}\overline{\zeta_{1}v_{1}}\big{)}\big{|}_{0}+\partial_{x}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\overline{\xi_{1}v_{1}}\,\mathrm{d}z+\partial_{t}\partial_{y}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\tfrac{1}{2}\overline{\eta_{1}^{2}}\mathrm{d}z\,.$
(29b)

The horizontal divergence of $\bm{T}^{\mathrm{\scriptscriptstyle{S}}}$ is
therefore

$\bm{\nabla}\\!\bm{\cdot}\\!\bm{T}^{\mathrm{\scriptscriptstyle{S}}}=\partial_{x}\big{(}\overline{\zeta_{1}u_{1}}\big{)}\big{|}_{0}+\partial_{y}\big{(}\overline{\zeta_{1}v_{1}}\big{)}\big{|}_{0}+\partial_{t}\Bigg{(}\tfrac{1}{2}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\bm{\xi}_{1,\alpha}\bm{\xi}_{1,\beta}\,\mathrm{d}z\Bigg{)}_{,\alpha\beta}\,.$
(30)

In the final term the indices $\alpha$ and $\beta$ are $1$ and $2$.
Eliminating
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{T}^{\mathrm{\scriptscriptstyle{S}}}$ between
(28) and (30)

$\displaystyle w^{\mathrm{\scriptscriptstyle{S}}}\big{|}_{0}$
$\displaystyle=-\partial_{x}\big{(}\overline{\zeta_{1}u_{1}}\big{)}\big{|}_{0}-\partial_{y}\big{(}\overline{\zeta_{1}v_{1}}\big{)}\big{|}_{0}+\partial_{t}\partial_{z}\tfrac{1}{4}\overline{|\bm{\xi}_{1}|^{2}}\big{|}_{0}$
$\displaystyle\qquad+\partial_{t}(\partial_{x}^{2}+\partial_{y}^{2})\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\tfrac{1}{4}|\bm{\xi}_{1}|^{2}\,\mathrm{d}z-\partial_{t}\Bigg{(}\tfrac{1}{2}\int^{0}_{-\infty}\\!\\!\\!\\!\\!\\!\bm{\xi}_{1,\alpha}\bm{\xi}_{1,\beta}\,\mathrm{d}z\Bigg{)}_{,\alpha\beta}\,.$
(31)

The results in (28), (30) and (31) are all more complicated than their
solenoidal cousins in (22b), (23a) and (23b).

###  3.3 A comment on approximations to the Stokes transport

A widely used expression for the Stokes transport is

$\bm{T}^{\mathrm{\scriptscriptstyle{S}}}\approx\overline{(\zeta_{1}\bm{u}_{1})}\Big{|}_{0}\,.$
(32)

Expression (32) is exact for a uniform ($\mu=0$) progressive wave and is a
leading-order approximation for a slowly modulated ($\mu\ll 1$) wave packet.
To see the connection with the more general and exact results above, align the
$x$-axis with the horizontal wavevector so that $\eta_{1}=v_{1}=0$; then
$\bm{\hat{y}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}=\bm{\hat{y}}\bm{\cdot}\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}=\chi=0$.
Thus, with an error by order $\epsilon^{2}\mu$, and in agreement with (32),

$\tfrac{1}{2}\big{(}\overline{\zeta_{1}u_{1}}-\overline{w_{1}\xi_{1}}\big{)}\Big{|}_{0}\approx\overline{(\zeta_{1}u_{1})}\Big{|}_{0}\,.$
(33)

In other words, (32) is a leading-order approximation and at this order
$\bm{T}^{\mathrm{\scriptscriptstyle{S}}}$ and
$\bm{T}^{\mathrm{\scriptscriptstyle{S}}}_{\text{sol}}$ are identical.

##  4 glm

We return to the definition of the Stokes velocity and show how the glm theory
of Soward & Roberts [24] rationalizes and generalizes the heuristic
construction of the solenoidal velocity
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$. In §22.1 we
emphasized the ambiguity in the decomposition (6a) of trajectories into a mean
part $\bm{x}^{\mathrm{\scriptscriptstyle{L}}}$ and a perturbation $\bm{\xi}$.
GLM resolves this ambiguity by imposing that

$\overline{\bm{\xi}}=0\,.$ (34)

While (34) is widely accepted, it is neither inevitable nor particularly
natural. It implies that the Lagrangian-mean trajectory of a particle is
defined by the equality
$\bm{x}^{\mathrm{\scriptscriptstyle{L}}}=\overline{\bm{x}}$ between the
coordinates of the Lagrangian-mean position and the average of the coordinates
of the particle. This indicates that the Lagrangian-mean trajectory and, as a
result, the Lagrangian-mean velocity $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$
depend on a choice of coordinates. This undesirable feature of GLM is remedied
by glm, while also ensuring that $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ is
non-divergent (see [6] for other alternatives to GLM).

###  4.1 Formulation

To introduce glm, it is convenient to rewrite equation (4) governing the fluid
trajectories in terms of the flow map $\bm{\varphi}(\bm{x},t,\alpha,\epsilon)$
giving the position at time $t$ of the fluid particle initially at $\bm{x}$.
Equation (4) becomes

$\dot{\bm{\varphi}}(\bm{x},t,\alpha,\epsilon)=\bm{u}(\bm{\varphi}(\bm{x},t,\alpha,\epsilon),t,\alpha,\epsilon)\,.$
(35)

The decomposition of trajectories into mean and perturbation is best written
as the composition

$\bm{\varphi}=\bm{\Xi}\circ\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}\,,$
(36)

or, more explicitly,
$\bm{\varphi}(\bm{x},t,\alpha,\epsilon)=\bm{\Xi}(\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}(\bm{x},t,\epsilon),t,\alpha,\epsilon)$.
Here $\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}$ is the
($\alpha$-independent) mean map, sending the initial position of particles to
their mean position, and $\bm{\Xi}$ is the perturbation map, sending the mean
position to the exact, perturbed position. The perturbation map can only be
represented in the familiar form
$\bm{x}\mapsto\bm{x}+\bm{\xi}(\bm{x},t,\alpha,\epsilon)$, as in (6a), in
Euclidean space, where positions $\bm{x}$ can be identified with vectors and
added, or, in more general geometries, once a specific coordinate system has
been chosen and $\bm{x}$ is interpreted as a triple of coordinates. The
smallness of the perturbation, usually stated at $|\bm{\xi}|\ll 1$, translates
into the requirement that $\bm{\Xi}$ is close to the identity map. We
emphasise that the mean map ${\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}}$
is not obtained from $\bm{\varphi}$ by applying an averaging operator:
averaging is a linear operation that applies to linear objects, such as vector
fields, but not to nonlinear maps such as $\bm{\varphi}$. Instead,
$\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}$ is defined by imposing a
condition on the perturbation map $\bm{\Xi}$. The form of the Lagrangian-mean
velocity $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$, defined by

$\dot{\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}}(\bm{x},t,\epsilon)=\epsilon^{2}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\big{(}\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}(\bm{x},t),t,\epsilon\big{)}\,,\quad\textrm{i.e.}\quad\epsilon^{2}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\dot{\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}}\circ(\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}})^{-1}\,,$
(37)

depends on this condition.

The defining condition of glm is expressed as follows. The small parameter
$\epsilon$ is regarded as a fictitious time, and the perturbation map
$\bm{\Xi}$ is constructed as the flow at ‘time’ $\epsilon$ of a vector field,
$\bm{q}$ say; that is, $\bm{\Xi}(\bm{x},t,\epsilon)$ is the solution of

$\partial_{\epsilon}\bm{\Xi}(\bm{x},t,\alpha,\epsilon)=\bm{q}(\bm{\Xi}(\bm{x},t,\alpha,\epsilon),t,\alpha,\epsilon)\quad\textrm{with}\
\ \bm{\Xi}(\bm{x},t,\alpha,\epsilon=0)=\bm{x}\,,$ (38)

and $t$ is treated as a fixed parameter. The glm condition is then

$\overline{\bm{q}}=0.$ (39)

Although superficially similar to the GLM condition (34), (39) is
fundamentally different in that it is an intrinsic statement, applicable to
any manifold and independent of any coordinate choice. Moreover, the glm
formulation defines an exactly divergence-free Lagrangian mean flow, by
requiring that

$\bm{\nabla}\bm{\cdot}\bm{q}=0$ (40)

to ensure that $\bm{\Xi}$ and $\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}$
preserve volume and hence
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=0$.

Eq. (39) generalizes the condition on $\overline{\bm{\xi}_{2}}$ in (15) which
leads to $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ as the
Stokes velocity. We show this by solving (38) by Taylor expansion. In
coordinates, we have that

$\bm{\Xi}(\bm{x},t,\alpha,\epsilon)=\bm{x}+\bm{\xi}(\bm{x},t,\alpha,\epsilon)=\bm{x}+\epsilon\bm{q}(\bm{x},t,\alpha,0)+\tfrac{1}{2}\epsilon^{2}(\partial_{\epsilon}\bm{q}+\bm{q}\bm{\cdot}\bm{\nabla}\bm{q})(\bm{x},t,\alpha,0)+\cdots.$
(41)

The power series expansions (6b) for $\bm{\xi}$ and

$\bm{q}=\bm{q}_{1}+\epsilon\bm{q}_{2}+\epsilon^{2}\bm{q}_{3}+\cdots$ (42)

for $\bm{q}$ then give

$\bm{\xi}_{1}=\bm{q}_{1}\quad\textrm{and}\quad\bm{\xi}_{2}=\bm{q}_{2}+\tfrac{1}{2}\bm{q}_{1}\bm{\cdot}\bm{\nabla}\bm{q}_{1}$
(43)

so that (39) implies (15). In the next section, we develop a systematic
computation of $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ and hence
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ order by order in
$\epsilon$ using Lie series. This removes the need to introduce $\bm{\xi}$ by
focussing the perturbation expansion on $\bm{q}$.

###  4.2 Lie series expansion

The glm formalism can be regarded as an instance of classical perturbation
theory, which approximates solutions to the ordinary differential equation
(35) by performing a variable transformation designed to eliminate fast time
dependences. In this interpretation, $\bm{\Xi}$ determines the variable
transformation and $\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}$ represents
the new variable. Lie series [15, 21] provide a powerful tool for the
systematic implementation of classical perturbation theory which we now apply
to glm.

Introducing the decomposition (36) into (35) and using (37) gives

$\bm{w}+\bm{\Xi}_{*}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\bm{u},$ (44)

where

$\bm{w}=\dot{\bm{\Xi}}\circ\bm{\Xi}^{-1}$ (45)

is the perturbation velocity and $\bm{\Xi}_{*}$ is the push-forward by
$\bm{\Xi}$, with
$\bm{\Xi}_{*}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=(\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\bm{\cdot}\bm{\nabla})\bm{\Xi}$
in Cartesian coordinates. Pulling back (44) gives

$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\bm{\Xi}^{*}\bm{u}-\hat{\bm{w}},\quad\textrm{where}\
\ \hat{\bm{w}}=\bm{\Xi}^{*}\bm{w}.$ (46)

We seek an $\epsilon$-dependent $\bm{\Xi}$ to eliminate fast time dependence
from $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ order-by-order in $\epsilon$,
and formulate the problem in terms of the vector field $\bm{q}$ that generates
$\bm{\Xi}$ according to (38). We impose that $\bm{\nabla}\bm{\cdot}\bm{q}=0$
to ensure that
$\bm{\nabla}\bm{\cdot}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=0$, and the glm
condition (39).

Expanding $\bm{q}$ as in (42) we relate the various terms by differentiating
(46) repeatedly with respect to $\epsilon$ and evaluating the results at
$\epsilon=0$. Two key identities turns this into a mechanical exercise. The
first, essentially the definition of the Lie derivative [5], is

$\partial_{\epsilon}\bm{\Xi}^{*}\bm{u}=\bm{\Xi}^{*}\mathcal{L}_{\bm{q}}\bm{u},$
(47)

where
$\mathcal{L}_{\bm{q}}\bm{u}=\bm{q}\bm{\cdot}\bm{\nabla}\bm{u}-\bm{u}\bm{\cdot}\bm{\nabla}\bm{q}$
is the Lie derivative of $\bm{u}$ along $\bm{q}$. The second,

$\partial_{\epsilon}\hat{\bm{w}}=\partial_{t}(\bm{\Xi}^{*}\bm{q})+\mathcal{L}_{\bm{\Xi}^{*}\bm{q}}\hat{\bm{w}},$
(48)

is established in Appendix A. Iterating (47) and (48) we find

$\displaystyle\bm{\Xi}^{*}\bm{u}$
$\displaystyle=\bm{u}+\epsilon\mathcal{L}_{\bm{q}_{1}}\bm{u}+\tfrac{1}{2}\epsilon^{2}\left(\mathcal{L}^{2}_{\bm{q}_{1}}+\mathcal{L}_{\bm{q}_{2}}\right)\bm{u}+\tfrac{1}{6}\epsilon^{3}\left(\mathcal{L}^{3}_{\bm{q}_{1}}+2\mathcal{L}_{\bm{q}_{1}}\mathcal{L}_{\bm{q}_{2}}+\mathcal{L}_{\bm{q}_{2}}\mathcal{L}_{\bm{q}_{1}}+2\mathcal{L}_{\bm{q}_{3}}\right)\bm{u}+\cdots$
(49) $\displaystyle\hat{\bm{w}}$
$\displaystyle=\epsilon\partial_{t}\bm{q}_{1}+\tfrac{1}{2}\epsilon^{2}\left(\partial_{t}\bm{q}_{2}+\mathcal{L}_{\bm{q}_{1}}\partial_{t}\bm{q}_{1}\right)+\tfrac{1}{3}\epsilon^{3}\left(\partial_{t}\bm{q}_{3}+\mathcal{L}_{\bm{q}_{1}}\partial_{t}\bm{q}_{2}+\tfrac{1}{2}\mathcal{L}_{\bm{q}_{2}}\partial_{t}\bm{q}_{1}+\tfrac{1}{2}\mathcal{L}^{2}_{\bm{q}_{1}}\partial_{t}\bm{q}_{1}\right)+\cdots$
(50)

Introducing (49)–(50) into (46) gives at the first two orders in $\epsilon$,

$\partial_{t}\bm{q}_{1}=\bm{u}_{1},\quad\textrm{i.e.}\quad\bm{q}_{1}=\bm{\xi}_{1}$
(51)

and

$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\overline{\bm{u}}_{2}+\tfrac{1}{2}\overline{\mathcal{L}_{\bm{q}_{1}}\bm{u}_{1}}$
(52)

on using (39). This provides the geometric formula

$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}=\tfrac{1}{2}\overline{\mathcal{L}_{\bm{q}_{1}}\bm{u}_{1}}$
(53)

for the solenoidal Stokes drift, equivalent to (14a) since
$\mathcal{L}_{\bm{q}_{1}}\bm{u}_{1}=\bm{q}_{1}\bm{\cdot}\bm{\nabla}\bm{u}_{1}-\bm{u}_{1}\bm{\cdot}\bm{\nabla}\bm{q}_{1}=\bm{\nabla}\times(\bm{u}_{1}\times\bm{\xi}_{1})$.
The expansion can be pursued to higher orders by choosing divergence-free
$\bm{q}_{n}$ that push the fast dependence on the right-hand side of (46) to
order $\epsilon^{n+1}$. This yields $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$
and hence $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ to
arbitrary order in $\epsilon$.

###  4.3 Lagrangian-mean momentum equation

A Lagrangian-mean momentum equation, analogous to Andrews & McIntyre Theorem I
[1], can be derived for the glm formalism order-by-order in $\epsilon$ [24,
6]. We sketch a derivation, focussing on the case of small-amplitude surface
gravity waves. This derivation is conveniently carried out using a
representation of the rotating Euler equation

$D_{t}\bm{u}+\bm{f}\times\bm{u}=-\bm{\nabla}p,$ (54)

with $D_{t}=\partial_{t}+\bm{u}\bm{\cdot}\bm{\nabla}$, in terms of the
absolute momentum $1$-form [6, 11, 12]

$\nu_{a}=\bm{u}\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}.$
(55)

It can be checked that (54) is equivalent to

$(\partial_{t}+\mathcal{L}_{\bm{u}})\nu_{a}=-\mathrm{d}\pi,$ (56)

where
$\pi=p-\tfrac{1}{2}|\bm{u}|^{2}-\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\bm{u}$,
using basic properties of the Lie derivative, namely Leibniz rule, commutation
with the differential $\mathrm{d}$, and that
$\mathcal{L}_{\bm{u}}=\bm{u}\bm{\cdot}\nabla$ when applied to scalars [5]; see
Appendix A for details. Equation (56) can be thought of as a local version of
Kelvin’s circulation theorem. This is readily obtained by integration along a
closed curve $\mathcal{C}(t)$ moving with $\bm{u}$ to find

$\oint_{\mathcal{C}(t)}\nu_{a}=\oint_{\mathcal{C}(t)}\left(\bm{u}+\tfrac{1}{2}(\bm{f}\times\bm{x})\right)\bm{\cdot}\mathrm{d}\bm{x}=\mathrm{const}.$
(57)

An advantage of (56) is that it directly leads to a Lagrangian-mean momentum
equation of a similar form [6],

$(\partial_{t}+\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}})\nu_{\mathrm{a}}^{\mathrm{\scriptscriptstyle{L}}}=-\mathrm{d}\pi^{\mathrm{\scriptscriptstyle{L}}}\,,$
(58)

and to the corresponding Lagrangian-mean Kelvin’s circulation theorem

$\oint_{\mathcal{C}^{\mathrm{\scriptscriptstyle{L}}}(t)}\nu_{\mathrm{a}}^{\mathrm{\scriptscriptstyle{L}}}=\mathrm{const},$
(59)

where the closed curve $\mathcal{C}^{\mathrm{\scriptscriptstyle{L}}}(t)$ moves
with the Lagrangian-mean velocity $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$.
Here the Lagrangian-mean momentum and effective pressure are given by

$\nu_{\mathrm{a}}^{\mathrm{\scriptscriptstyle{L}}}=\overline{\bm{\Xi}^{*}\nu_{a}}\quad\textrm{and}\quad\pi^{\mathrm{\scriptscriptstyle{L}}}=\overline{\bm{\Xi}^{*}\pi}\,.$
(60)

The pull-back $\bm{\Xi}^{*}$ by the perturbation map acts on scalars as a
composition, e.g.
$(\bm{\Xi}^{*}\pi)(\bm{x},t,\alpha,\epsilon)=\pi(\bm{\Xi}(\bm{x},t,\alpha,\epsilon),t,\alpha,\epsilon)$,
and commutes with the differential so that

$\left(\bm{\Xi}^{*}(\bm{u}\bm{\cdot}\mathrm{d}\bm{x})\right)(\bm{x})=\bm{u}(\bm{\Xi}(\bm{x}))\bm{\cdot}\mathrm{d}\bm{\Xi}(\bm{x})=\bm{u}(\bm{\Xi}(\bm{x}))\bm{\cdot}(\mathrm{d}\bm{x}\bm{\cdot}\bm{\nabla})\bm{\Xi}(\bm{x}).$
(61)

We show in Appendix A that (58), together with the coordinate representation
$\bm{\Xi}(\bm{x},t,\alpha,\epsilon)=\bm{x}+\bm{\xi}(\bm{x},t,\alpha,\epsilon)$
and GLM condition (34) recovers Andrews & McIntyre’s Lagrangian-mean momentum
equation [1, Theorem I]. For glm, we can use the Lie-series expression (49) to
obtain an expansion of $\nu_{\mathrm{a}}^{\mathrm{\scriptscriptstyle{L}}}$ as

$\nu_{\mathrm{a}}^{\mathrm{\scriptscriptstyle{L}}}=\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}+\epsilon^{2}\,\overline{\bm{u}_{2}\bm{\cdot}\mathrm{d}\bm{x}+\mathcal{L}_{\bm{q}_{1}}(\bm{u}_{1}\bm{\cdot}\mathrm{d}\bm{x})+\tfrac{1}{4}\mathcal{L}_{\bm{q}_{1}}^{2}\left((\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\right)}+O(\epsilon^{3})$
(62)

using that $\bm{u}=\epsilon\bm{u}_{1}+\epsilon^{2}\bm{u}_{2}+\cdots$ and
$\overline{\bm{q}_{1}}=\overline{\bm{q}_{2}}=0$. Introducing (62) into (58)
and expanding the Lie derivative gives an evolution equation for the Eulerian
mean velocity $\bm{u}^{\mathrm{\scriptscriptstyle{E}}}=\overline{\bm{u}_{2}}$
supplemented by the incompressibility condition
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{E}}}=0$.

The Lagrangian-mean momentum equation greatly simplifies for surface gravity
waves: typical wave frequencies $\sigma$ satisfy $\sigma/f\gg 1$ and we can
therefore neglect the term involving $\mathcal{L}_{\bm{q}_{1}}^{2}$ in (62)
against the other two in the average (recall that
$\partial_{t}\bm{q}_{1}=\bm{u}_{1}$). Moreover, using the irrotationality
condition (16), we compute

$\displaystyle\overline{\mathcal{L}_{\bm{q}_{1}}(\bm{u}_{1}\bm{\cdot}\mathrm{d}\bm{x})}$
$\displaystyle=\overline{(\bm{q}_{1}\bm{\cdot}\bm{\nabla})\bm{u}_{1}\bm{\cdot}\mathrm{d}\bm{x}+\bm{u}_{1}\bm{\cdot}\mathrm{d}\bm{q}_{1}}=\overline{(q_{1j}u_{1i,j}+u_{1j}q_{1j,i})}\mathrm{d}x_{i}$
$\displaystyle=\overline{(q_{1j}u_{1j,i}+u_{1j}q_{1j,i})}\mathrm{d}x_{i}=\partial_{t}\overline{(q_{1j}q_{1j,i})}\mathrm{d}x_{i}$
$\displaystyle=\tfrac{1}{2}\partial_{t}\bm{\nabla}\overline{|\bm{q}_{1}|^{2}}\cdot\mathrm{d}\bm{x}=\tfrac{1}{2}\mathrm{d}(\partial_{t}\overline{|\bm{q}_{1}|^{2}})=O(\mu\epsilon^{2}).$
(63)

Therefore, to leading order, the glm mean momentum equation reduces to

$(\partial_{t}+\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}})\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\right)=-\mathrm{d}\pi^{\mathrm{\scriptscriptstyle{L}}},$
(64)

with the glm Lagrangian-mean velocity in (15) (see [11, 12] for an analogous
formulation of the GLM CL equation). Using Cartan’s formula in the form

$\mathcal{L}_{\bm{u}}(\bm{v}\bm{\cdot}\mathrm{d}\bm{x})=\left((\bm{\nabla}\\!\times\\!\bm{v})\times\bm{u}\right)\bm{\cdot}\mathrm{d}\bm{x}+\mathrm{d}(\bm{u}\cdot\bm{v})\,,$
(65)

the Lie derivatives in (64) can be written as

$\displaystyle\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}}\left((\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\right)$
$\displaystyle=2\left(\bm{f}\times\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\right)\bm{\cdot}\mathrm{d}\bm{x}+\mathrm{d}(\cdot),$
(66)
$\displaystyle\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}}\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\mathrm{d}\bm{x}\right)$
$\displaystyle=\mathcal{L}_{\bm{u}^{\mathrm{E}}}\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\mathrm{d}\bm{x}\right)+\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}}\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\mathrm{d}\bm{x}\right)$
$\displaystyle=(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\bm{\nabla}\bm{u}^{\mathrm{\scriptscriptstyle{E}}})\bm{\cdot}\mathrm{d}\bm{x}+\left((\bm{\nabla}\\!\times\\!\bm{u}^{\mathrm{\scriptscriptstyle{E}}})\times\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}\right)\bm{\cdot}\mathrm{d}\bm{x}+\mathrm{d}(\cdot),$
(67)

where we do not detail the exact differentials. This makes it possible to
rewrite (64) as

$\partial_{t}\bm{u}^{\mathrm{\scriptscriptstyle{E}}}+\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\\!\bm{\cdot}\\!\bm{\nabla}\right)\,\bm{u}^{\mathrm{\scriptscriptstyle{E}}}+\bm{f}\times(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}+\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}})=-\bm{\nabla}\varpi+\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}\times(\bm{\nabla}\\!\times\\!\bm{u}^{\mathrm{\scriptscriptstyle{E}}})\,,$
(68)

for a suitable definition of the effective pressure $\varpi$. Equation (68)
can be recognized as the CL equation [4, 14, 10] with the solenoidal Stokes
velocity $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ replacing
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$. This is not surprising since the
$O(\epsilon^{2}\mu)$ difference between
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ and
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ is of the same order as terms
neglected in the derivation of the CL equation. (CL geared $\mu$ to $\epsilon$
by taking $\mu=\epsilon^{2}$. In our derivation this gearing is not
necessary.) In fact the assumption $\mu\ll 1$ is only used to neglect the term
(63) from the Lagrangian-mean momentum equation to obtain (64) and hence (68).
Restoring this term leads to the generalisation

$(\partial_{t}+\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}})\left(\bm{u}^{\mathrm{\scriptscriptstyle{E}}}\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}\mathrm{d}(\partial_{t}\overline{|\bm{q}_{1}|^{2}})+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\right)=-\mathrm{d}\pi^{\mathrm{\scriptscriptstyle{L}}}\,,$
(69)

of the CL equation valid for $\mu=O(1)$. Since
$\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}}$ and $\mathrm{d}$
commute,
$\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}}\mathrm{d}(\partial_{t}\overline{|\bm{q}_{1}|^{2}})=\mathrm{d}(\cdots)$
and the additional terms involving $|\bm{q}_{1}|^{2}$ can be absorbed in the
differential of the pressure-like term on the right-hand side. We conclude
that the CL equation (68) with
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ as Stokes velocity
holds for $\mu=O(1)$ provided that the effective pressure $\varpi$ is suitably
redefined.

##  5 Conclusion

This paper examines a problematic aspect of the Stokes velocity
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$, namely its divergence
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and non-
zero vertical component $w^{\mathrm{\scriptscriptstyle{S}}}$. Some confusion
has arisen because
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and
$w^{\mathrm{\scriptscriptstyle{S}}}$ are small when the wave field has an
amplitude that varies on scales longer and slower than the wavelength and
period. This scale-separation approximation corresponds to the existence of a
small parameter $\mu\ll 1$, as is the case for slowly varying wavepackets. The
distinction between approximate and exact results in the literature is not
always made plain. Beyond this, there are no irretrievable difficulties with
the familiar form (2) of $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$: when
$\bm{\nabla}\\!\bm{\cdot}\\!\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$ and
$w^{\mathrm{\scriptscriptstyle{S}}}$ cannot be neglected, they are readily
computed in terms of first-order fields.

The main point of the paper, however, is that the Stokes velocity, understood
as the difference
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}=\bm{u}^{\mathrm{\scriptscriptstyle{L}}}-\bm{u}^{\mathrm{\scriptscriptstyle{E}}}$
between Lagrangian- and Eulerian-mean velocities, is not uniquely defined and
that an alternative version,
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ in (14), that is
exactly solenoidal can serve as a convenient substitute for the familiar (2).
The non-uniqueness arises because there is no single definition of the
Lagrangian-mean velocity $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$, which is
only constrained to serve as a good ‘guiding-centre’ representative of the
motion of an ensemble of fluid particles. The ambiguity is usually resolved by
imposing the GLM condition $\overline{\bm{\xi}}=0$, which is restricted to
Euclidean geometry or depends on coordinates, but alternatives which have the
benefit of both coordinate-independence and non-divergent
$\bm{u}^{\mathrm{\scriptscriptstyle{L}}}$ exist [6]. One of these, Soward &
Roberts’ glm [24], leads to
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ in (14) or, in a more
geometric form, (53). We give only the form of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ to leading order in
the wave amplitude parameter $\epsilon$; as we describe, higher-order
corrections can be computed systematically using a Lie series expansion. We
further show that the Craik–Leibovich (CL) equations [4, 14], widely used to
model the feedback of surface gravity waves on the flow, apply in the glm
framework, with $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$
simply replacing $\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$.

We stress that, in general, the coordinates associated with a Lagrangian-mean
trajectory
$\bm{x}^{\mathrm{\scriptscriptstyle{L}}}=\bm{\varphi}^{\mathrm{\scriptscriptstyle{L}}}(\bm{x}_{0},t,\epsilon)$
of given fluid particle (identified by initial position $\bm{x}_{0}$) are not
the average of the coordinates of this particle. This is a feature specific to
GLM, intimately connected to its coordinate dependence and the decision to
make $\overline{\bm{\xi}}=0$.

We conclude by noting that the assumption $\mu\ll 1$ which makes it possible
to ignore the divergence and vertical component of
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}$, is restrictive. The assumption
clearly holds for single wavepackets, but may fail for more complex wave
fields. When it does, using the solenoidal Stokes velocity
$\bm{u}^{\mathrm{\scriptscriptstyle{S}}}_{\textrm{sol}}$ brings considerable
simplifications, e.g. for the computation of the Eulerian mean flow from the
CL or other wave-averaged equations.

Acknowledgments.JV is supported by the UK Natural Environment Research
Council, grant NE/W002876/1. WRY is supported by National Science Foundation
Award OCE-2048583. We thank Oliver Bühler and an anonymous referee for
comments that improved this work.

##  Appendix A Computational details for §4

###  A.1 Derivation of (48)

We define $\hat{\bm{q}}=\bm{\Xi}^{*}\bm{q}$ in analogy with
$\hat{\bm{w}}=\bm{\Xi}^{*}\bm{w}$ and note that both are the negative of
velocities associated with the inverse map $\bm{\Xi}^{-1}$ in the sense that

$\partial_{t}\bm{\Xi}^{-1}(\bm{x},t,\epsilon)=-\hat{\bm{w}}(\bm{\Xi}(\bm{x},t,\epsilon),t,\epsilon)\quad\textrm{and}\quad\partial_{\epsilon}\bm{\Xi}^{-1}(\bm{x},t,\epsilon)=-\hat{\bm{q}}(\bm{\Xi}(\bm{x},t,\epsilon),t,\epsilon),$
(70)

as follows from the differentiation of the identity
$\bm{\Xi}^{-1}(\bm{\Xi}(\bm{x},t,\epsilon),t,\epsilon)=\bm{x}$ with respect to
$t$ and $\epsilon$. (We do not make the dependence on $\alpha$ explicit for
simplicity.) Equating the derivative of the first equality with respect to
$\epsilon$ with the derivative of the second with respect to $t$ gives

$-\partial_{\epsilon}\hat{\bm{w}}+\hat{\bm{q}}\bm{\cdot}\bm{\nabla}\hat{\bm{w}}=-\partial_{t}\hat{\bm{q}}+\hat{\bm{w}}\bm{\cdot}\bm{\nabla}\hat{\bm{q}}$
(71)

and, on rearranging,

$\partial_{\epsilon}\hat{\bm{w}}=\partial_{t}\hat{\bm{q}}+\hat{\bm{q}}\bm{\cdot}\bm{\nabla}\hat{\bm{w}}-\hat{\bm{w}}\bm{\cdot}\bm{\nabla}\bm{q}=\partial_{t}\hat{\bm{q}}+\mathcal{L}_{\hat{\bm{q}}}\hat{\bm{w}},$
(72)

i.e. (48).

###  A.2 Equivalence between (54) and (56)

We show that (56) is equivalent to the familiar form (56) of the Euler
equation. We first write (56) explicitly as

$(\partial_{t}+\mathcal{L}_{\bm{u}})\left(\bm{u}\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\right)=-\mathrm{d}\left(p-\tfrac{1}{2}|\bm{u}|^{2}-\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\bm{u}\right).$
(73)

Expanding, we find

$\displaystyle(D_{t}\bm{u})\bm{\cdot}\mathrm{d}\bm{x}+\bm{u}\bm{\cdot}\mathrm{d}\bm{u}+\tfrac{1}{2}(\bm{f}\times\bm{u})\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{u}$
$\displaystyle=-\nabla
p\bm{\cdot}\mathrm{d}\bm{x}+\bm{u}\bm{\cdot}\mathrm{d}\bm{u}+\tfrac{1}{2}(\bm{f}\times\mathrm{d}\bm{x})\bm{\cdot}\bm{u}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{u},$
(74)

which recovers (54) since
$(\bm{f}\times\mathrm{d}\bm{x})\bm{\cdot}\bm{u}=-(\bm{f}\times\bm{u})\bm{\cdot}\mathrm{d}\bm{x}$.

###  A.3 GLM from (58)

We show that (58), together with the defining property of GLM
$\bar{\bm{\xi}}=0$ is equivalent to Theorem I of Andrews & McIntyre [1] in the
case of an incompressible fluid. We need to apply the pull-back $\bm{\Xi}^{*}$
to

$\nu_{a}=\bm{u}\bm{\cdot}\mathrm{d}\bm{x}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}\quad\textrm{and}\quad\pi=p-\tfrac{1}{2}|\bm{u}|^{2}-\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\bm{u},$
(75)

then average. Here $\bm{u}$ stands for the 3 components of $\bm{u}$ rather
than for a vector, hence the pull-back is simply composition with $\bm{\Xi}$,
that is, $\bm{\Xi}^{*}\bm{u}=\bm{u}\circ\bm{\Xi}=:\bm{u}^{\xi}$. We also have
$\bm{\Xi}^{*}\bm{x}=\bm{x}+\bm{\xi}(\bm{x},t,\epsilon)$ and
$\bm{\Xi}^{*}\mathrm{d}\bm{x}=\mathrm{d}(\bm{\Xi}^{*}\bm{x})=\mathrm{d}\bm{x}+\mathrm{d}\bm{\xi}$.
Using that $\bar{\bm{\xi}}=0$, we compute

$\nu_{a}^{\mathrm{L}}=\overline{\bm{u}^{\xi}\bm{\cdot}\mathrm{d}\bm{x}+\bm{u}^{\xi}\bm{\cdot}\mathrm{d}\bm{\xi}+\tfrac{1}{2}(\bm{f}\times\bm{\xi})\bm{\cdot}\mathrm{d}\bm{\xi}}+\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\mathrm{d}\bm{x}.$
(76)

The averaged term can be written as
$(\bm{u}^{\mathrm{\scriptscriptstyle{L}}}-\bm{\mathsf{p}})\bm{\cdot}\mathrm{d}\bm{x}$
with $\bm{u}^{\mathrm{\scriptscriptstyle{L}}}=\overline{\bm{u}^{\xi}}$ and

$\mathsf{p}_{i}=-\overline{\xi_{j,i}(u^{\xi}_{j}+\tfrac{1}{2}\left(\bm{f}\times\bm{\xi})_{j}\right)}$
(77)

the pseudomomentum. Similarly,

$\pi^{\mathrm{\scriptscriptstyle{L}}}=\overline{{p}^{\xi}-\tfrac{1}{2}|\bm{u}^{\xi}|^{2}-\tfrac{1}{2}\left(\bm{f}\times\bm{\xi}\right)\bm{\cdot}\bm{u}^{\xi}}-\tfrac{1}{2}\left(\bm{f}\times\bm{x}\right)\bm{\cdot}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}.$
(78)

Introducing (76) and (78) into (58), we compute

$\displaystyle(\partial_{t}+\mathcal{L}_{\bm{u}^{\mathrm{\scriptscriptstyle{L}}}})\left((\bm{u}^{\mathrm{\scriptscriptstyle{L}}}-\bm{\mathsf{p}})\bm{\cdot}\,\mathrm{d}\bm{x}\right)$
$\displaystyle+\left(\left(\bm{\nabla}\\!\times\\!\tfrac{1}{2}(\bm{f}\times\bm{x})\right)\times\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\right)\bm{\cdot}\mathrm{d}\bm{x}+\mathrm{d}\left(\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\right)$
$\displaystyle=-\nabla\tilde{\pi}^{\mathrm{\scriptscriptstyle{L}}}\bm{\cdot}\mathrm{d}\bm{x}+\mathrm{d}\left(\tfrac{1}{2}(\bm{f}\times\bm{x})\bm{\cdot}\bm{u}^{\mathrm{\scriptscriptstyle{L}}}\right),$
(79)

where
$\tilde{\pi}^{\mathrm{\scriptscriptstyle{L}}}=\overline{{p}^{\xi}-\tfrac{1}{2}|\bm{u}^{\xi}|^{2}-\tfrac{1}{2}(\bm{f}\times\bm{\xi})\cdot\bm{u}^{\xi}}$
and we have used Cartan’s formula given below (64). Since
$\bm{\nabla}\\!\times\\!(\bm{f}\times\bm{x})=2\bm{f}$ we finally obtain

$D_{t}(u^{\mathrm{L}}_{i}-\mathsf{p}_{i})+(u^{\mathrm{L}}_{j}-\mathsf{p}_{j})u^{\mathrm{L}}_{j,i}+(\bm{f}\times\bm{u}^{\mathrm{L}})_{i}=-\tilde{\pi}^{\mathrm{L}}_{,i},$
(80)

that is, Theorem I of [1].

## References

  * [1] D. G. Andrews and M. E. McIntyre.  An exact theory of nonlinear waves on a Lagrangian-mean flow.  J. Fluid Mech., 89:609–646, 1978.
  * [2] O. Bühler.  Waves and Mean Flows.  Cambridge University Press, 2nd edition, 2014.
  * [3] Oliver Bühler and Miranda Holmes-Cerfon.  Particle dispersion by random waves in rotating shallow water.  Journal of Fluid Mechanics, 638:5–26, 2009.
  * [4] A. D. Craik and S. Leibovich.  A rational model for langmuir circulations.  J. Fluid Mech., 73:401–426, 1976.
  * [5] T. Frankel.  The Geometry of Physics.  Cambridge University Press, 2011.
  * [6] A. D. Gilbert and J. Vanneste.  Geometric generalised Lagrangian-mean theories.  J. Fluid Mech., 839:95–134, 2018.
  * [7] S. Haney and William R. Young.  Radiation of internal waves from groups of surface gravity waves.  J. Fluid Mech., 829:280–303, 2017.
  * [8] K. Hasselmann.  Wave‐driven inertial oscillations.  Geophys. Astrophys. Fluid Dynam., 1:463–502, 1970.
  * [9] C. Higgins, T. S. van den Bremer, and J. Vanneste.  Lagrangian transport by deep-water surface gravity wavepackets: effects of directional spreading and stratification.  J. Fluid Mech., 883:A42, 2020.
  * [10] D. D. Holm.  The ideal Craik–Leibovich equations.  Physica D, 98:415–449, 1996.
  * [11] D. D. Holm.  Stochastic closures for wave–current interaction dynamics.  J. Nonlinear Sci., 29:2987–3031, 2019.
  * [12] D. D. Holm.  Stochastic variational formulations of fluid wave–current interaction.  J. Nonlinear Sci., 31:4, 2021.
  * [13] N. E. Huang.  On surface drift currents in the ocean.  J. Fluid Mech., 91:191–208, 1979.
  * [14] S. Leibovich.  On wave-current interaction theories of langmuir circulations.  J. Fluid Mech., 99:715–724, 1980.
  * [15] A.J. Lichtenberg and M.A. Lieberman.  Regular and Stochastic Motion.  Springer–Verlag, 2nd edition, 1992.
  * [16] M. E. McIntyre.  A note on the divergence effect and the Lagrangian-mean surface elevation in periodic water waves.  J. Fluid Mech., 189:235–242, 1988.
  * [17] J. C. McWilliams and J. M. Restrepo.  The wave-driven ocean circulation.  J. Phys. Oceanogr., 29:2523–2540, 1999.
  * [18] J. C. McWilliams, J. M. Restrepo, and E. M. Lane.  An asymptotic theory for the interaction of waves and currents in coastal waters.  J. Fluid Mech., 511:135–178, 2004.
  * [19] G. Mellor.  On theories dealing with the interaction of surface waves and ocean circulation.  J. Geophys. Res. Oceans, 121:4474–4486, 2016.
  * [20] S. G. Monismith.  Stokes drift: theory and experiments.  J. Fluid Mech., 884:F1, 2019.
  * [21] A. H. Nayfeh.  Perturbation Methods.  Wiley VCH, 2000.
  * [22] O. M. Phillips.  The Dynamics of the Upper Ocean.  Cambridge University Press, 1977.
  * [23] Jerome A Smith.  Observed variability of ocean wave Stokes drift, and the Eulerian response to passing groups.  J. Phys. Oceanogr., 36(7):1381–1402, 2006.
  * [24] A. M. Soward and P. H. Roberts.  The hybrid euler–lagrange procedure using an extension of Moffatt’s method.  J. Fluid Mech., 661:45–72, 2010.
  * [25] G. G. Stokes.  On the theory of oscillatory waves.  Trans. Camb. Philos. Soc., 8:441–455, 1847.
  * [26] T. S. van den Bremer and Ø Breivik.  Stokes drift.  Phil. Trans. R. Soc. A, 376:20170104, 2018.
  * [27] T. S. van den Bremer, C. Whittaker, R. Calvert, A. Raby, and P. H. Taylor.  Experimental study of particle trajectories below deep-water surface gravity wave groups.  J. Fluid Mech., 879:168–186, 2019.
  * [28] Gregory L Wagner, Gregory P Chini, Ali Ramadhan, Basile Gallet, and Raffaele Ferrari.  Near-inertial waves and turbulence driven by the growth of swell.  J. Phys. Oceanogr., 51(5):1337–1351, 2021.