2. The equations of motion

In this section we derive the equations governing a mixture of widely separated vortex patches and surface-wave packets. In the wide-separation limit, the vortex patches correspond to point vortices and the wave packets correspond to ‘point dipoles’. We obtain our equations by coupling the Lagrangian for the wave field in the form proposed by Whitham (Reference Whitham1965, Reference Whitham1974) to the Lagrangian for two-dimensional incompressible flow representing the surface current. We use the Doppler term in the dispersion relation to couple the two Lagrangians together. This appears to be a simple and powerful method for deriving equations governing the interactions between waves and mean flows. Further details of the method are given by Salmon (Reference Salmon2020).

For the waves by themselves, the Lagrangian proposed by Whitham is

(2.1)\begin{equation} L_w[\theta, \mathcal{A}]=\iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x}\,( \omega-\omega_r (\boldsymbol{k}) -\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{k}) \mathcal{A}, \end{equation}

where the integral is over time and the ocean surface; the frequency $\omega = -\theta _t$ and wave vector $\boldsymbol {k}=\boldsymbol {\nabla } \theta$ are abbreviations for the derivatives of the wave phase $\theta (\boldsymbol {x},t)$;

(2.2)\begin{equation} \mathcal{A}=\frac{E}{\omega_r} \end{equation}

is the wave action; $E$ is the wave energy per unit area; and $\omega _r (\boldsymbol {k})$ is the prescribed relative frequency of the waves – the frequency measured in a reference frame moving at the mean flow velocity $\boldsymbol {U}(\boldsymbol {x},t)$. Our notation is $\boldsymbol {k}=(k,l)$, $\boldsymbol {x}=(x,y)$ and $\boldsymbol {\nabla }=(\partial _x,\partial _y)$. For surface waves,

(2.3)\begin{equation} \omega_r (\boldsymbol{k})= \sqrt{g\vert\boldsymbol{k}\vert}. \end{equation}

At this stage, we consider the mean flow to be prescribed. For the sake of completeness, Appendix A provides a systematic derivation of (2.1) following Whitham's averaged-Lagrangian method. Variations of $\mathcal {A}$ yield the dispersion relation:

(2.4)\begin{equation} \omega = \sqrt{g\vert\boldsymbol{k}\vert} +\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{k}. \end{equation}

From variations of $\theta$ we obtain

(2.5)\begin{align} \delta L_w[\theta, \mathcal{A}]&=\iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x} \left( -(\delta \theta)_t -\frac{\partial \omega_r}{\partial\boldsymbol{k}}\boldsymbol{\cdot} \boldsymbol{\nabla} (\delta \theta) -\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} (\delta \theta) \right) \mathcal{A} \nonumber\\ &=\iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x}\,( \mathcal{A}_t +\boldsymbol{\nabla} \boldsymbol{\cdot} [(\boldsymbol{c}_g+\boldsymbol{U}) \mathcal{A}] ) \delta \theta, \end{align}

where $\boldsymbol {c}_g(\boldsymbol {k})=\partial \omega _r / \partial \boldsymbol {k}$ is the relative group velocity. Thus we obtain the action conservation equation:

(2.6)\begin{equation} \mathcal{A}_t +\boldsymbol{\nabla} \boldsymbol{\cdot} [(\boldsymbol{c}_g+\boldsymbol{U}) \mathcal{A}] = 0. \end{equation}

We emphasize that the Lagrangian (2.1) applies to any truly two-dimensional system with arbitrary dispersion relation $\omega _r(\boldsymbol {k})$. However, the mixed layer of the ocean is not a two-dimensional system. To obtain a closed, two-dimensional model of the waves, one must do more than to assert (2.3); one must isolate the waves from their deeper surroundings. The most straightforward way to do this is to assume the existence of a rigid lower mixed-layer boundary at depth $H_0$. The need for such a strong assumption, and its implications, becomes more obvious as we proceed.

The Lagrangian for two-dimensional incompressible flow is

(2.7)\begin{equation} L_{m}[\alpha,\beta,\psi]=\iiint \textrm{d}t \,\textrm{d}\kern0.07em \boldsymbol{x}\, H_0 \left( -\alpha \beta_t +\psi \frac{\partial(\alpha,\beta)}{\partial(x,y)} +\frac{1}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \right), \end{equation}

where the subscript $m$ stands for ‘mean flow’ and

(2.8)\begin{equation} \frac{\partial(A,B)}{\partial(x,y)} \equiv [A,B] = A_x B_y - B_x A_y \end{equation}

is the Jacobian, defined for any two functions $A(x,y)$ and $B(x,y)$. The variables $\alpha (\boldsymbol {x},t)$, $\beta (\boldsymbol {x},t)$ and $\psi (\boldsymbol {x},t)$ represent averages over the constant mixed-layer depth $H_0$. We assume that the mean flow is depth-independent in this range. Stationarity of $L_{m}$ implies

(2.9)\begin{gather} \delta \alpha: \ \beta_t + [\psi,\beta]=0, \end{gather}

(2.10)\begin{gather}\delta \beta: \ \alpha_t + [\psi,\alpha]=0, \end{gather}

(2.11)\begin{gather}\delta \psi: \ [\alpha,\beta]= \nabla^2 \psi. \end{gather}

We see that $\alpha$ and $\beta$ are vorticity labels in the following sense: first, by (2.9) and (2.10), they are conserved following the fluid motion, and hence each fluid particle is identified by its two labels $(\alpha ,\beta )$; second, by (2.11), the vorticity in an arbitrary area of the flow is given by

(2.12)\begin{equation} \iint \textrm{d}\kern0.07em x\,\textrm{d}y \,\nabla^2 \psi = \iint \textrm{d}\kern0.07em x\,\textrm{d}y\, \frac{\partial(\alpha,\beta)}{\partial(x,y)} =\iint \textrm{d}\alpha \,\textrm{d}\beta, \end{equation}

where the integration is over the arbitrary area in physical space and the corresponding area in label space. Taking the time derivative of (2.11) and using the Jacobi identity,

(2.13)\begin{equation} [A,[B,C]]+[B,[C,A]]+[C,[A,B]]=0, \end{equation}

we obtain the vorticity equation,

(2.14)\begin{equation} \nabla^2 \psi _t+[\psi,\nabla^2 \psi]=0, \end{equation}

for the mean flow by itself.

We couple $L_w$ to $L_{m}$ by replacing the mean velocity $\boldsymbol {U}$ in (2.1) with $\boldsymbol {u}_\psi \equiv (-\psi _y,\psi _x)$, and by assuming that the Lagrangian for the entire system is the sum

(2.15)\begin{equation} L[\theta,\mathcal{A},\alpha,\beta,\psi]=L_w+L_{m}=\iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x}\, ( -\theta_t\mathcal{A} -H_0 \alpha \beta_t ) -\int \textrm{d}t \, H \end{equation}

of (2.1) and (2.7), where

(2.16)\begin{equation} H[\theta,\mathcal{A},\alpha,\beta,\psi] =\iint \textrm{d}\kern0.07em \boldsymbol{x} \left( \omega_r \mathcal{A} -H_0 \psi \frac{\partial(\alpha,\beta)}{\partial(x,y)} -\frac{H_0}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \psi +\boldsymbol{\nabla} \psi \times \boldsymbol{k} \mathcal{A} \right) \end{equation}

is the Hamiltonian, and $(\theta ,\mathcal {A})$ and $(\alpha ,\beta )$ are canonical pairs. Using (2.2), (2.11) and integrations by parts to evaluate (2.16) we find that

(2.17)\begin{equation} H=\iint \textrm{d}\kern0.07em \boldsymbol{x} \left( E +\frac{H_0}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \right). \end{equation}

Thus, as expected, our dynamics conserves the sum of the wave energy and the kinetic energy of the mean flow. The equations corresponding to $\delta L=0$ are

(2.18)\begin{gather} \delta \mathcal{A}: \ \omega \equiv{-} \theta_t =\omega_r(k,l)+\boldsymbol{u}_\psi\boldsymbol{\cdot} \boldsymbol{k}, \end{gather}

(2.19)\begin{gather}\delta \theta: \ \mathcal{A}_t +\boldsymbol{\nabla} \boldsymbol{\cdot} [(\boldsymbol{c}_g+\boldsymbol{u}_\psi) \mathcal{A}] = 0, \end{gather}

(2.20)\begin{gather}\delta \psi: \ H_0[\alpha,\beta]=H_0\nabla^2 \psi- \boldsymbol{\nabla} \times (\boldsymbol{k} \mathcal{A}) , \end{gather}

(2.21)\begin{gather}\delta \alpha: \ \beta_t + [\psi,\beta]=0, \end{gather}

(2.22)\begin{gather}\delta \beta: \ \alpha_t + [\psi,\alpha]=0, \end{gather}

where $\boldsymbol {\nabla } \times (A,B) \equiv B_x-A_y$ will be our notation for the vertical component of the curl of a horizontal vector. By the Jacobian identity, (2.20)–(2.22) imply

(2.23)\begin{equation} q_t + [\psi,q]=0, \end{equation}

where

(2.24)\begin{equation} q = H_0\nabla^2 \psi-\boldsymbol{\nabla} \times \boldsymbol{p}, \end{equation}

and

(2.25)\begin{equation} \boldsymbol{p} = \boldsymbol{k} \mathcal{A} \end{equation}

is the pseudomomentum. The wave action equation (2.5) is unchanged, but now $\boldsymbol {U}=\boldsymbol {u}_\psi$. That is, the previously arbitrary mean flow is now specifically identified with the velocity field $(-\psi _y,\psi _x)$ induced by the point vortices and wave packets.

The most interesting effect of the coupling and summation of Lagrangians is the generalization of (2.14) to (2.23)–(2.24). By these equations, the quantity $H_0\nabla ^2 \psi -\boldsymbol {\nabla } \times \boldsymbol {p}$ is conserved following the mean motion of fluid particles. Consider waves propagating into a region of fluid that is initially at rest. Before the arrival of the waves, $\nabla ^2 \psi =\boldsymbol {p}=0$, and hence

(2.26)\begin{equation} H_0\nabla^2 \psi = \boldsymbol{\nabla} \times \boldsymbol{p}. \end{equation}

By (2.23), (2.26) applies at all times, including when waves are present. Equation (2.26) is a concise definition of Bretherton flow, the flow generated by a wave packet in a formerly quiescent fluid. If wave breaking destroys the pseudomomentum $\boldsymbol {p}$ before the broad mean flow represented by $\psi$ has time to react, then real, actual, vorticity is created and remains behind after the remaining wave energy propagates away.

Taking the gradient of the dispersion relation (2.4) and using $\boldsymbol {\nabla } \omega = - \boldsymbol {\nabla } \theta _t= -\boldsymbol {k}_t$, we obtain the refraction equation

(2.27)\begin{equation} \frac{\partial \boldsymbol{k}}{\partial t} +( (\boldsymbol{c}_g+\boldsymbol{U})\boldsymbol{\cdot} \boldsymbol{\nabla} )\boldsymbol{k} ={-}k \boldsymbol{\nabla} U-l\boldsymbol{\nabla} V. \end{equation}

The refractive change in $\boldsymbol {k}$ predicted by (2.27) causes a change in $\omega _r(\boldsymbol {k})$ that can be determined from (2.3). If the waves do not break, then the action $\mathcal {A}=E/\omega _r(\boldsymbol {k})$ is conserved. If $\omega _r(\boldsymbol {k})$ increases, then the wave energy $E$ must also increase to keep their ratio constant. We anticipate that wave vector stretching, which increases $\vert \boldsymbol {k} \vert$ and hence $\omega _r$, is typical for the same reason that fluid particles typically move apart, and hence wave-mean interactions typically transfer energy from surface currents to waves. On the other hand, wave breaking always transfers energy from waves to currents.

