Hostname: page-component-cd9895bd7-8ctnn Total loading time: 0 Render date: 2024-12-14T14:53:19.444Z Has data issue: false hasContentIssue false

Particle description of the interaction between wave packets and point vortices

Published online by Cambridge University Press:  27 August 2021

Nick Pizzo*
Affiliation:
Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA92037, USA
Rick Salmon
Affiliation:
Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA92037, USA
*
Email address for correspondence: npizzo@ucsd.edu

Abstract

This paper explores an idealized model of the ocean surface in which widely separated surface-wave packets and point vortices interact in two horizontal dimensions. We start with a Lagrangian which, in its general form, depends on the fields of wave action, wave phase, stream function and two additional fields that label and track the vertical component of vorticity. By assuming that the wave action and vorticity are confined to infinitesimally small, widely separated regions of the flow, we obtain model equations that are analogous to, but significantly more general than, the familiar system consisting solely of point vortices. We analyse stable and unstable harmonic solutions, solutions in which wave packets eventually coincide with point vortices (violating our assumptions), and solutions in which the wave vector eventually blows up. Additionally, we show that a wave packet induces a net drift on a passive vortex in the direction of wave propagation which is equivalent to Darwin drift. Generalizing our analysis to many wave packets and vortices, we examine the influence of wave packets on an otherwise unstable vortex street and show analytically, according to linear stability analysis, that the wave-packet-induced drift can stabilize the vortex street. The system is then numerically integrated for long times and an example is shown in which the configuration remains stable, which may be particularly relevant for the upper ocean.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2021. Published by Cambridge University Press

1. Introduction

This paper explores an idealized model of the ocean surface in which widely separated surface-wave packets and point vortices interact in two horizontal dimensions. Each wave packet $p$ is defined by its location $\boldsymbol {x}_p(t)$, its wave action $\mathcal {A}_p$ and its wave vector $\boldsymbol {k}_p(t)$. Each point vortex $i$ is defined by its location $\boldsymbol {x}_i(t)$ and its strength $\varGamma _i$. In reality, wave breaking converts wave action into vorticity, and vorticity is destroyed by viscosity. However, in this initial study we consider only the ideal case, in which $\mathcal {A}_p$ and $\varGamma _i$ are conserved.

The velocity field attached to the wave packets is dipolar; it is sometimes called a ‘Bretherton flow’. Wave packets advect the point vortices by their Bretherton flow. Point vortices advect wave packets and other point vortices, and change the wave vector of the wave packets by refraction. For simplicity, we omit the interactions between wave packets, which are expected to be weak.

In § 2 we derive the equations governing $\boldsymbol {x}_p(t)$, $\boldsymbol {k}_p(t)$ and $\boldsymbol {x}_i(t)$ from a Lagrangian which, in its general form, depends on the fields of wave action, wave phase, stream function and two additional fields that label and track the vertical component of vorticity. In our application, the Lagrangian couples Whitham's Lagrangian for surface waves to the Langrangian for two-dimensional, incompressible flow. Coupling is achieved by replacing the ‘mean velocity’ in the Doppler term of the dispersion relationship with the velocity field corresponding to the stream function of the vortical flow. We obtain our final equations by assuming that the wave action and vorticity are confined to infinitesimally small, widely separated regions of the flow. To leading order, each wave packet induces a dipolar horizontal flow, and each vortex patch induces a monopolar flow. In its general formulation (Salmon Reference Salmon2020), the method applies to any type of wave and any type of mean flow in two or three dimensions. It seems easier to apply than other, apparently equivalent methods that do not employ a Lagrangian.

In § 3 we consider the system consisting of a single wave packet and a single point vortex. We analyse harmonic solutions in which the two particles move in circular orbits. For these configurations, we show that solutions in which the vortex orbit lies outside the orbit of the wave packet are stable, whereas solutions in which the vortex orbit lies inside that of the wave packet are unstable. We also investigate solutions in which the vortex and wave packet eventually coincide, violating the assumption of our model, and solutions in which the wave vector grows without bound.