Now we specialize the dynamics (2.15)–(2.16) to the case in which the vorticity and wave action are concentrated at widely separated points. This specialization is motivated by a desire to produce equations amenable to analytical and numerical solution. We assume that the mean flow vorticity consists solely of point vortices. Then

(2.28)\begin{equation} [\alpha,\beta]=\sum_i \varGamma_i \delta(\boldsymbol{x}-\boldsymbol{x}_i(t)), \end{equation}

where $\boldsymbol {x}_i(t)$ is the location at time $t$ of a point vortex with strength $\varGamma _i$. The subscript $i$ replaces the continuous vorticity labels $\alpha$ and $\beta$. The Hamiltonian (2.16) becomes

(2.29)\begin{align} H[\theta,\mathcal{A},\boldsymbol{x}_i,\psi] &={-}\sum_i H_0\varGamma_i \psi(\boldsymbol{x}_i(t))\nonumber\\ &\quad +\iint \textrm{d}\kern0.07em \boldsymbol{x} \left( \omega_r(\theta_x,\theta_y) \mathcal{A} -\frac{H_0}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \psi +[\psi,\theta] \mathcal{A} \right). \end{align}

To fully convert from $(\alpha ,\beta )$ to $\boldsymbol {x}_i$ we must transform the term

(2.30)\begin{equation} \iiint \textrm{d}\kern0.07em \boldsymbol{x} \,\textrm{d}t \,\alpha \beta_t \end{equation}

in (2.15). It becomes

(2.31)\begin{align} &\iiint {\textrm{d}\kern0.07em x} \,\textrm{d} y\, \textrm{d}t\, \alpha \frac{\partial (x,y,\beta)}{\partial(x,y,t)} =\iiint \textrm{d}\alpha \,\textrm{d}\beta \,\textrm{d}\tau \,\alpha \frac{\partial (x,y,\beta)}{\partial(\alpha,\beta,\tau)} \nonumber\\ &\quad =\iiint \textrm{d}\alpha \,\textrm{d}\beta \,\textrm{d}\tau -x \frac{\partial (\alpha,y,\beta)}{\partial(\alpha,\beta,\tau)} =\iiint \textrm{d}\alpha \,\textrm{d}\beta \,\textrm{d}\tau x\frac{\partial y}{\partial \tau} \nonumber\\ &\quad =\int \textrm{d}t \sum_i \varGamma_i x_i \frac{\textrm{d}y_i}{\textrm{d} t}. \end{align}

Thus, when point vortices replace continuous vorticity, the Lagrangian (2.15) becomes

(2.32)\begin{align} L[\theta,\mathcal{A},\boldsymbol{x}_i,\psi]&= \iiint \textrm{d}\kern0.07em \boldsymbol{x} \,\textrm{d}t \left( -\theta_t\mathcal{A} -\omega_r(\theta_x,\theta_y) \mathcal{A} +\frac{H_0}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} \psi -[\psi,\theta] \mathcal{A} \right) \nonumber\\ &\quad +\int \textrm{d}t \left( -\sum_i H_0\varGamma_i x_i \frac{\textrm{d}y_i}{\textrm{d}t} + \sum_i H_0\varGamma_i \psi(\boldsymbol{x}_i(t)) \right). \end{align}

Instead of (2.20) we now have

(2.33)\begin{equation} \delta \psi: \ H_0\nabla^2 \psi = \sum_i H_0\varGamma_i \delta(\boldsymbol{x}-\boldsymbol{x}_i(t))+[\mathcal{A},\theta], \end{equation}

with solution

(2.34)\begin{equation} \psi(\boldsymbol{x},t)=\frac{1}{2{\rm \pi}}\sum_i \varGamma_i \ln \vert\boldsymbol{x}-\boldsymbol{x}_i(t)\vert +\psi_w(\boldsymbol{x},t), \end{equation}

where (suppressing the time dependence)

(2.35)\begin{equation} \psi_w(\boldsymbol{x})=\frac{1}{H_0}\iint \textrm{d}\kern0.07em \boldsymbol{x}' \rho(\boldsymbol{x}') \frac{1}{2{\rm \pi}} \ln \vert \boldsymbol{x} - \boldsymbol{x}' \vert \end{equation}

and

(2.36)\begin{equation} \rho =[\mathcal{A},\theta] = \boldsymbol{\nabla} \mathcal{A} \times \boldsymbol{k}. \end{equation}

We now assume that the wave field consists solely of isolated wave packets. The stream function field generated by a single wave packet at $\boldsymbol {x}_p$ is given by (2.35) and (2.36) with $\boldsymbol {k}=\boldsymbol {k}_p$, where $\boldsymbol {k}_p(t)$ is the wave vector associated with the wave packet. We assume that $\boldsymbol {k}_p$ depends only on time; its variation within the wave packet is assumed negligible. The integration in (2.35) is over the area of the wave packet, the region of the flow in which $\mathcal {A}\neq 0$. In Appendix B we show that, far from $\boldsymbol {x}_p$, the stream function generated by a wave packet at $\boldsymbol {x}_p$ takes the form of a dipole:

(2.37)\begin{equation} \psi(\boldsymbol{x})= \frac{1}{2 {\rm \pi}H_0} \frac {(\boldsymbol{x}-\boldsymbol{x}_p) \times \boldsymbol{k}_p}{\vert \boldsymbol{x}-\boldsymbol{x}_p \vert^2} \mathcal{A}_p, \end{equation}

where $\mathcal {A}_p = \iint \textrm{d}\kern1.5pt\boldsymbol{x} \mathcal {A}$ is the total action of the wave packet.

The stream function response to many point vortices and many wave packets is clearly

(2.38)\begin{equation} \psi(\boldsymbol{x},t)=\sum_i \varGamma_i \psi_m(\boldsymbol{x},\boldsymbol{x}_i) + \sum_p \mathcal{A}_p \psi_d(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p), \end{equation}

where

(2.39)\begin{equation} \psi_m(\boldsymbol{x},\boldsymbol{x}_i)\equiv \frac{1}{2{\rm \pi}}\ln \vert \boldsymbol{x}-\boldsymbol{x}_i(t)\vert \end{equation}

is the response to a monopole at $\boldsymbol {x}_i$ and

(2.40)\begin{equation} \psi_d(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p)\equiv \frac{1}{2{\rm \pi} H_0}\frac {(\boldsymbol{x}-\boldsymbol{x}_p) \times \boldsymbol{k}_p}{\vert \boldsymbol{x}-\boldsymbol{x}_p \vert^2} \end{equation}

is the response to a dipole with wave vector $\boldsymbol {k}_p$ at $\boldsymbol {x}_p$. The constants $\varGamma _i$ and $\mathcal {A}_p$ measure the strength of the monopole and the dipole, respectively. Constant $\mathcal {A}_p$ is always positive but $\varGamma _i$ can have either sign. Until dissipation occurs $\mathcal {A}_p$ and $\varGamma _i$ remain constant.

Since our aim is to produce a Lagrangian that depends only on the point vortex locations $\boldsymbol {x}_i(t)$, the wave packet locations $\boldsymbol {x}_p(t)$ and their wave vectors $\boldsymbol {k}_p(t)$, we must transform all of the terms in (2.32). If we integrate the first term in (2.32) over the $p$th wave packet, we obtain

(2.41)\begin{align} -\iiint \textrm{d}\kern0.07em \boldsymbol{x} \,\textrm{d}t \theta_t\mathcal{A} &=\iiint \textrm{d}\kern0.07em \boldsymbol{x}\,\textrm{d}t \theta \mathcal{A}_t ={-}\iiint \textrm{d}\kern0.07em \boldsymbol{x} \,\textrm{d}t \theta \frac{\textrm{d}\kern0.07em \boldsymbol{x}_p}{\textrm{d}t}\boldsymbol{\cdot} \boldsymbol{\nabla} \mathcal{A}\nonumber\\ &={-}\int \textrm{d}t \frac{\textrm{d}\kern0.07em \boldsymbol{x}_p}{\textrm{d}t} \boldsymbol{\cdot} \iint \textrm{d}\kern0.07em \boldsymbol{x} \theta \boldsymbol{\nabla} \mathcal{A} =\int \textrm{d}t \frac{\textrm{d}\kern0.07em \boldsymbol{x}_p}{\textrm{d}t} \boldsymbol{\cdot} \iint \textrm{d}\kern0.07em \boldsymbol{x} \mathcal{A} \boldsymbol{\nabla} \theta \nonumber\\ &= \int \textrm{d}t \frac{\textrm{d}\kern0.07em \boldsymbol{x}_p}{\textrm{d}t} \boldsymbol{\cdot} \boldsymbol{k}_p \mathcal{A}_p, \end{align}

where we have used integrations by parts and the relation

(2.42)\begin{equation} \left( \frac{\partial}{\partial t} + \frac{\textrm{d}\kern0.07em \boldsymbol{x}_p}{\textrm{d}t} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \mathcal{A}(\boldsymbol{x},t)=0, \end{equation}

which follows from the definition of $\boldsymbol {x}_p(t)$: $\textrm {d}\kern1.5pt\boldsymbol {x}_p(t)/\textrm {d}t$ is the velocity of the wave envelope. The second term in (2.32) becomes

(2.43)\begin{equation} -\int \textrm{d}t \omega_r(\boldsymbol{k}_p) \mathcal{A}_p. \end{equation}

The three terms in (2.32) containing $\psi$ combine as

(2.44)\begin{align} & \iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x} \left( \frac{H_0}{2} \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla}\psi -[\psi,\theta] \mathcal{A} \right) + \int \textrm{d}t \sum_i H_0\varGamma_i \psi(\boldsymbol{x}_i(t),t) \nonumber\\ &\quad = \iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x} \left( -\frac{H_0}{2} \psi \nabla^2 \psi -[\psi,\theta] \mathcal{A} + \sum_i H_0\varGamma_i \psi \delta(\boldsymbol{x}-\boldsymbol{x}_i) \right) \nonumber\\ &\quad =\iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x} \left( -\frac{H_0}{2} \psi \nabla^2 \psi +[\mathcal{A},\theta] \psi + \psi ( H_0\nabla^2 \psi - [\mathcal{A},\theta] ) \right) \nonumber\\ &\quad = \frac{H_0}{2} \iiint \textrm{d}t\,\textrm{d}\kern0.07em \boldsymbol{x} \psi \nabla^2 \psi, \end{align}

where we have used (2.33).

Our final step is to substitute (2.38) back into the Lagrangian, removing its dependence on $\psi$. The last integral in (2.44) becomes

(2.45)\begin{align} \iint \textrm{d}\kern0.07em \boldsymbol{x} \psi \nabla^2 \psi &= \iint \textrm{d}\kern0.07em \boldsymbol{x} \left[ \left( \sum_i \varGamma_i \psi_m(\boldsymbol{x},\boldsymbol{x}_i) + \sum_p \mathcal{A}_p \psi_d(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p) \right) \right. \nonumber\\ &\quad \times \left.\nabla^2 \left( \sum_j \varGamma_j \psi_m(\boldsymbol{x},\boldsymbol{x}_j) + \sum_q \mathcal{A}_q \psi_d(\boldsymbol{x},\boldsymbol{x}_q,\boldsymbol{k}_q) \right) \right]. \end{align}

We simplify (2.45) by neglecting the dipole–dipole interactions, which are expected to be weak: the velocity field associated with the monopoles falls off like $1/r$, whereas the velocity field associated with Bretherton dipoles falls off like $1/r^2$. Dropping these terms from (2.45) gives us