In § 4 we consider the case of a wave packet encountering a pair of counter-rotating point vortices. The highly symmetric arrangement permits thorough analysis, which is confirmed by numerical solutions. This solution is very similar to that discussed by Bühler & McIntyre (Reference Bühler and McIntyre2005) and invites a comparison with their method of analysis. We also show that, in the limit that the circulation of the vortices is much weaker than the wave action, the equations are equivalent to those diagnosing the motion of a particle in the presence of a uniformly translating cylinder. Following classical analysis (Maxwell Reference Maxwell1870; Darwin Reference Darwin1953) it is shown that the wave packet induces a net displacement on the vortices.

In § 5 we study a solution in which $N>1$ wave packets are equidistant from, and symmetrically arranged about, a single vortex. The wave packets circle the vortex at a uniform angular velocity, while the vortex remains stationary at the centre of the pattern.

In § 6 we generalize our system to be periodic in one dimension, and investigate the motion of a periodic array of weak point vortices in the presence of a periodic array of wave packets. We find asymptotic solutions in which the wave packets induce a net drift on the vortices.

In § 7 we use the analysis in § 6 to investigate the linear stability of a vortex street in the presence of a wave packet. We find that the wave packet can change the stability of the vortex street. Numerical analysis demonstrates that vortex streets can be stable for long times in the presence of a wave packet.

Section 8 concludes with an assessment of our results and their oceanographic implications.

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.30ad)\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.30ad). 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.1ac)\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.1ac) 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.1ac), 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$. (ad) 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.9ad)\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.

8. Conclusion

It has long been recognized in the oceanographic community that surface waves may freely exchange momentum and energy with underlying currents. The dynamics is governed by wave action conservation (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962; Whitham Reference Whitham1965), while the evolution of the phase obeys the equations of geometrical optics. These equations are valid for small-amplitude inviscid waves that are slowly varying. Despite the maturity of this theory (Phillips Reference Phillips1966), there are still many fundamental questions regarding wave–current interaction, including a need to better understand their two-way coupling (McWilliams, Restrepo & Lane Reference McWilliams, Restrepo and Lane2004). This is thrown into particularly stark relief by numerical models of climate, which are beginning to resolve the submesoscale (of the order of 1–10 km), where these interactions may be especially pronounced (McWilliams Reference McWilliams2016; Romero, Lenain & Melville Reference Romero, Lenain and Melville2017). At even smaller scales, wave breaking in deep water occurs for waves with finite crest length, which implies that at the free surface the breaking-induced flow is characterized by a dipole structure (Peregrine Reference Peregrine1999; Pizzo & Melville Reference Pizzo and Melville2013). The interaction of this flow with the wave field is thought to be significant for establishing Langmuir circulations (Leibovich Reference Leibovich1983), a crucial process for mixing the upper ocean. However, these classical theories do not take into account finite bandwidth effects, nor do they account for the two-way coupling between the wave and current fields.

Recently, there have been efforts to better describe two-way coupling effects (e.g. Phillips Reference Phillips2002; McWilliams et al. Reference McWilliams, Restrepo and Lane2004; Suzuki Reference Suzuki2019). Although these theories provide crucial insight, they are complicated and often obscure simple underlying physical constraints. Here, we have provided a simplified framework to examine wave–current interaction by assuming that the wave packets are compact, and that the currents are a collection of point vortices. Since this simplified system is derived from a variational principle, conservation laws arise naturally from the symmetries of the Lagrangian.

A central assumption made in this study is that the Doppler-shifted dispersion relationship serves as a faithful starting point to model wave–current interaction. That is, no additional terms are needed in the dispersion relationship to account for the vortical nature of the currents (see, for example, Stewart & Joy (Reference Stewart and Joy1974) for how vertical shear may modulate the dispersion relationship). Additionally, we assume that the currents (in the form of point vortices) and the wave packets are compact and widely spaced. A further simplifying but unnecessary assumption is that the wave packets do not interact with each other.

We have examined several solutions for the case of one wave packet and one vortex. These include stable bound orbits and unstable configurations. The wave packet and vortex may capture one another. We have also examined situations in which the wave packet and vortex collapse, occupying the same location at the same time. When collapse occurs, our theory breaks down. We also considered blow-up solutions, in which the modulus of the wavenumber grows without bound. In reality this growth would be arrested by wave breaking.