(2.46)\begin{align} \iint \textrm{d}\kern0.07em \boldsymbol{x} \psi \nabla^2 \psi\, &\approx \iint \textrm{d}\kern0.07em \boldsymbol{x} \sum_i \varGamma_i \nabla^2 \psi_m(\boldsymbol{x},\boldsymbol{x}_i) \left( \sum_j \varGamma_j \psi_m(\boldsymbol{x},\boldsymbol{x}_j) +2 \sum_p \mathcal{A}_p \psi_d(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p) \right) \nonumber\\ &=\iint \textrm{d}\kern0.07em \boldsymbol{x} \sum_i \varGamma_i \delta (\boldsymbol{x}-\boldsymbol{x}_i) \left( \sum_j \varGamma_j \psi_m(\boldsymbol{x},\boldsymbol{x}_j) +2 \sum_p \mathcal{A}_p\psi_d(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p) \right) \nonumber\\ &= \sum_i \varGamma_i \left( \sum_j \varGamma_j \psi_m(\boldsymbol{x_i},\boldsymbol{x}_j) +2 \sum_p \mathcal{A}_p \psi_d(\boldsymbol{x_i},\boldsymbol{x}_p,\boldsymbol{k}_p) \right). \end{align}

Putting all this together, we obtain the Lagrangian

(2.47)\begin{equation} L[ \boldsymbol{x}_i,\boldsymbol{x}_p,\boldsymbol{k}_p]= \int \textrm{d}t \left( \sum_p \mathcal{A}_p \boldsymbol{k}_p \boldsymbol{\cdot} \dot{\boldsymbol{x}}_p - \sum_i H_0\varGamma_i x_i \dot{y}_i -H[\boldsymbol{x}_i,\boldsymbol{x}_p,\boldsymbol{k}_p] \right), \end{equation}

where

(2.48)\begin{align} H[\boldsymbol{x}_i,\boldsymbol{x}_p,\boldsymbol{k}_p] &= \sum_p \mathcal{A}_p \omega_r(\boldsymbol{k}_p) - \frac{H_0}{2{\rm \pi}} \sum_i \sum_{j>i} \varGamma_i \varGamma_j \ln \vert \boldsymbol{x}_i -\boldsymbol{x}_j \vert \nonumber\\ &\quad -\frac{1}{2{\rm \pi}} \sum_i \sum_p \varGamma_i \mathcal{A}_p \frac{(\boldsymbol{x}_i-\boldsymbol{x}_p) \times \boldsymbol{k}_p} {\vert \boldsymbol{x}_i-\boldsymbol{x}_p \vert ^2} \end{align}

is the Hamiltonian. For every wave packet there are two canonical pairs, $(x_p,k_p)$ and $(l_p,y_p)$, and for every point vortex there is one canonical pair, $(x_i,y_i)$. Again, $\varGamma _i$ and $\mathcal {A}_p$ are constants. The Hamiltonian (2.48) contains $\varGamma \varGamma$ terms and $\varGamma \mathcal {A}$ terms. If we had not dropped the dipole–dipole interactions, it would also contain $\mathcal {A} \mathcal {A}$ terms.

We remark that it is generally quite wrong to substitute an equation resulting from the variational principle back into the Lagrangian. If, for example, we substitute the dispersion relation back into (2.1), the Lagrangian vanishes. However, it is legitimate to use the equation obtained by varying a particular field to eliminate that same field from the Lagrangian; see Appendix C. Thus it is acceptable to use (2.33) to eliminate $\psi$ from (2.32).

The equations corresponding to (2.47) and (2.48) are

(2.49)\begin{gather} \delta \boldsymbol{k}_p: \ \dot{\boldsymbol{x}}_p=\frac{1}{\mathcal{A}_p} \frac{\partial H}{\partial \boldsymbol{k}_p} =\boldsymbol{c}_g(\boldsymbol{k}_p) +\boldsymbol{U}_m(\boldsymbol{x}_p), \end{gather}

(2.50)\begin{gather}\delta \boldsymbol{x}_p: \ \dot{\boldsymbol{k}}_p={-}\frac{1}{\mathcal{A}_p} \frac{\partial H}{\partial \boldsymbol{x}_p} ={-}k_p \boldsymbol{\nabla} U_m(\boldsymbol{x}_p)-l_p\boldsymbol{\nabla} V_m(\boldsymbol{x}_p), \end{gather}

(2.51)\begin{gather}\delta \boldsymbol{x}_i: \ \dot{\boldsymbol{x}}_i=\frac{1}{\varGamma_i} \left( \frac{\partial H}{\partial y_i}, -\frac{\partial H}{\partial x_i} \right) =\boldsymbol{U}_m(\boldsymbol{x}_i)+\boldsymbol{U}_d(\boldsymbol{x}_i), \end{gather}

where

(2.52)\begin{align} \boldsymbol{U}_m(\boldsymbol{x})=(U_m(\boldsymbol{x}),V_m(\boldsymbol{x})) &=\sum_i \varGamma_i \left( -\frac{\partial \psi_m}{\partial y}(\boldsymbol{x},\boldsymbol{x}_i), \frac{\partial \psi_m}{\partial x}(\boldsymbol{x},\boldsymbol{x}_i) \right) \nonumber\\ &=\frac{1}{2{\rm \pi}} \sum_i \varGamma_i \frac{(y_i-y,x-x_i)}{\vert \boldsymbol{x}_i-\boldsymbol{x} \vert ^2} \end{align}

is the velocity field induced by the point vortices and

(2.53)\begin{equation} \boldsymbol{U}_d(\boldsymbol{x})=\sum_p \mathcal{A}_p \left( -\frac{\partial \psi_d}{\partial y}(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p), \frac{\partial \psi_d}{\partial x}(\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p) \right) \end{equation}

is the velocity field induced by the wave packets. The total velocity is $\boldsymbol {U}(\boldsymbol {x})=\boldsymbol {U}_m(\boldsymbol {x})+\boldsymbol {U}_d(\boldsymbol {x})$. In our approximation, the wave packets talk to point vortices but not to one another, while the point vortices talk to both point vortices and wave packets. We can add the missing physics if necessary; it would, for example, add the term $\boldsymbol {U}_d(\boldsymbol {x}_p)$ to (2.49).

Equations (2.49)–(2.51) are the fundamental equations of our model. If we were to regard $\boldsymbol {U}_m$ as a prescribed mean flow, then (2.49) and (2.50) would be the standard equations of ray theory (e.g. Bühler Reference Bühler2014). Similarly, if we omit $\boldsymbol {U}_d$, then (2.51) is the standard equation of point vortex dynamics (Kirchhoff Reference Kirchhoff1883). The new feature of our derivation is that $\boldsymbol {U}_m$ is not prescribed, but rather is determined by the locations of the point vortices. Similarly, the dipolar velocity field of the wave packets is not dropped, but rather contributes to the advection of the point vortices. Again, if the relatively weak interactions between the wave packets had not been dropped, then $\boldsymbol {U}_d$ would also appear in (2.49) and (2.50). Tchieu, Kanso & Newton (Reference Tchieu, Kanso and Newton2012) consider a system consisting solely of interacting point dipoles. In our context, their system corresponds to adding dipole–dipole interactions but completely omitting the point vortices.

The derivation of (2.49)–(2.51) from a Lagrangian guarantees that our dynamics maintains important conservation laws. The conservation of energy (2.48) corresponds to the time-translation symmetry of (2.47)–(2.48). The conservation of momentum

(2.54)\begin{equation} \mathcal{M}=\sum_p \mathcal{A}_p \boldsymbol{k}_p + H_0\sum_i\varGamma_i(y_i,-x_i) \end{equation}

corresponds to space-translation symmetry and is proved by considering variations of the form

(2.55)\begin{equation} \delta \boldsymbol{x}_i =\delta \boldsymbol{x}_p= \boldsymbol{\epsilon}(t), \end{equation}

where $\boldsymbol {\epsilon }(t)$ is an arbitrary infinitesimal vector. If we think of the interactions between the dipoles and point vortices as the sum of pair interactions between each dipole–vortex pair, then pairwise conservation of (2.54) shows that the refraction of wave packet $p$ (i.e. the change in $\boldsymbol {k}_p$) caused by vortex $i$ is accompanied by a change in the position of vortex $i$. Bühler & McIntyre (Reference Bühler and McIntyre2003) refer to this as ‘remote recoil’. Conservation of (2.54) also governs wave breaking in the following sense. If the $p$th wave packet is completely destroyed by wave breaking, then $\mathcal {A}_p$ is suddenly replaced by two counter-rotating vortices with a dipole moment equal to $\varGamma D$, where $D$ is the separation between counter-rotating vortices of strength $\pm \varGamma$. See also Bühler & McIntyre (Reference Bühler and McIntyre2005), Bühler & Jacobson (Reference Bühler and Jacobson2001) and Bühler (Reference Bühler2014). Our dynamics also conserves the angular momentum:

(2.56)\begin{equation} \mathcal{L}=\sum_p \mathcal{A}_p (\boldsymbol{k}_p \times \boldsymbol{x}_p) + \frac{H_0}{2} \sum_i \varGamma_i (x_i^2 + y_i^2), \end{equation}

which can be proved by considering variations of the form $(x_i+\textrm {i}y_i) \rightarrow (x_i+\textrm {i}y_i)\, \textrm {e}^{\textrm {i}\delta \theta (t)}$, $(x_p+\textrm {i}y_p) \rightarrow (x_p+\textrm {i}y_p)\, \textrm {e}^{\textrm {i}\delta \theta (t)}$ and $(k_p+\textrm {i}l_p) \rightarrow (k_p+\textrm {i}l_p)\,\textrm {e}^{\textrm {i}\delta \theta (t)}$, where $\delta \theta (t)$ is an infinitesimal angle.

The velocity field induced by the point vortices falls off like $r^{-1}$, where $r$ is the distance from the vortex, whereas the velocity field induced by the wave packets falls off at the faster rate $r^{-2}$. In reality, the vorticity associated with a surface wave packet resides in a horseshoe-shaped vortex tube whose surface manifestation is the vortex pair represented by our dipole. In semi-infinite three-dimensional geometry, the velocity induced by the horseshoe-shaped vortex falls off like $r^{-3}$, where $r$ is the distance from the wave packet. This follows from the Biot–Savart law relating vorticity to velocity. However, in our model, the assumption of a wave layer of finite thickness, with a rigid lower boundary at depth $H_0$, converts the $r^{-3}$ fall-off to $r^{-2}$, because the rigid mixed-layer bottom confines the Biot–Savart response to the layer. The need to invoke a rigid lower mixed-layer boundary demotes our model to something of a toy model. However, some such strong assumption is unavoidable if one intends to model the wave layer as a two-dimensional system. Despite the severity of this assumption, we believe that the two-dimensional model captures enough of the physics to be a useful thinking tool and that some of our conclusions will survive generalization to a more inclusive three-dimensional model.

Onsager (Reference Onsager1949) considered the equilibrium statistical mechanics of a system of point vortices. Our system reduces to Onsager's system when no waves are present ($\mathcal {A}_p\equiv 0$). Our phase space is larger than that considered by Onsager because it contains dimensions corresponding to the wave packet locations $\boldsymbol {x}_p$ and their wave vectors $\boldsymbol {k}_p$. However, the difference is not merely a matter of extra dimensions. In Onsager's problem the volume of the phase space is finite, because the point vortices are confined to a box. In our problem the phase space has infinite volume because $-\infty <\boldsymbol {k}_p <\infty$. We therefore expect an ultraviolet catastrophe in which energy spreads to ever larger $\vert \boldsymbol {k}_p \vert$ by the process of wave vector stretching. If wave vector stretching increases the first term in (2.48), as would be the case for surface waves, this increase must be compensated by a decrease in the other two terms.

Our method is easily adapted to other types of waves and mean flows. For example, to investigate internal waves interacting with a quasigeostrophic mean flow, we need only replace (2.4) with the dispersion relation for internal waves, and (2.7) with the Lagrangian for quasigeostrophic flow. This approach offers advantages of simplicity and transparency over the more formal approaches followed by Bühler & McIntyre (Reference Bühler and McIntyre2005), Wagner & Young (Reference Wagner and Young2015) and Salmon (Reference Salmon2016). For many further details, see Salmon (Reference Salmon2020). In the remainder of this paper we investigate the dynamics (2.49)–(2.51).

3. One vortex, one wave packet

We begin by considering the system consisting of a single vortex of strength $\varGamma$ located at $\boldsymbol {x}(t)$, and a single wave packet with action $\mathcal {A}_p$ and wave vector $\boldsymbol {k}(t)$ located at $\boldsymbol {x}(t)+\boldsymbol {\xi }(t)$. This system exhibits a much more complicated range of behaviour than the more familiar system consisting of two point vortices. The Langrangian (2.47) takes the form

(3.1)\begin{equation} L[ \boldsymbol{x},\boldsymbol{\xi},\boldsymbol{k}]= \int \textrm{d}t ( \mathcal{A}_p \boldsymbol{k} \boldsymbol{\cdot} ( \dot{\boldsymbol{\xi}} +\dot{ \boldsymbol{x}}) - H_0 \varGamma x \dot{y} -H(\boldsymbol{\xi},\boldsymbol{k}) ), \end{equation}

with Hamiltonian

(3.2)\begin{equation} H(\boldsymbol{\xi},\boldsymbol{k})= \mathcal{A}_p \sqrt{g \vert \boldsymbol{k} \vert}- \frac{1}{2{\rm \pi}} \varGamma \mathcal{A}_p \frac{\boldsymbol{k} \times \boldsymbol{\xi}} {\vert \boldsymbol{\xi} \vert^2}. \end{equation}

The equations of motion become

(3.3)\begin{gather} \delta x : \ \frac{\textrm{d}}{\textrm{d}t}(\mathcal{A}_pk+H_0 \varGamma y)=0, \end{gather}

(3.4)\begin{gather}\delta y : \ \frac{\textrm{d}}{\textrm{d}t}(\mathcal{A}_pl-H_0 \varGamma x)=0, \end{gather}

(3.5)\begin{gather}\delta \boldsymbol{k} : \ \frac{\textrm{d}}{\textrm{d}t} (\boldsymbol{\xi}+\boldsymbol{x}) =\boldsymbol{c}_g(\boldsymbol{k})+ \frac{\varGamma}{2{\rm \pi} \vert \boldsymbol{\xi} \vert^2}(-\eta,\xi), \end{gather}

(3.6)\begin{gather}\delta \xi : \ \frac{\textrm{d}k}{\textrm{d}t}=\frac{\varGamma}{2{\rm \pi} \vert \boldsymbol{\xi} \vert^4} (l(\xi^2-\eta^2)-2k \xi \eta), \end{gather}

(3.7)\begin{gather}\delta \eta : \ \frac{\textrm{d}l}{\textrm{d}t}=\frac{\varGamma}{2{\rm \pi} \vert \boldsymbol{\xi} \vert^4} ( k(\xi^2-\eta^2)+2l \xi \eta ), \end{gather}

where $\boldsymbol {\xi }=(\xi ,\eta )$. We simplify notation by taking $g=1$ and choosing a characteristic wavenumber $k_0=1$ so that $\omega _0^2=gk_0=1$. We also assume $H_0=k_0^{-1}=1$, while we take $\varGamma = 2{\rm \pi} \omega _0 k_0^{-2}=2{\rm \pi}$ and $\mathcal {A}_p=2{\rm \pi} \omega _0k_0^{-5} = 2{\rm \pi}$. Then the Hamiltonian (3.2) becomes

(3.8)\begin{equation} H = 2{\rm \pi}\left((k^2+l^2)^{1/4}- \frac{k\eta - l \xi}{\xi^2+\eta^2}\right). \end{equation}

The system (3.3)–(3.7) conserves the energy (3.8), the angular momentum

(3.9)\begin{equation} \mathcal{L} = 2{\rm \pi} (\boldsymbol{k}\times (\boldsymbol{x}+\boldsymbol{\xi})) + \frac{2{\rm \pi}}{2}(x^2+y^2) \end{equation}

(cf. (2.56)) and the momentum $\mathcal {M}\equiv (\mathcal {M}_x,\mathcal {M}_y)$ (cf. (2.54)), where

(3.10a,b)\begin{equation} \mathcal{M}_x = 2{\rm \pi} (k+ y); \quad \mathcal{M}_y=2{\rm \pi}(l- x). \end{equation}

We use the conserved momenta (3.10a,b) to eliminate the variables $(x,y)$ in favour of $(k,l,\xi ,\eta )$. The resulting system conserves the energy (3.8) and the quantity

(3.11)\begin{equation} \mathcal{R}_0 \equiv \frac{1}{2}\left(\left(\frac{\mathcal{M}_x}{2{\rm \pi}}\right)^2+\left(\frac{\mathcal{M}_y}{2{\rm \pi}}\right)^2\right) -\frac{\mathcal{L}}{2{\rm \pi}} = \frac{1}{2}(k^2+l^2)- (k\eta - l\xi) \end{equation}

obtained by eliminating $(x,y)$ between (3.9) and (3.10a,b). We also define

(3.12)\begin{equation} \mathcal{H}_0 = \frac{H}{2{\rm \pi}}. \end{equation}

The reduced dynamics takes the form of four coupled ordinary differential equations:

(3.13)\begin{gather} \dot{\xi}={-}\frac{ k(\xi^2-\eta^2)+2l\xi \eta}{(\xi^2+\eta^2)^2}+\frac{k}{2(k^2+l^2)^{3/4}}- \frac{\eta }{\xi^2 + \eta^2}, \end{gather}

(3.14)\begin{gather}\dot{\eta} = \frac{ l(\xi^2-\eta^2)-2k\xi \eta}{(\xi^2+\eta^2)^2}+\frac{l}{2(k^2+l^2)^{3/4}}+ \frac{\xi }{\xi^2 + \eta^2}, \end{gather}

(3.15)\begin{gather}\dot{k} = \frac{ l(\xi^2-\eta^2)-2k\xi \eta}{(\xi^2+\eta^2)^2}, \end{gather}

(3.16)\begin{gather}\dot{l} = \frac{ k(\xi^2-\eta^2)+2l\xi \eta}{(\xi^2+\eta^2)^2}, \end{gather}

with the two conserved quantities (3.11) and (3.12).

Define

(3.17a,b)\begin{equation} k+\mathrm{i} l = \kappa\, \textrm{e}^{\mathrm{i}\phi}, \quad \xi+ \mathrm{i} \eta = a \,\textrm{e}^{\mathrm{i}\theta}. \end{equation}

We shall obtain a single, closed equation for the wavenumber magnitude $\kappa (t)$. First, using (3.8), (3.11) and (3.12), we obtain an expression for $a^2$ in terms of $\kappa$:

(3.18)\begin{equation} a^2= \frac{\mathcal{R}_0-\frac{1}{2}\kappa ^2}{\mathcal{H}_0-\sqrt{\kappa}}. \end{equation}

Then using (3.8) we obtain the constraint

(3.19)\begin{equation} \sin (\phi-\theta) = \frac{a}{\kappa} (\mathcal{H}_0-\sqrt{\kappa}) \end{equation}

on the phases. From (3.15) and (3.16), we find an evolution equation for $\phi$, namely

(3.20)\begin{equation} \dot{\phi}=\frac{-l\dot{k}+k\dot{l}}{k^2+l^2}=\frac{1}{a^2} \cos 2(\phi-\theta). \end{equation}

Equations (3.15) and (3.16) also imply

(3.21)\begin{equation} \dot{\kappa} = \frac{\dot{k}k+\dot{l}l}{\sqrt{k^2+l^2}} = \frac{\kappa}{a^2} \sin 2(\phi-\theta). \end{equation}

Combining (3.20) and (3.21), we obtain

(3.22)\begin{equation} \dot{\kappa}^2+\kappa^2 \dot{\phi}^2 = \frac{\kappa^2}{a^4} = \kappa^2 \left(\frac{\mathcal{H}_0-\sqrt{\kappa}}{\mathcal{R}_0-\frac{1}{2}\kappa ^2}\right)^2. \end{equation}

Our final step is to eliminate $\dot {\phi }$ to arrive at an equation involving only $\dot {\kappa }$ and $\kappa$. We use the identity

(3.23)\begin{equation} \cos 2(\phi - \theta) =1- 2\sin^2 (\phi-\theta) = 1-\frac{2 a^2}{\kappa^2} (\mathcal{H}_0-\sqrt{\kappa})^2, \end{equation}

where the last substitution is via (3.18) and (3.19). Then by (3.18) and (3.23) we have

(3.24)\begin{equation} \dot{\phi} = \frac{\mathcal{H}_0-\sqrt{\kappa}}{\mathcal{R}_0-\frac{1}{2}\kappa ^2}-\frac{2}{\kappa^2} (\mathcal{H}_0-\sqrt{\kappa})^2. \end{equation}

Substituting (3.24) back into (3.22) we obtain the closed evolution equation

(3.25)\begin{equation} \tfrac{1}{2}\dot{\kappa}^2+\varPi(\kappa)=0 \end{equation}

for $\kappa (t)$, where

(3.26)\begin{equation} \varPi(\kappa)=\frac{1}{2}\left(\frac{\kappa (\mathcal{H}_0-\sqrt{\kappa})}{(\mathcal{R}_0-\frac{1}{2}\kappa^2)}\right)^2\left({-}1+ \left(1-\frac{2}{ \kappa^2}\left(\mathcal{R}_0-\frac{1}{2}\kappa^2\right)(\mathcal{H}_0-\sqrt{\kappa})\right)^2 \right). \end{equation}

Equation (3.25) takes the form of a particle moving in a potential $\varPi (\kappa )$. This permits a qualitative analysis of system behaviour based on the form of (3.26). Solutions may be written out in implicit form, as in Tur & Yanovsky (Reference Tur and Yanovsky2017), but a qualitative analysis offers better physical insight.

3.1. Circular motion

We begin by seeking solutions that exhibit simple harmonic motion. Thus we take $\dot {\kappa }=0$ and look for the $\kappa _i$ that satisfy

(3.27)\begin{equation} \varPi(\kappa_i) = 0. \end{equation}

Let $\kappa _i\equiv 1$. This implies a simple relation between $\mathcal {H}_0$ and $\mathcal {R}_0$. Its solutions are $\mathcal {H}_0=1$ and $\mathcal {R}_0$ a free parameter; or

(3.28)\begin{equation} \mathcal{H}_0=\frac{1+2\mathcal{R}_0}{-1+2\mathcal{R}_0}. \end{equation}

If we take $\kappa _i= 1$ to be a critical point, then $\varPi '(1)=0$. It may be shown that when $\mathcal {H}_0=1$, $\varPi ''(1)=0$. Therefore, in order to exhibit unstable and stable solutions, we consider the set of solutions described by (3.28). This leads to the two possibilities

(3.29)\begin{equation} \mathcal{R}^{{\pm}}_0 \equiv \tfrac{1}{2}({-}3\pm 2\sqrt{2}). \end{equation}

The solutions take the form

(3.30a–d)\begin{equation} k=\cos \phi , \quad l=\sin \phi, \quad \xi =a \cos \theta , \quad \eta= a \sin \theta , \end{equation}

where $a$ is given by (3.18) and $\phi$ and $\theta$ are found from (3.19) and (3.20).

An example of this behaviour is shown in figure 1, where the two potentials and corresponding solutions are shown. In one of these solutions the wave packet orbit lies inside the orbit of the vortex. In the other solution, the opposite occurs.