After examining a solution with two point vortices and one wave packet, we considered the motion induced by a wave packet on a weak point vortex in a horizontally periodic domain. It was shown that a net drift is induced by the wave packet, which may have implications for the advection of jetsam, flotsam and pollution at the ocean surface. This drift can also change the stability of a vortex street. Numerical calculations confirm that this system can be stable for long times, suggesting that this phenomenon might occur in nature.

This work motivates several future studies. In particular, the addition of the generation of waves by wind, wave packet–wave packet interaction and wave breaking, which creates a pair of oppositely signed vortices, would make this a more realistic description of the upper ocean.

Acknowledgements

We thank B. Young for helpful discussions. Additionally, we thank an anonymous referee for insightful comments on the far-field velocity of wave packets that led to a significantly improved manuscript.

Declaration of interests

The authors report no conflicts of interest.

Appendix A. Justification of Whitam's Lagrangian

This appendix offers a derivation of (2.1) following Whitham's averaged-Lagrangian method. Our starting point is the linear approximation to the Lagrangian of Miles (Reference Miles1977) – see also Luke (Reference Luke1967) and Zakharov (Reference Zakharov1968) – namely

(A1)\begin{equation} L[\phi,\eta]=\int \textrm{d}t \int \textrm{d}\kern0.07em x ( \phi \eta_t - H[\phi,\eta] ), \end{equation}

where

(A2)\begin{equation} H[\phi,\eta]= \tfrac{1}{2}g\eta^2 + \int_{-\infty}^0 \textrm{d}z \tfrac{1}{2}( \phi_x^2+\phi_z^2). \end{equation}

Here $\phi$ is the velocity potential and $\eta$ is the surface elevation. The integral in (A1) is over the sea surface $z=0$, and, for simplicity of notation, we ignore the $y$ direction. Variations $\delta \phi$, $\delta \eta$ yield the familiar linear equations and boundary conditions. A solution is

(A3a,b)\begin{equation} \eta(x,t)=A\cos (kx-\omega t), \quad \phi(x,z,t)=\frac{A\omega}{k}\,\textrm{e}^{kz}\sin (kx-\omega t), \end{equation}

where $A$ and $k$ are constants and $\omega ^2=gk$. Following Whitham, we substitute

(A4a,b)\begin{equation} \eta(x,t)=A(x,t)\cos (\theta(x,t)), \quad \phi(x,z,t)=\frac{A(x,t)\theta_t}{\theta_x}\,\textrm{e}^{\theta_x z}\sin (\theta(x,t)) \end{equation}

back into (A1) and (A2), obtaining

(A5)\begin{align} L[A,\theta]&=\iint \textrm{d}t\,\textrm{d}\kern0.07em x \frac{\omega^2 A^2 }{k}\sin^2\theta -\int \textrm{d} t \,H[A,\theta] \nonumber\\ &=\frac{1}{2}\iint \textrm{d}t\,\textrm{d}\kern0.07em x \frac{\omega^2 A^2 }{k} -\int \textrm{d}t \,H[A,\theta] \end{align}

and

(A6)\begin{align} H[A,\theta]&=\int {\textrm{d}\kern0.07em x} \, A^2 \left( \frac{1}{2}g\cos^2\theta +\frac{1}{4}\frac{\omega^2}{k}\cos^2\theta +\frac{1}{4}\frac{\omega^2}{k}\sin^2\theta \right) \nonumber\\ &= \frac{1}{4}\int {\textrm{d}\kern0.07em x} \, A^2 \left( g + \frac{\omega^2}{k} \right), \end{align}

where now $\omega =-\theta _t$ and $k=\theta _x$. In the final step of (A5) and (A6), we average over the fast dependence of $\theta (x,t)$. Combining results, we have

(A7)\begin{equation} L[A,\theta]= \frac{1}{4}\iint \textrm{d}t\,\textrm{d}\kern0.07em x \left( \frac{\omega^2}{k}-g \right)A^2 =\frac{1}{2}\iint \textrm{d}t\,\textrm{d}\kern0.07em x \left( \omega-\frac{gk}{\omega} \right) \mathcal{A}, \end{equation}