Figure 1. (a) The potentials $\varPi (\kappa )$ for the two cases ($\mathcal {R}^{\pm }_0$) in which a critical point is present at $\kappa =1$. The inset enlarges the potentials near $\kappa =1$. In one case (dashed) the critical point coincides with a maximum of $\varPi (\kappa )$ (implying unstable motion) while in the other case (solid) the critical point is a minimum. (b,c) Numerical solutions showing the locations of the wave packet (black line) and vortex (red line) in the two cases. The initial wave vector is indicated by the black arrow.

3.2. Stability of orbits

The $\varPi (\kappa )$ graphed in figure 1(a) suggest that the circular orbits shown there may not be locally stable (in a spectral sense) to perturbations. Therefore we study solutions in the neighbourhood of $\kappa _0 = 1$. We take $\kappa =1 +\epsilon \kappa _1$ and expand $\varPi (\kappa )$ about the critical point $\kappa =1$, to find

(3.31)\begin{equation} \varPi(\kappa) \approx \varPi(1)+\epsilon \kappa_1 \varPi'(1) + \epsilon^2 \kappa_1^2 \varPi''(1)+\cdots. \end{equation}

By construction, $\varPi (1) = \varPi '(1) = 0$. Taking $\kappa _1 = \kappa ^0_1\, \textrm {e}^{\lambda t}$, we find that the spectral stability of the system will be set by the sign of $\epsilon ^2 (\kappa _1^0)^2 \varPi ''(1)$. From figure 1 we see that $\mathcal {R}^+_0$ corresponds to stable orbits and $\mathcal {R}^-_0$ corresponds to unstable orbits. This is demonstrated in figure 2, which shows that the stable orbits are confined to the neighbourhood of their initial trajectories, whereas the unstable orbits deviate considerably.

Figure 2. (a) Orbit of wave packet and vortex with $\mathcal {R}_0 = \mathcal {R}_0^+$ as defined in (3.30a–d). For this case perturbations are stable, and the orbits remain close to the unperturbed solution. (b) Orbit of wave packet and vortex with $\mathcal {R}_0 = \mathcal {R}_0^-$. This orbit is unstable, and the solutions deviate considerably from circles.

3.3. Collapse

The above analysis addresses local spectral stability at a critical point. There are other solutions in which $a \to 0$ so that the wave packet and the vortex overlap. We call this phenomenon ‘collapse’. Collapsed solutions violate the assumption of our model that the wave packets and vortices remain far apart. Nonetheless, collapse is a real property of our equations that demands investigation. Vortex collapse for three point vortices has been extensively studied (see Aref (Reference Aref1983) and references therein). The case of overlapping vorticity and wave action has been analysed by McIntyre (Reference McIntyre2019).

The conditions for collapse are clear from (3.18). Collapse occurs at the time $t^*$ at which

(3.32)\begin{equation} \kappa^2 \to 2\mathcal{R}_0. \end{equation}

As an example we suppose that $\mathcal {H}_0=\mathcal {R}_0=0$. Then the system collapses as $\kappa \to 0$. Under these assumptions, the governing equation for $\kappa$ reduces to

(3.33)\begin{equation} \dot{\kappa}={\pm} 2\sqrt{\frac{2-\sqrt{\kappa}}{\sqrt{\kappa}}}. \end{equation}

Define $\vartheta$ by

(3.34)\begin{equation} \tan \vartheta ={\pm} \frac{\kappa^{1/4}}{\sqrt{2-\sqrt{\kappa}}}. \end{equation}

Then

(3.35a,b)\begin{equation} \cos \vartheta ={\pm} \frac{1}{\sqrt{2}}\sqrt{2-\sqrt{\kappa}}, \quad \sin \vartheta = \frac{1}{\sqrt{2}}\kappa^{1/4} \end{equation}

and, solving for $\kappa$, we obtain

(3.36)\begin{equation} \kappa = \tfrac{1}{2}(3-4\cos 2\vartheta +\cos 4\vartheta). \end{equation}

To find $t=t(\vartheta )$, we note that

(3.37)\begin{equation} \frac{\textrm{d}t}{\textrm{d}\vartheta}=\frac{\textrm{d}t}{\textrm{d}\kappa}\frac{\textrm{d}\kappa}{\textrm{d}\vartheta} ={\pm}\tan \vartheta(2\sin 2\vartheta -\sin 4\vartheta ). \end{equation}

This can be integrated, and we arrive at

(3.38)\begin{equation} t=t_0 \pm \tfrac{1}{4}(12\vartheta -8\sin 2\vartheta + \sin 4\vartheta ). \end{equation}

Collapse occurs when

(3.39)\begin{equation} \frac{\textrm{d} \kappa}{\textrm{d} \vartheta}=\frac{\textrm{d}t}{\textrm{d}\vartheta} = 0. \end{equation}

Thus the collapse corresponds to the formation of a cusp in $\kappa$.

Figure 3 confirms these results. Figure 3(b) shows the convergence of the wave packet and the vortex. Figure 3(c) compares the theoretical prediction of $\kappa (t)$ (where we have taken the negative branch of the solution corresponding to $\dot {\kappa }<0$) with the numerical result. The two curves are indistinguishable.

Figure 3. An example of collapse, in which the wave packet and the vortex converge, violating model assumptions. (a) The potential $\varPi (\kappa )$ with $\mathcal {R}_0=0$ and $\mathcal {H}_0=0$. Note the singularity at $\kappa =0$. (b) Converging particle paths. (c) The evolution of $\kappa (t)$, which vanishes in a cusp. The analytic solution is shown by the red dashed line and is indistinguishable from the numerical result.

In our model the wave action $\mathcal{A}_{p}$ is fixed. Therefore, the wave energy $\mathcal{A}_{p} \; \omega (\kappa )$ vanishes as $\kappa \to 0$ since $\omega (\kappa ) \propto \sqrt {\kappa }$ for surface gravity waves. The energy lost by the wave packet appears as an increase in the ‘interaction energy’ between the wave packet and the vortex – an increase in the last term in (3.2) – but, again, the whole theory breaks down when the two particles finally converge.

3.4. Blow-up

There are also solutions in which $\kappa$, the wavenumber modulus, grows without bound. We call these blow-up solutions. They correspond to wave packets that steepen. In reality, wave breaking limits the steepness of waves, and could be added to our model to extend its validity. For example, wave packets that exceed a prescribed steepness could be replaced by counter-rotating vortices with a dipole moment determined by momentum conservation (2.54) as in Bühler & Jacobson (Reference Bühler and Jacobson2001). In this paper we consider only ideal solutions, and we do not include wave breaking.

The particle in a potential well analogy implies that $\kappa (t)$ may grow without bound when $\varPi$ is a monotonically non-increasing function for large $\kappa$. As an example, we take $\mathcal {H}_0=-1$. Then for large $\kappa$:

(3.40)\begin{equation} \varPi \sim 2+O(1/\sqrt{\kappa}), \end{equation}

which implies that the blow-up solution takes the form

(3.41)\begin{equation} \kappa = \kappa_0 + \kappa_1 t. \end{equation}

We examine this numerically in figure 4, and find agreement with the theoretical prediction.

Figure 4. An example of a solution that ‘blows up’, meaning $\kappa \to \infty$. (a) The trajectories of the wave packet and vortex. (b) Plot of $\kappa$. The asymptotic form of the growth is predicted to go like $t$, which is shown by the dashed red line and is seen to agree well with the numerical integration.

4. Two vortices and one wave packet

We now examine the system comprising a single wave packet with action $\mathcal {A}_p$ and wave vector $(k_p, l_p)$ located at $(x_p, y_p)$; a point vortex of strength $-\varGamma$ located at $x_1,y_1$; and a second point vortex of strength $+\varGamma$ located at $x_2,y_2$. Refer to figure 5. Initially,

(4.1a–c)\begin{equation} y_p=l_p=0, \quad x_2=x_1, \quad y_2={-}y_1 \end{equation}

and, by symmetry, these conditions hold at all later times. The Lagrangian is (with $H_0=1$)

(4.2)\begin{align} &L[x_p,y_p,k_p,l_p,x_1,y_1,x_2,y_2] \nonumber\\ &\quad =\int \textrm{d}t \left[ \mathcal{A}_p ( k_p \dot{x}_p+l_p \dot{y}_p -\omega_r(k_p,l_p) ) +\varGamma ( {x}_1\dot{y}_1-{x}_2\dot{y}_2 ) -\frac{\varGamma^2}{2{\rm \pi}} \ln\vert \boldsymbol{x}_1-\boldsymbol{x}_2 \vert\right. \nonumber\\ &\qquad \left.+\frac{\mathcal{A}_p}{2{\rm \pi}} \left( -\varGamma \frac{(\boldsymbol{x}_1-\boldsymbol{x}_p) \times \boldsymbol{k}_p} {\vert \boldsymbol{x}_1-\boldsymbol{x}_p \vert ^2} +\varGamma \frac{(\boldsymbol{x}_2-\boldsymbol{x}_p) \times \boldsymbol{k}_p} {\vert \boldsymbol{x}_2-\boldsymbol{x}_p \vert ^2} \right) \right]. \end{align}

We vary all the dependent variables, and then apply the symmetry conditions (4.1a–c) to obtain a closed set of four equations. (It is not legitimate to apply the symmetry condition before taking the variations.) The reduced set of equations is

(4.3)\begin{gather} \delta k_p: \ \dot{x}_p= c_g(k_p,0)-\frac{\varGamma y_1}{{\rm \pi} d^2}, \end{gather}

(4.4)\begin{gather}\delta x_p: \ \dot{k}_p=\frac{2\varGamma k_p}{\rm \pi}\frac{(x_1-x_p)y_1}{d^4} , \end{gather}

(4.5)\begin{gather}\delta x_1: \ \dot{y}_1=\frac{\mathcal{A}_p k_p}{\rm \pi} \frac{(x_1-x_p)y_1}{d^4}, \end{gather}

(4.6)\begin{gather}\delta y_1: \ \dot{x}_1={-}\frac{\varGamma}{4{\rm \pi} y_1} +\frac{\mathcal{A}_p k_p}{2{\rm \pi}} \frac{(x_1-x_p)^2-y_1^2}{d^4}, \end{gather}

where $c_g$ is the $x$ component of the group velocity and

(4.7)\begin{equation} d^2 \equiv (x_1-x_p)^2+y_1^2 \end{equation}

is the squared distance between the wave packet and either vortex. Because of the symmetry conditions (4.1a–c), we do not need the evolution equations for $y_p$, $l_p$, $x_2$ and $y_2$.

Figure 5. (a) A right-moving wave packet (in black), with its wave vector denoted by the straight arrow and its associated dipolar flow indicated by circles, collides with a left-moving pair of counter-rotating vortices (in red). As the wave packet approaches the vortices, the flow induced by the vortex pair squeezes the wave packet in the $x$ direction, stretching its wave vector. The dipolar flow induced by the wave packet pushes the vortices apart (b). After passage of the wave packet (c) the solution ‘unwinds’, and all three particles return to their original configurations. (d) Partition of the energy (4.15) into wave energy $H_w$ (the first term on the left-hand side of (4.15)), vortex energy $H_v$ (the second term) and interaction energy $H_{int}$. Energies $H_w$ and $H_v$ increase during the interaction, while $H_{int}$ decreases.

Equations (4.3)–(4.6) conserve energy in the form

(4.8)\begin{equation} E=\omega_r(k_p) \mathcal{A}_p+\frac{\varGamma^2}{2{\rm \pi}} \ln y_1 -\frac{\mathcal{A}_p\varGamma}{\rm \pi} \frac{y_1 k_p}{d^2} \end{equation}

and momentum in the form

(4.9)\begin{equation} \mathcal{M}=\mathcal{A}_p k_p -2\varGamma y_1. \end{equation}

The angular momentum vanishes. Defining

(4.10)\begin{equation} X(t)\equiv x_p(t)-x_1(t), \end{equation}

we rewrite (4.3)–(4.6) as three equations, namely

(4.11)\begin{gather} \dot{X}= c_g(k_p)+\frac{\varGamma}{4{\rm \pi} d^2 y_1}(X^2-3y_1^2) -\frac{\mathcal{A}_p k_p}{2{\rm \pi} d^4} (X^2-y_1^2), \end{gather}

(4.12)\begin{gather}\dot{k}_p={-}\frac{2\varGamma k_p}{{\rm \pi} d^4}Xy_1, \end{gather}

(4.13)\begin{gather}\dot{y}_1={-}\frac{\mathcal{A}_p k_p}{{\rm \pi} d^4} Xy_1, \end{gather}

in the three unknowns $k_p$, $y_1$ and $X$, where now $d^2= X^2+y_1^2$. The two conserved quantities, (4.8) and (4.9), make this an integrable system. Eliminating $y_1$ between (4.8) and (4.9), we obtain an expression for the energy in terms of $X$ and $k_p$. The motion is confined to curves of constant $E(X,k_p)$. We can determine the solution by considering $E(X,k_p)$ or, even more conveniently, by considering

(4.14)\begin{equation} E(k_p,d^2)=\omega_r(k_p) \mathcal{A}_p +\frac{\varGamma^2}{2{\rm \pi}} \ln (\mathcal{A}_pk_p -\mathcal{M}) -\frac{\mathcal{A}_pk_p}{2{\rm \pi}} \frac{(\mathcal{A}_pk_p-\mathcal{M})}{d^2}, \end{equation}

in which we have dropped additive constants. Only the last term in (4.14) involves $d^2$.

Consider a gravity wave packet, initially at $X=-\infty$ with $k_p>0$, approaching the vortex pair from the left, as shown in figure 5(a). While the wave packet is still far from the vortex pair ($d^2$ very large) the last term in (4.14) is negligible. According to (4.13), $k_p$ increases with time on $X<0$. This increase in $k_p$ occurs because the velocity field associated with the vortices squeezes the wave packet in the $x$ direction. Since $c_g (k_p)>0$ the wave energy $\omega _r \mathcal {A}_p$ and the vortex-interaction energy – the middle term in (4.14) – both increase with $k_p$. The increase in the latter corresponds to the two vortices being pushed apart by the velocity field associated with the dipole. The increase in these two terms must be balanced by the last term in (4.14), which represents the energy stored in the superposed velocity fields of the wave packets and vortices. These superposed fields tend to cancel as the wave packet approaches the vortex pair. The value of $k_p$ reaches its maximum at $X=0$, where

(4.15)\begin{equation} d^2=y_1^2=\left( \frac{\mathcal{A}_pk_p -\mathcal{M}}{2\varGamma} \right)^2. \end{equation}

Substituting (4.15) into (4.14), we obtain an equation for this maximum value of $k_p$. After passing $X=0$, the solution ‘unwinds’, and $k_p$ returns to its original value as $X \rightarrow \infty$. The numerical solution shown in figure 5 confirms this analysis.

4.1. Wave-packet-induced drift

If $\varGamma \ll \mathcal{A}_{p}$ and we define $\tilde {X}= -X$, then (4.9)–(4.13) imply

(4.16)\begin{gather} \dot{\tilde{X}} ={-}c_g(k_p) +\frac{\mathcal{A}_pk_p}{2{\rm \pi} d^4} (\tilde{X}^2-y_1^2), \end{gather}

(4.17)\begin{gather}\dot{k}_p = 0, \end{gather}

(4.18)\begin{gather}\dot{y}_1 = \frac{\mathcal{A}_p k_p}{{\rm \pi} d^4}\tilde{X}y_1. \end{gather}

Hence $k_p$ and $c_g$ are constants. The governing equations (with $k_{p} \mathcal{A}_{p}=2{\rm \pi} a^2$ and $c_{g0} = 1$) reduce to

(4.19)\begin{gather} \dot{\tilde{X}} ={-}c_{g0} +\frac{a^2}{d^4} (\tilde{X}^2-y_1^2), \end{gather}

(4.20)\begin{gather}\dot{y}_1 = \frac{2 a^2}{d^4}\tilde{X}y_1. \end{gather}

In this limit, the point vortices are passive; their motion is the same as that of fluid particles in the presence of a uniformly translating cylinder. This problem was examined by Maxwell (Reference Maxwell1870; see also Morton Reference Morton1913; Darwin Reference Darwin1953). In the reference frame moving with the wave packet, the stream function is an integral of motion. Hence

(4.21)\begin{equation} Y_0=y_1\left(1-\frac{a^2}{d^2}\right) \end{equation}

is constant. Define $\dot {\mathcal {X}} = \dot {\tilde {X}}-c_{g0}$ and $\theta = \tan ^{-1}(y_1/\tilde {X})$. Using (4.21) and following Maxwell (Reference Maxwell1870) and Darwin (Reference Darwin1953) we find that

(4.22)\begin{equation} \mathcal{X} = \int \frac{a^2\cos 2\theta }{\sqrt{Y_0^2+4a^2\sin^2\theta}}\, \textrm{d} \theta. \end{equation}

Let $\cos \theta =- \text {sn} (\nu )$ with the suppressed modulus of the Jacobi elliptic function understood to be

(4.23)\begin{equation} m=\kappa^2=\frac{4a^2}{Y_0^2+4a^2}. \end{equation}

It is tedious but straightforward to show that

(4.24)\begin{equation} \mathcal{X}=\frac{a}{\kappa}\left(\left(1-\frac{\kappa^2}{2}\right) \nu -E(\nu)\right). \end{equation}

Similar expressions may be found for $y_1(\nu )$ and $t(\nu )$ but fall outside the scope of our discussion (Darwin Reference Darwin1953).

From (4.24), the total drift in the $x$ direction is

(4.25)\begin{equation} \Delta \mathcal{X} = \frac{2a}{\kappa} \left(\left(1-\frac{\kappa^2}{2}\right) K-E\right), \end{equation}

where $E$ and $K$ are the complete elliptical integrals of the first and second kind, respectively.

The drift volume, $D$, is defined as

(4.26)\begin{equation} D=\int_{-\infty}^{\infty}\mathcal{X} \, \textrm{d} Y_0 = {\rm \pi}a^2. \end{equation}

Note that the connection between Stokes drift (specifically the motion in the vertical plane) and Darwin drift has been examined by Eames & McIntyre (Reference Eames and McIntyre1999).

5. One vortex, $N$ wave packets

We now search for simple harmonic motion in a system comprising one vortex and $N>1$ wave packets (see also the related discussion of a ring of geostrophic vortices in Morikawa & Swenson (Reference Morikawa and Swenson1971)). The single vortex of strength $\varGamma = 2{\rm \pi}$ remains stationary at $\boldsymbol {x}=\boldsymbol {0}$. The $N$ wave packets at $\boldsymbol {x}_p=(x_p,y_p)$ have equal actions $\mathcal {A}_p=\mathcal {A}_0$, and wave vectors of equal magnitude $|\boldsymbol {k}_p|=\kappa$. They lie symmetrically on the circle $|\boldsymbol {x}_p|=\chi$. Both $\kappa$ and $\chi$ are constants. We take the depth parameter to be unity, $H_0=1$.

A vortex at $\boldsymbol {x}$ induces the velocity field