where, by (A6),

(A8)\begin{equation} E=\frac{1}{2} gA^2=\frac{1}{2} \frac{\omega^2}{k}A^2 \end{equation}

is the wave energy per unit area and $\mathcal {A}=E/\omega$ is the action. Independent variations of $A$ and $\theta$ are equivalent to independent variations of $\mathcal {A}$ and $\theta$, and either choice of variation yields the dispersion relation and the action conservation equation for surface waves. The Lagrangian

(A9)\begin{equation} L[\mathcal{A},\theta]=\iint \textrm{d}t\,\textrm{d}\kern0.07em x ( \omega-\sqrt{gk}) \mathcal{A} \end{equation}

yields the same two equations and is equivalent to (2.1).

Appendix B. Stream function generated by a single wave packet

Let the wave packet be located at $\boldsymbol {x}=0$. For $r\equiv \vert \boldsymbol {x} \vert \gg \vert \boldsymbol {x}' \vert$,

(B1)\begin{equation} \ln \vert \boldsymbol{x} - \boldsymbol{x}' \vert \approx \ln r - \frac {\boldsymbol{x} \boldsymbol{\cdot} \boldsymbol{x}'}{r^2} \end{equation}

and (2.35) becomes

(B2)\begin{equation} \psi_w(\boldsymbol{x})\approx \frac{\ln r}{2 {\rm \pi}} \iint \textrm{d}\kern0.07em \boldsymbol{x}' \, \rho(\boldsymbol{x}) - \frac{1}{2 {\rm \pi}} \frac {\boldsymbol{x}}{r^2} \boldsymbol{\cdot} \iint \textrm{d}\kern0.07em \boldsymbol{x}' \, \boldsymbol{x}' \rho(\boldsymbol{x}). \end{equation}

The first term in (B2) vanishes, because $\mathcal {A}=0$ at the boundary of the wave packet. In the second term,

(B3)\begin{equation} \iint \textrm{d}\kern0.07em \boldsymbol{x} \, \boldsymbol{x} \rho(\boldsymbol{x}) =\iint \textrm{d}\kern0.07em \boldsymbol{x} \, \boldsymbol{x} ( \mathcal{A}_x l_p -\mathcal{A}_y k_p) = \mathcal{A}_p ({-}l_p,k_p) \end{equation}

after integrations by parts, where

(B4)\begin{equation} \mathcal{A}_p = \iint \textrm{d}\kern0.07em \boldsymbol{x} \mathcal{A}. \end{equation}

Appendix C. Elimination of $\psi$ from the Lagrangian

In this appendix we justify the step of using the equations obtained by varying a particular field $\psi (\boldsymbol {x},t)$ to eliminate that same field from a Lagrangian that depends on several fields. We show that variations of the modified Lagrangian yield the correct equations for the remaining fields.

First we consider the related problem of finding the stationary points – maxima, minima or inflection points – of an ordinary function of two variables. Let the function be $f(\phi ,\psi )$. The stationary points are found by solving the set

(C1)\begin{equation} \frac{\partial}{\partial \phi} f(\phi,\psi)=f_1(\phi,\psi)=0 \end{equation}

and

(C2)\begin{equation} \frac{\partial}{\partial \psi} f(\phi,\psi)=f_2(\phi,\psi)=0, \end{equation}

where $f_1$ denotes the derivative of $f$ with respect to its first argument and $f_2$ denotes the derivative of $f$ with respect to its second argument. Suppose that (C2) can be solved explicitly for $\psi$ in the form

(C3)\begin{equation} \psi=g(\phi). \end{equation}

Then, substituting (C3) into (C1) we obtain a single equation for $\phi$, namely

(C4)\begin{equation} f_1(\phi,g(\phi))=0. \end{equation}

Our contention is that (C4) is equivalent to

(C5)\begin{equation} \frac{\partial}{\partial \phi} f(\phi,g(\phi))=0. \end{equation}

Clearly (C5) is equivalent to

(C6)\begin{equation} f_1(\phi,g(\phi))+f_2(\phi,g(\phi))g'(\phi)=0, \end{equation}

and the fact that (C3) solves (C2) means that $f_2(\phi ,g(\phi ))=0$. Thus (C5) is indeed equivalent to (C1).

To see that this proves our contention about the Lagrangian, we replace the integral over space and time by a sum over gridded values, and replace the derivatives of the field variables by finite differences. Then the Lagrangian becomes an ordinary function of many variables, namely the gridded values of the fields. We again regard this function as $f(\phi ,\psi )$ where now $\psi$ stands for a vector whose components are all the gridded values of $\psi$, and $\phi$ stands for a vector whose components are all the gridded values of $\phi$. If there are $N$ space–time grid points, then (C1) and (C2) each represent $N$ equations, but the essence of the proof is the same as that given above. It is easy to invent examples that show that the use of (C3) to eliminate some but not all of the $\psi$ terms in $f(\phi ,\psi )$ leads to erroneous results. The results given here border on the trivial, but the strategy of completely eliminating a field using the equations that result from the variations of that same field is important, because it seems to be one of the few legitimate methods of using the results of a variational principle to simplify the variational principle itself.

References

REFERENCES

Aref, H. 1983 Integrable, chaotic, and turbulent vortex motion in two-dimensional flows. Annu. Rev. Fluid Mech. 15, 345389.CrossRefGoogle Scholar
Bühler, O. 2014 Waves and Mean Flows. Cambridge University Press.CrossRefGoogle Scholar
Bühler, O. & Jacobson, T.E. 2001 Wave-driven currents and vortex dynamics on barred beaches. J. Fluid Mech. 449, 313339.CrossRefGoogle Scholar
Bühler, O. & McIntyre, M.E. 2003 Remote recoil: a new wave–mean interaction effect. J. Fluid Mech. 492, 207230.CrossRefGoogle Scholar
Bühler, O. & McIntyre, M.E. 2005 Wave capture and wave-vortex duality. J. Fluid Mech. 534, 6795.CrossRefGoogle Scholar
Darwin, C. 1953 Note on hydrodynamics. Math. Proc. Camb. Phil. Soc. 49 (2), 342354.CrossRefGoogle Scholar
Domm, U. 1956 Über die Wirbelstrassen von geringster Instabilität. Z. Angew. Math. Mech. 36 (9–10), 367371.CrossRefGoogle Scholar
Eames, I. & McIntyre, M.E. 1999 On the connection between Stokes drift and Darwin drift. In Mathematical Proceedings of the Cambridge Philosophical Society, vol. 126, pp. 171–174. Cambridge Philosophical Society.CrossRefGoogle Scholar
von Kármán, T. 1911 Über den Mechanismus des Widerstandes, den ein bewegter Körper in einer Flüssigkeit erfährt. Nachr. Ges. Wiss. Göttingen 1911, 509517.Google Scholar
Kirchhoff, G. 1883 Vorlesungen Über Mathematische Physik. Druck und Verlag von BG Teubner.Google Scholar
Lamb, H. 1932 Hydrodynamics. Cambridge University Press.Google Scholar
Leibovich, S. 1983 The form and dynamics of Langmuir circulations. Annu. Rev. Fluid Mech. 15 (1), 391427.CrossRefGoogle Scholar
Longuet-Higgins, M.S. & Stewart, R.W. 1962 Radiation stress and mass transport in gravity waves, with application to surf beats. J. Fluid Mech. 13 (4), 481504.CrossRefGoogle Scholar
Love, A.E.H. 1893 On the stability of certain vortex motions. Proc. Lond. Math. Soc. 1 (1), 1843.CrossRefGoogle Scholar
Luke, J.C. 1967 A variational principle for a fluid with a free surface. J. Fluid Mech. 27 (2), 395397.CrossRefGoogle Scholar
Maxwell, J.C. 1870 On the displacement in a case of fluid motion. Proc. Lond. Math. Soc. 3, 8287.Google Scholar
McIntyre, M.E. 2019 Wave–vortex interactions, remote recoil, the Aharonov–Bohm effect and the Craik–Leibovich equation. J. Fluid Mech. 881, 182217.CrossRefGoogle Scholar
McWilliams, J.C. 2016 Submesoscale currents in the ocean. In Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472 (2189), 20160117.CrossRefGoogle ScholarPubMed
McWilliams, J.C., Restrepo, J.M. & Lane, E.M. 2004 An asymptotic theory for the interaction of waves and currents in coastal waters. J. Fluid Mech. 511, 135178.CrossRefGoogle Scholar
Miles, J.W. 1977 On Hamilton's principle for surface waves. J. Fluid Mech. 83 (1), 153158.CrossRefGoogle Scholar
Morikawa, G.K. & Swenson, E.V. 1971 Interacting motion of rectilinear geostrophic vortices. Phys. Fluids 14 (6), 10581073.CrossRefGoogle Scholar
Morton, W.B. 1913 On the displacements of the particles and their paths in some cases of two-dimensional motion of a frictionless liquid. Proc. R. Soc. Lond. A 89 (608), 106124.Google Scholar
Onsager, L. 1949 Statistical hydrodynamics. Il Nuovo Cimento (1943–1954) 6 (2), 279287.CrossRefGoogle Scholar
Peregrine, D.H. 1999 Large-scale vorticity generation by breakers in shallow and deep water. Eur. J. Mech. B/Fluids 18 (3), 403408.CrossRefGoogle Scholar
Phillips, O.M. 1966 The Dynamics of the Upper Ocean. Cambridge University Press.Google Scholar
Phillips, W.R.C. 2002 Langmuir circulations beneath growing or decaying surface waves. J. Fluid Mech. 469, 317342.CrossRefGoogle Scholar
Pizzo, N.E. & Melville, W.K. 2013 Vortex generation by deep-water breaking waves. J. Fluid Mech. 734, 198218.CrossRefGoogle Scholar
Romero, L., Lenain, L. & Melville, W.K. 2017 Observations of surface wave–current interaction. J. Phys. Oceanogr. 47 (3), 615632.CrossRefGoogle Scholar
Salmon, R. 2016 Variational treatment of inertia–gravity waves interacting with a quasi-geostrophic mean flow. J. Fluid Mech. 809, 502529.CrossRefGoogle Scholar
Salmon, R. 2020 More lectures on geophysical fluid dynamics. Available at: http://pordlabs.ucsd.edu/rsalmon/More.Lecture.pdf.Google Scholar
Stewart, R.H. & Joy, J.W. 1974 HF radio measurements of surface currents. Deep-Sea Res. Oceanogr. Abstracts 21 (12), 10391049.CrossRefGoogle Scholar
Suzuki, N. 2019 On the physical mechanisms of the two-way coupling between a surface wave field and a circulation consisting of a roll and streak. J. Fluid Mech. 881, 906950.CrossRefGoogle Scholar
Tchieu, A.A., Kanso, E. & Newton, P.K. 2012 The finite-dipole dynamical system. Proc. R. Soc. A 468 (2146), 30063026.CrossRefGoogle Scholar
Tophøj, L. & Aref, H. 2013 Instability of vortex pair leapfrogging. Phys. Fluids 25 (1), 014107.CrossRefGoogle Scholar
Tur, A. & Yanovsky, V. 2017 Coherent Vortex Structures in Fluids and Plasmas. Springer.CrossRefGoogle Scholar
Wagner, G.L. & Young, W.R. 2015 Available potential vorticity and wave-averaged quasi-geostrophic flow. J. Fluid Mech. 785, 401424.CrossRefGoogle Scholar
Whitham, G.B. 1965 A general approach to linear and non-linear dispersive waves using a Lagrangian. J. Fluid Mech. 22 (2), 273283.CrossRefGoogle Scholar
Whitham, G.B. 1974 Linear and Nonlinear Waves. John Wiley & Sons.Google Scholar
Zakharov, V.E. 1968 Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys. 9 (2), 190194.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

Figure 2. (a) Orbit of wave packet and vortex with $\mathcal {R}_0 = \mathcal {R}_0^+$ as defined in (3.30ad). 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.

Figure 2

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.

Figure 3

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.

Figure 4

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.

Figure 5

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.

Figure 6

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$. (ad) 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$.

Figure 7

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.