(5.1)\begin{equation} \boldsymbol{U}_m(\boldsymbol{x}', \boldsymbol{x} )= \frac{(y-y',x'-x)}{(x-x')^2+(y-y')^2} \end{equation}

at $\boldsymbol {x}'$. Hence,

(5.2)\begin{equation} \boldsymbol{U}_m(\boldsymbol{x}_p, \boldsymbol{0} ) = \frac{({-}y_p,x_p)}{\chi^2}. \end{equation}

We also need

(5.3)\begin{equation} \boldsymbol{\nabla} U_m(\boldsymbol{x}_p, \boldsymbol{0} ) = \frac{1}{\chi^4} ( 2x_p y_p , x_p^2- y_p^2 ) \end{equation}

and

(5.4)\begin{equation} \boldsymbol{\nabla} V_m(\boldsymbol{x}_p, \boldsymbol{0} ) = \frac{1}{\chi^4} ({-}x_p^2+ y_p^2 , -2x_p y_p ). \end{equation}

The wave packet at $\boldsymbol {x}_p$ induces the velocity

(5.5)\begin{equation} \boldsymbol{U}_d={-} \frac{\mathcal{A}_0}{2{\rm \pi}\chi^4} ( 2l_px_py_p+k_p(x_p^2-y_p^2), \ 2k_px_py_p-l_p(x_p^2-y_p^2) ) \end{equation}

at the vortex.

The equations of motion reduce to

(5.6)\begin{gather} \dot{\boldsymbol{x}}_p=\frac{1}{2}\sqrt{\frac{g}{\kappa^3}}(k_p,l_p)+\frac{({-}y_p,x_p)}{\chi^2}, \end{gather}

(5.7)\begin{gather}\chi^4\dot{\boldsymbol{k}}_p=({-}2k_px_py_p+l_p(x_p^2-y_p^2),2l_px_py_p+k_p(x_p^2-y_p^2)) \end{gather}

and

(5.8)\begin{equation} \sum_p \mathcal{A}_0({-}2l_px_py_p-k_p(x_p^2-y_p^2),-2k_px_py_p+l_p(x_p^2-y_p^2))=\boldsymbol{0}. \end{equation}

The last equation is the condition that the vortex remains stationary.

We look for solutions exhibiting simple harmonic motion with $\boldsymbol {k}_p=\kappa (\sin \theta _p,-\cos \theta _p)$ and $\boldsymbol {x}_p=\chi (\cos \theta _p, \sin \theta _p)$. This leads to the constraints

(5.9)\begin{gather} \kappa = \frac{\chi^2}{16}, \end{gather}

(5.10)\begin{gather}\frac{\textrm{d}\theta_{p}}{\textrm{d}t}={-} 1/\chi^2 \end{gather}

and

(5.11)\begin{equation} \mathcal{A}_0\sum (\cos \theta_p,\sin \theta_p) = \boldsymbol{0}. \end{equation}

We observe that solutions to this system are related to the $N$th roots of unity. This sets the initial phase of $\theta _p$. We find that

(5.12a,b)\begin{equation} \mathcal{M}=\boldsymbol{0}; \quad \mathcal{H} = \frac{3\chi }{16}\sum \frac{\mathcal{A}_0}{2{\rm \pi}}. \end{equation}

Figure 6 shows a solution of this type with four wave packets.

Figure 6. Trajectories exhibiting circular motion for a system with four wave packets and one vortex. The wave packets are shown in black, with their intensity increasing with time. The red circle represents the stationary vortex.

6. The Lagrangian motion of a particle near a periodic wave packet

We now consider the motion of a very weak point vortex near a periodic array of wave packets. In the limit of vanishing circulation, the vortex acts as a passive tracer, and sets the foundation for the stability analysis performed in § 7. This is a generalization of the motion considered by Maxwell (Reference Maxwell1870) and Darwin (Reference Darwin1953). Unlike the analysis there, closed-form solutions for the motion of a particle are not found, but asymptotic analysis reveals interesting features of the induced flow.

The system we now consider is unbounded in $x$ and $y$, and infinitely periodic in the $x$ direction. Initially, the wave packet propagates along the $x$ axis. Each vortex or wave packet at $(x,y)$ sees its images at $(x+2{\rm \pi} n,y)$, where $n$ is any integer. Again we take $H_0=1$.

The governing equations are (2.49)–(2.51) with $\boldsymbol {U}_m$ and $\boldsymbol {U}_d$ now calculated from the stream functions

(6.1)\begin{equation} \psi_m(\boldsymbol{x},\boldsymbol{x}_i) = \frac{1}{4{\rm \pi}} \sum_n \ln [(x-x_i+2{\rm \pi} n)^2 + (y-y_i)^2] \end{equation}

and

(6.2)\begin{equation} \psi_d (\boldsymbol{x},\boldsymbol{x}_p,\boldsymbol{k}_p) = \frac{1}{2{\rm \pi}}\sum_n \frac{(\boldsymbol{x}-\boldsymbol{x}_p+2{\rm \pi}(n,0))\times \boldsymbol{k}_p}{|\boldsymbol{x}_p-\boldsymbol{x}+2{\rm \pi}(n,0)|^2}. \end{equation}

6.1. Limit of weak point vortices

We begin by considering the wave-packet-induced motion of the vortices when $|\varGamma _i| \ll 1$. As in § 4 – see particularly § 4.1 – this motion takes a non-trivial form.

We consider a wave packet travelling along the $x$ axis so that $l=y_p=0$. To $O(1)$, (2.49)–(2.51) reduce to

(6.3)\begin{gather} \dot{x}_p = \frac{1}{2}\sqrt{\frac{g}{|\boldsymbol{k}_p|}} \frac{\boldsymbol{k}_p}{|\boldsymbol{k}_p|}, \end{gather}

(6.4)\begin{gather}\dot{\boldsymbol{k}}_p=0, \end{gather}

(6.5)\begin{gather}\dot{\boldsymbol{x}}_i = \boldsymbol{U}_d(\boldsymbol{x}_i, \boldsymbol{x}_p, \boldsymbol{k}_p ), \end{gather}

where $\boldsymbol {U}_d$ is computed from (6.2). As we are considering one (weak) vortex per period, $i=1$ and we take $y_1=y$.

It is convenient to define

(6.6)\begin{equation} \chi = x_1 -x_p. \end{equation}

The governing equations become

(6.7)\begin{gather} \dot{\chi} ={-}c_{g0}+\mu \sum_{n={-}\infty}^{\infty} \frac{(\chi-2{\rm \pi} n)^2-y^2}{((\chi-2{\rm \pi} n)^2+y^2)^2}, \end{gather}

(6.8)\begin{gather}\dot{y} = 2\mu \sum_{n={-}\infty}^{\infty} \frac{(\chi-2{\rm \pi} n)y}{((\chi-2{\rm \pi} n)^2+y^2)^2}, \end{gather}

where we have set $\mu \equiv k_0 \mathcal {A}_p/2{\rm \pi}$ and $c_{g0}= 1/2\sqrt {g/k_0}$.

Defining $z = \chi +\mathrm {i} y$, we have

(6.9)\begin{equation} \dot{z} ={-}c_{g0}+\mu\sum_{n={-}\infty}^{\infty} \frac{(z-2{\rm \pi} n)^2}{(z-2{\rm \pi} n)^2(z^*-2{\rm \pi} n)^2}={-}c_{g0}+\mu\sum_{n={-}\infty}^{\infty} \frac{1}{(z^*-2{\rm \pi} n)^2}. \end{equation}

Define the complex-valued velocity potential $w=\phi +\mathrm {i} \psi$, where $(w_z)^*\equiv \dot {z}$. This implies (Lamb Reference Lamb1932, § 64)

(6.10)\begin{equation} w={-}c_{g0}\chi -\mu \sum_{n={-}\infty}^{\infty}\frac{1}{2{\rm \pi} n-z} ={-}c_{g0}\chi - \mu \cot\frac{z}{2}, \end{equation}

so that

(6.11)\begin{equation} \phi ={-} c_{g0}\chi + \mu\frac{\sin \chi }{ \cos \chi - \cosh y} \end{equation}

and

(6.12)\begin{equation} \psi ={-} c_{g0}y- \mu \frac{\sinh y }{\cos \chi - \cosh y}. \end{equation}

From the relation $(w_z)^*=\dot {z}$, we find

(6.13)\begin{gather} \dot{\chi} ={-}c_{g0}+ \frac{\mu}{2} \frac{1-\cos \chi \cosh y}{(\cos \chi - \cosh y)^2}, \end{gather}

(6.14)\begin{gather}\dot{y} =\frac{\mu}{2} \frac{\sin \chi \sinh y}{(\cos \chi - \cosh y)^2}. \end{gather}

The stream function is a material contour. As in § 4.1, this provides an additional conserved quantity. When $(\chi _i,y_i)$ are small – vortex $i$ very close to the wave packet – our equations reduce to the equations of Maxwell discussed in § 4.1.

Although we have been unable to find closed-form solutions, asymptotic analysis reveals interesting properties of this system. We rewrite the transcendental pieces of (6.13) and (6.14) as

(6.15)\begin{equation} \frac{\sinh y\sin \chi}{(\cosh y - \cos \chi)^2} = 2 \sum_{n=1}^{\infty} n\,\textrm{e}^{{-}ny}\sin n\chi \end{equation}

and

(6.16)\begin{equation} \frac{1-\cosh y\cos \chi}{(\cosh y - \cos \chi)^2} ={-}2 \sum_{n=1}^{\infty} n\,\textrm{e}^{{-}n y}\cos n\chi. \end{equation}

This expansion assumes $y > 0$; a similar formula holds for $y<0$. The equations of motion may then be written as

(6.17)\begin{gather} \dot{\chi}={-}c_{g0}-\mu \sum_{n=1}^{\infty} n\,\textrm{e}^{{-}ny}\cos n\chi, \end{gather}

(6.18)\begin{gather}\dot{y}= \mu \sum_{n=1}^{\infty} n \,\textrm{e}^{{-}ny}\sin n\chi, \end{gather}

or, more compactly, as

(6.19)\begin{equation} \dot{z} ={-}c_{g0}-\mu \sum_{n=1}^{\infty} n \,\textrm{e}^{- \mathrm{i} n z^*}. \end{equation}

Let $b>0$ be the initial vertical coordinate of the particle and (for reasons that become clear in § 7) let ${\rm \pi} /2$ be its horizontal coordinate. Let the initial location of the wave packet be $({\rm \pi} ,0)$.

As the wave packets and vortex are sufficiently far apart, a natural small parameter arises, namely

(6.20)\begin{equation} \epsilon \equiv \textrm{e}^{{-}b}\ll 1. \end{equation}

We expand $z$ as

(6.21)\begin{equation} z = z_{0}+\epsilon z_{1}+\epsilon^2 z_{2}+\cdots. \end{equation}

To $O(\epsilon ^2)$, we find

(6.22)\begin{equation} \dot{z} =\dot{z}_{0}+\epsilon \dot{z}_{1}+\epsilon^2 \dot{z}_{2}={-}c_{g0} - \mu (\epsilon \,\textrm{e}^{-\mathrm{i} \chi_0} +\epsilon^2(2\textrm{e}^{{-}2\mathrm{i} \chi_0}-\mathrm{i} z_{1}^*\,\textrm{e}^{-\mathrm{i} \chi_0})). \end{equation}

We solve this system iteratively. To lowest order,

(6.23)\begin{equation} z_0 = z_0^0-c_{g0}t, \end{equation}

where $z_0^0=-{\rm \pi} /2+ \mathrm {i} b$. This implies

(6.24)\begin{equation} z_1 ={-} \frac{\mathrm{i} \mu}{c_{g0}}(\textrm{e}^{-\mathrm{i} \chi_0}-\mathrm{i}) ={-}\frac{\mu}{c_{g0}}(\textrm{e}^{\mathrm{i} \theta}-1), \end{equation}

where $\chi _0 = -{\rm \pi} /2-c_{g0}t$ and we define

(6.25)\begin{equation} \theta \equiv c_{g0}t. \end{equation}

Thus $z_1$ takes the form of a simple harmonic oscillator. The constant is taken to ensure that $z_1=0$ at $t=0$.

At second order we find

(6.26)\begin{equation} \dot{z}_2 = \mu (2 \textrm{e}^{2\mathrm{i} \theta} - z_1^* \,\textrm{e}^{\mathrm{i} \theta}), \end{equation}

so that upon substitution of the result for $z_1$, we find

(6.27)\begin{equation} z_2 = \frac{\mu^2}{c_{g0}}t-\mathrm{i} \frac{\mu}{ c_{g0}}(\textrm{e}^{2\mathrm{i}\theta }-1) + \mathrm{i} \frac{\mu^2}{ c_{g0}^2} (\textrm{e}^{\mathrm{i} \theta}-1). \end{equation}

At this order there is a mean drift in the direction of wave propagation. This drift has important implications for the stability of vortex streets, as discussed in § 7.

Figure 7 compares our approximate analytical solutions with numerical integrations of the full equations for $\epsilon = 1/10$ (figure 7a,c,e) and $\epsilon =1/3$ (figure 7b,d,f), where we have taken $\mathcal {A}_p=k_0=1$. We see that the asymptotic theory works relatively well for small values of $\epsilon$.

Figure 7. A comparison of the numerical integration of the full equations of motion (black lines) and the second-order asymptotic solutions (dashed red lines). (a,c,e) Results when $\epsilon =1/10$. (b,d,f) Results for $\epsilon =1/3$. (a–d) Plots of $y_1$ and $x_1$, the vertical and horizontal motion of the vortex, as a function of time. (e,f) The behaviour of the wave packet. We see that the asymptotic theory describes the numerical results well for the case of $\epsilon =1/10$ but begins to break down at longer times when $\epsilon = 1/3$.

7. The stability of a symmetric vortex street in the presence of a wave packet

Numerical solutions (not here described in detail) suggest that there are cases in which the wave packets organize the vortices into patterns. As a first step in understanding this phenomenon we study the stability of a vortex street in the presence of a single wave packet (per period) in the semi-periodic domain. This system conserves energy and momentum, but angular momentum is not conserved because the periodicity breaks rotational symmetry.

Within the (periodic) domain, we have four vortices, arranged symmetrically about $y=0$, with $y=\pm b$, and spaced ${\rm \pi}$ apart in the $x$ direction. Refer to figure 8. This is the minimal system to illustrate the instability of a periodic vortex street (Domm Reference Domm1956). The stability of vortex streets was first considered by von Kármán (Reference von Kármán1911), and discussed in detail by Lamb (Reference Lamb1932, § 156). Domm (Reference Domm1956) examined the nonlinear stability, reducing the system of four vortices to two dependent variables. The symmetric vortex street proved to be linearly unstable in all of parameter space, whereas the asymmetric vortex street is linearly stable at a single value of the ratio of vertical to horizontal spacing of the vortices. However, even the linearly stable asymmetric vortex street is nonlinearly unstable (Domm Reference Domm1956).

Figure 8. The vortex street configuration considered in this section. The vortices in the top row have strength $\varGamma$ while those in the bottom row have strength $-\varGamma$. The domain is periodic in the $x$ direction, and there is a wave packet travelling along the $x$ axis. We study the stability of this vortex street in the presence of the wave packet.

The equations of motion are (2.49)–(2.51) with the stream functions given by (6.1) and (6.2). We take $\varGamma _i = \epsilon ^2 \gamma$ so that $\varGamma _1 = \varGamma _2 = \epsilon ^2\gamma$ and $\varGamma _3 = \varGamma _4 = - \epsilon ^2\gamma$. We seek equations of motion valid to $O(\epsilon ^2)$. We begin by expanding the variables describing the wave packet as

(7.1a,b)\begin{gather} x_p = x_{p0}+c_{g0}t+\epsilon^2 x_{p2},\quad y_p = 0, \end{gather}

(7.2a,b)\begin{gather}k_p = k_{p0},\quad l_p=0. \end{gather}

If there are no second-order corrections to the wavenumbers initially, and if the momentum is conserved, then the second-order wavenumbers must remain zero for all time.

The governing second-order equation for the horizontal motion of the wave packet is

(7.3)\begin{equation} \dot{x}_{p2} ={-}\frac{1}{4{\rm \pi}} \sum_i \gamma_i \frac{\sinh y_{i0}}{\cos \chi_{i0}-\cosh y_{i0}}. \end{equation}

Expanding hyperbolic functions, we find that

(7.4)\begin{equation} \dot{x}_{p2} = \frac{\gamma}{ {\rm \pi}}+O(\epsilon), \end{equation}

where $\gamma \equiv \gamma _1$. Equation (7.4) implies a correction to the wave packet speed, due to the velocity induced by the vortex street, which has the magnitude of the mean velocity in the plane of symmetry of the vortex street (Lamb Reference Lamb1932, § 156). Expanding

(7.5)\begin{equation} \dot{y}_{p2} = \frac{1}{4{\rm \pi}}\sum_i \frac{\sin \chi_{i0} }{\cos \chi_{i0}-\cosh y_{i0}} = 0+O(\epsilon), \end{equation}

we find that there is no correction to the vertical speed of the wave packet.

We now solve for the motion of the point vortices. It will be a combination of the motion induced by the wave packet (as calculated in § 6) and the uniform self-advection of the vortex street. Lamb (Reference Lamb1932, § 156) found that the self-advection of the vortex street leads to a uniform translation of the vortices with speed

(7.6)\begin{equation} \epsilon^2 \frac{\gamma}{2{\rm \pi}} \coth b = \epsilon^2 \frac{\gamma}{2{\rm \pi}} + O(\epsilon^3). \end{equation}

Then, from (6.27), we obtain the second-order behaviour of the vortices as

(7.7)\begin{equation} z_{2i} = \left(\frac{\gamma}{2{\rm \pi}} + \frac{\mu^2}{c_{g0}}\right)t-\frac{\mu }{ c_{g0}}(\textrm{e}^{{-}2\mathrm{i}\chi_{i0} }-1) + \frac{ \mu^2}{c_{g0}^2}( \textrm{e}^{-\mathrm{i} \chi_{i0}}-1), \end{equation}

where, from our initial conditions, $\chi _{10} =\chi _{30}= -{\rm \pi} /2 -c_{g0}t$ and $\chi _{20}=\chi _{40}={\rm \pi} /2 - c_{g0}t$. Note the presence of two mean flows. One is induced by the wave packet; the other represents the self-advection of the vortex street. These two competing mean flows are now shown to have implications for the stability of the vortex street.

7.1. Linear stability analysis

Now we perturb our system to examine its linear stability. We expand the vortex locations as

(7.8)\begin{gather} z_i = z_{i0}+\epsilon z_{i1}+\epsilon^2 z_{i2} + \delta z_{\delta i}, \end{gather}

(7.9a–d)\begin{gather}x_p = x_{p0}+\epsilon^2 x_{p2}+\delta x_{p\delta}, \quad y_p=\delta y_{p\delta},\quad k_p = k_{p0}+\epsilon^2\delta k_{p\delta}, \quad l_p =\epsilon^2\delta l_{p\delta}, \end{gather}

where $\delta$ is the amplitude of the perturbations. The goal is to expand the equations of motion to $O(\epsilon ^2 \delta )$. Starting with the motion of the wave packet, and using the results found in the previous subsection, we find that at $O(\delta )$

(7.10)\begin{gather} \dot{x}_{p\delta} = c_{gx\delta}, \end{gather}

(7.11)\begin{gather}\dot{y}_{p\delta} = c_{gy\delta}, \end{gather}

where $(c_{gx\delta },c_{gy\delta })$ are the $O(\delta )$ expansions of the group velocity. The wavenumbers evolve according to

(7.12)\begin{equation} \dot{k}_{p\delta}=\dot{l}_{p\delta}=0. \end{equation}

If our initial conditions are such that these perturbation wavenumbers vanish, they will vanish for all time. This implies that $(x_{p\delta },y_{p\delta })$ are constants which we take to be zero.

The non-trivial part of the analysis comes from the evolution of the vortices. We must solve for $\boldsymbol {U}_m$ and $\boldsymbol {U}_d$ to $O(\epsilon ^2 \delta )$. At this order, the vortex-induced flow is the same as it would be in the absence of the wave packet; hence, from Lamb (Reference Lamb1932), we have the $O(\epsilon ^2 \delta )$ result:

(7.13)\begin{gather} \boldsymbol{U}_{m}(z_1) ={-}\mathrm{i}\frac{\gamma}{8{\rm \pi}}(- z_{\delta 1}+z_{\delta 2}), \end{gather}

(7.14)\begin{gather}\boldsymbol{U}_{m}(z_2) = \mathrm{i}\frac{\gamma}{8{\rm \pi}}(- z_{\delta 1}+z_{\delta 2}). \end{gather}

Also, at this order $\boldsymbol {U}_m(z_1)$ and $\boldsymbol {U}_m(z_2)$ only depend on vortices 1 and 2; hence we omit $\boldsymbol {U}_m(z_3)$ and $\boldsymbol {U}_m(z_4)$ for clarity of presentation. That is, we only need to consider the evolution of vortices 1 and 2, or equivalently vortices 3 and 4.

The $O(\delta )$ contributions to $\boldsymbol {U}_d(z_{1,2})$ are found by substituting the expansion given by (7.8) into (6.19) so that

(7.15)\begin{equation} \boldsymbol{U}_d(z_i) = \mu \epsilon \delta \,\textrm{e}^{{-}2\mathrm{i} \chi_{i0}}z_{\delta i}^* (\mathrm{i} \textrm{e}^{\mathrm{i} \chi_{i0}}+\epsilon (4\mathrm{i} +\textrm{e}^{\mathrm{i} \chi_{i0}}z_{i1}^*)). \end{equation}

To solve this system of equations, we must also expand the perturbations in $\epsilon$, so that

(7.16)\begin{equation} z_{\delta i} = z_{\delta i0} +\epsilon z_{\delta i1}+\epsilon^2 z_{\delta i2}. \end{equation}

We assume the time dependence of the order-zero terms goes like $\textrm {e}^{\varLambda t}$, where $\varLambda =\epsilon ^2\lambda$ (which may be inferred from the classical stability analysis, which implies that the growth rates are proportional to the strength of the vortices; see Lamb (Reference Lamb1932)). Additionally, we assume that the first-order terms have fast oscillations, so that their time derivatives have the same order as the original term. We need not solve explicitly for the second-order terms to conduct the stability analysis.

The restrictions on the first-order terms are found to imply

(7.17)\begin{equation} z_{\delta 1 1} =\mathrm{i}\frac{\mu}{2c_{g0}}z_{\delta 10}^*(\textrm{e}^{\mathrm{i} \theta}-1) \end{equation}

and

(7.18)\begin{equation} z_{\delta 2 1} ={-}\mathrm{i}\frac{\mu}{2c_{g0}}z_{\delta 20}^*(\textrm{e}^{\mathrm{i} \theta}-1). \end{equation}

We now have sufficient information to solve for the stability of the vortex street configuration, which will be dictated by the evolution of $\{z_{\delta 10},z_{\delta 20}\}$. Substituting the lower-order expansions into $\boldsymbol {U}_d$ and $\boldsymbol {U}_m$, we find that the dynamics of $\{z_{\delta 10},z_{\delta 20}\}$ is given by the phase-averaged equations of motion.

The evolution equations are given by

(7.19)\begin{gather} \dot{z}_{\delta 1 0} = \epsilon^2\left(-\frac{\mathrm{i} \nu}{2}(z_{\delta 10}^*-z_{\delta 20}^*)+\frac{\mathrm{i} \sigma }{2}(z_{\delta 10}-z_{\delta 10}^*) \right), \end{gather}

(7.20)\begin{gather}\dot{z}_{\delta 2 0} =\epsilon^2\left(\frac{\mathrm{i} \nu}{2}(z_{\delta 10}^*-z_{\delta 20}^*)+\frac{\mathrm{i} \sigma }{2}(z_{\delta 20}-z_{\delta 20}^*) \right), \end{gather}

where

(7.21a,b)\begin{equation} \sigma = \frac{\mu^2}{2c_{g0}}; \quad \nu = \frac{\gamma}{4{\rm \pi}}. \end{equation}

Define $z_{\delta 10} = \alpha \,\textrm {e}^{ \epsilon ^2 \lambda t},z_{\delta 20} = \beta \,\textrm {e}^{ \epsilon ^2 \lambda t}$. Then, (7.19)–(7.20), together with their complex conjugates, imply the following eigenvalue problem:

(7.22)\begin{equation} \begin{pmatrix} \sigma-2\mathrm{i}\lambda & 0 & -\nu-\sigma & \nu \\ 0 & \sigma-2\mathrm{i}\lambda & \nu & -\nu-\sigma \\ \nu+\sigma & -\nu & -\sigma-2\mathrm{i}\lambda & 0 \\ -\nu & \nu+\sigma & 0 & -\sigma-2\mathrm{i}\lambda \end{pmatrix} \begin{pmatrix} \alpha\\ \beta\\ \alpha^*\\ \beta^* \end{pmatrix}=\boldsymbol{0}. \end{equation}

The eigenvalues $\lambda$ are found to be

(7.23)\begin{equation} \lambda^2 =\nu(\nu+\sigma). \end{equation}

In the limit that $\sigma =0$, we recover the result of Lamb (Reference Lamb1932, § 156) that the symmetric vortex street is linearly unstable. It follows from (7.23) that the system is stable when

(7.24)\begin{equation} \nu(\nu+\sigma)<0. \end{equation}

We note that $\nu$ can be positive or negative depending on the sign of $\gamma$, so that the inequality may hold only if the sign of $\nu$ and $\sigma$ are different. Physically, the stability depends on the sign of the induced motion of the vortices. If the drift induced by the wave packets is larger than the self-advection of the vortices, and $\nu <0$, the system remains linearly stable.

7.2. Long-time behaviour

In the previous subsection, we analysed the linear stability of the vortex street to perturbations. However, as was shown by Domm (Reference Domm1956), even a linearly stable vortex street might be nonlinearly unstable. This may be examined analytically, but the algebra becomes considerably involved, so we instead perform a numerical investigation. We integrate our equations of motion using a fourth-order Runge–Kutta scheme, ensuring that the Hamiltonian and momentum are conserved. We also increase the number of adjacent periodic extensions until convergence is found (here we take our total domain to be of length $5001\times 2{\rm \pi}$). We take $\epsilon = 1/3$, $k_0 = 1$, $\mathcal {A}_p=1$, $\varGamma = \epsilon ^2$, and we integrate the system for $3\times 10^4$ s. This is shown in figure 9(a), where it is seen that the vortex street remains bound to a small neighbourhood over a long integration time. In the absence of the wave packet $(\mathcal {A}_p=0)$, shown in figure 9(b), the system is unstable, and the vortex street eventually dissolves, analogous to the behaviour found for two vortex pairs (Love Reference Love1893; Tophøj & Aref Reference Tophøj and Aref2013).

Figure 9. (a) Stability of periodic (in the $x$ direction) vortex street, with initial conditions depicted in figure 8. The wave packet motion is shown in blue, while the vortices are shown in black ($\varGamma _i<0$) and red $(\varGamma _i>0)$. The vortex street remains bound to a relatively small neighbourhood over this long-time integration. (b) The same initial configuration of the vortices with no wave packet, and we see the system is unstable, with vortices propagating far away from their initial locations.