Particle description of the interaction between wave packets and point vortices

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 \emph{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 analyze 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.


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 x p (t), its wave action A p , and its wave vector k p (t). Each point vortex i is defined by its location x i (t) and its strength Γ 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 A p and Γ i are conserved.
The velocity field attached to the wave packets is dipolar; it is sometimes called '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 x p (t), k p (t), and 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 2020), 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 analyze 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 symmetrical arrangement permits thorough analysis, which is confirmed by numerical solutions. This solution is very similar to the one discussed by Bühler & McIntyre (2005) 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 1870;Darwin 1953) 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 center 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.
§8 concludes with an assessment of our results and their oceanographic implications.

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 (1965Whitham ( , 1974 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 (2020).
For the waves by themselves the Lagrangian proposed by Whitham is where the integral is over time and the ocean surface; the frequency ω = −θ t and wave vector k = ∇θ are abbreviations for the derivatives of the wave phase θ(x, t); is the wave action; E is the wave energy per unit area; and ω r (k) is the prescribed relative frequency of the waves-the frequency measured in a reference frame moving at the mean flow velocity U(x, t). Our notation is k = (k, l), x = (x, y), and ∇ = (∂ x , ∂ y ). For surface waves, At this stage, we consider the mean flow to be prescribed. For sake of completeness, Appendix A provides a systematic derivation of (2.1) following Whitham's averaged-Lagrangian method. Variations of A yield the dispersion relation, where c g (k) = ∂ω r /∂k is the relative group velocity. Thus we obtain the action conservation equation, The Lagrangian for two-dimensional incompressible flow is where the subscript m stands for 'mean flow', and is the Jacobian, defined for any two functions A(x, y) and B(x, y). The variables α(x, t), β(x, t) and ψ(x, t) represent averages over the constant depth H 0 to which wave-mean interactions occur. We identify H 0 with the decay depth of the surface waves. We assume that the mean flow is depth-independent in this range. Stationarity of L m implies δα : β t + [ψ, β] = 0, (2.9) δβ : α t + [ψ, α] = 0, (2.10) We see that α and β are vorticity labels in the following sense: First, by (2.9) and (2.10), they are conserved following the fluid motion, hence each fluid particle is identified by its two labels (α, β). Second, by (2.11), the vorticity in an arbitrary area of the flow is given by 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, we obtain the vorticity equation, for the mean flow by itself. We couple L w to L m by replacing the mean velocity U in (2.1) with u ψ ≡ (−ψ y , ψ x ), and by assuming that the Lagrangian for the entire system is the sum, is the Hamiltonian, and (θ, A) and (α, β) are canonical pairs. Using (2.2), (2.11) and integrations by parts to evaluate (2.16) we find that (2.17) Thus, as expected, our dynamics conserves the sum of the wave energy and the kinetic energy of the mean flow. The equations corresponding to δL = 0 are where ∇ × (A, B) ≡ 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 is the pseudomomentum. The wave action equation (2.5) is unchanged, but now U = u ψ . That is, the previously arbitrary mean flow is now specifically identified with the velocity field (−ψ y , ψ 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 ∇ 2 ψ − ∇× 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, ∇ 2 ψ = p = 0, and hence (2.26) 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 p before the broad mean flow represented by ψ 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 ∇ω = −∇θ t = −k t , we obtain the refraction equation The refractive change in k predicted by (2.27) causes a change in ω r (k) that can be determined from (2.3). If the waves do not break, then the action A = E/ω r (k) is conserved. If ω r (k) increases, then the wave energy E must also increase to keep their ratio constant. We anticipate that wave-vector stretching, which increases |k| and hence ω 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) where x i (t) is the location at time t of a point vortex with strength Γ i . The subscripts i replace the continuous vorticity labels α and β. The Hamiltonian (2.16) becomes (2.29) To fully convert from (α, β) to x i we must transform the term dxdt αβ t (2.30) in (2.15). It becomes dxdydt α ∂(x, y, β) ∂(x, y, t) = dαdβdτ α ∂(x, y, β) ∂(α, β, τ ) (2.31) Thus, when point vortices replace continuous vorticity, the Lagrangian (2.15) becomes (2.32) Instead of (2.20) we now have where (suppressing the time-dependence) We now assume that the wave field consists solely of isolated wave packets. The stream function field generated by a single wave packet at x p is given by (2.35) and (2.36) with k = k p , where k p (t) is the wave vector associated with the wave packet. We assume that 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 A = 0. In Appendix B we show that, far from x p , the streamfunction generated by a wavepacket at x p takes the form of a dipole, where A p = dxA is the total action of the wave packet. The streamfunction response to many point vortices and many wave packets is clearly is the response to a monopole at x i , and is the response to a dipole with wavevector k p at x p . The constants Γ i and A p measure the strength of the monopole and the dipole, respectively. A p is always positive but Γ i can have either sign. Until dissipation occurs A p and Γ i remain constant. Since our aim is to produce a Lagrangian that depends only on the point vortex locations x i (t), the wave packet locations x p (t), and their wave vectors 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 where we have used integrations by parts and the relation which follows from the definition of x p (t): dx p (t)/dt is the velocity of the wave envelope. The second term in (2.32) becomes The three terms in (2.32) containing ψ combine as where we have used (2.33).
Our final step is to substitute (2.38) back into the Lagrangian, removing its dependence on ψ. The last integral in (2.44) becomes (2.45) 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 Putting all this together, we obtain the Lagrangian 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, Γ i and A p are constants. The Hamiltonian (2.48) contains Γ Γ terms and Γ A terms. If we had not dropped the dipole/dipole interactions it would also contain AA 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 legal to use (2.33) to eliminate ψ from (2.32).
The equations corresponding to (2.47)-(2.48) are is the velocity field induced by the point vortices, and is the velocity field induced by the wave packets. The total velocity is U( . 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 U d (x p ) to (2.49). Equations (2.49-2.51) are the fundamental equations of our model. If we were to regard 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 2014). Similarly, if we omit U d , then (2.51) is the standard equation of point vortex dynamics (Kirchhoff 1883). The new feature of our derivation is that 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 U d would also appear in (2.49) and (2.50). Tchieu et al. (2012) 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) corresponds to space-translation symmetry and is proved by considering variations of the form where ǫ(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 k p ) caused by vortex i is accompanied by a change in the position of vortex i. Bühler & McIntyre (2003) 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 A p is suddenly replaced by two counter-rotating vortices with a dipole moment equal to Γ D where D is the separation between counterrotating vortices of strength ±Γ . See also Bühler & McIntyre (2005); Bühler & Jacobson (2001), and Bühler (2014). Our dynamics also conserves the angular momentum, which can be proved by considering variations of the form ( Our primary motivation in deriving (2.49-2.51) is to obtain relatively simple equations, amenable to analytic and numerical solution, governing the interactions between ocean surface waves and the vorticity created by breaking waves. One could avoid the assumption of widely spaced wave packets and vortex patches by solving the more general equations (2.18-2.27) as coupled equations for the fields of A(x, t), k(x, t), and ψ(x, t). However, even these more general equations are idealized in the sense that they embody our treatment of the vertical dimension. 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. The three-dimensional structure of this vortex tube is locally important, but only the nearly vertical portions of the vortex tube induce a significant surface flow far from the tube. The arbitrary constant depth H 0 to which the tube's contribution extends is an artificial component of our model and could easily be absorbed into other parameters. We prefer to retain it as a constant reminder of the very idealized nature of our model. Models, to be useful, must be much simpler than reality. Onsager (1949) considered the equilibrium statistical mechanics of a system of point vortices. Our system reduces to Onsager's system when no waves are present (A p ≡ 0). Our phase space is larger than the one considered by Onsager because it contains dimensions corresponding to the wave packet locations x p and their wave vectors 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 −∞ < k p < ∞. We therefore expect an ultraviolet catastrophe in which energy spreads to ever larger |k p | 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 (2005), Wagner & Young (2015), and Salmon (2016). For many further details, see Salmon (2020). In the remainder of this paper we investigate the dynamics (2.49)-(2.51).

1 vortex, 1 wave packet
We begin by considering the system consisting of a single vortex of strength Γ located at x(t), and a single wavepacket with action A and wave vector k(t) located at x(t)+ξ(t). This system exhibits a much more complicated range of behavior than the more familiar system consisting of two point vortices. The Langrangian (2.47) takes the form with Hamiltonian The equations of motion become where ξ = (ξ, η). We simplify notation by taking g = 1 and choosing a characteristic wavenumber k 0 = 1 so that ω 2 0 = gk 0 = 1. We also assume H 0 = k −1 0 = 1, while we take (3.8) The system (3.3)-(3.7) conserves the energy (3.8), the angular momentum We use the conserved momenta (3.10) to eliminate the variables (x, y) in favor of (k, l, ξ, η). The resulting system conserves the energy (3.8) and the quantity obtained by eliminating (x, y) between (3.9) and (3.10). We also define The reduced dynamics takes the form of four coupled ordinary differential equations, ( 3.16) with the two conserved quantities, (3.11) and (3.12). Define We shall obtain a single, closed equation for the wavenumber magnitude κ(t). First, using (3.8) and (3.12), we obtain an expression for a 2 in terms of κ, (3.18) Then using (3.8) we obtain the constraint on the phases. From (3.15) and (3.16), we find an evolution equation for φ, namelẏ Equations (3.15) and (3.16) also implẏ Combining (3.20) and (3.21), we obtaiṅ Our final step is to eliminateφ to arrive at equation involving onlyκ and κ. We use the identity Substituting (3.24) back into (3.22) we obtain the closed evolution equation (3.26) Equation (3.25) takes the form of a particle moving in a potential Π(κ). This permits a qualitative analysis of system behavior based on the form of (3.26). Solutions may be written out in implicit form, as in Tur & Yanovsky (2017), but a qualitative analysis offers better physical insight.

Circular motion
We begin by seeking solutions that exhibit simple harmonic motion. Thus we takė κ = 0 and look for the κ i that satisfy Π(κ i ) = 0.
(3.27) Let κ i ≡ 1. This implies a simple relation between H 0 and R 0 . Its solutions are H 0 = 1 and R 0 a free parameter; or (3.28) If we take κ i = 1 to be a critical point, then Π ′ (1) = 0. It may be shown that when H 0 = 1, Π ′′ (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 The solutions take the form k = cos φ, l = sin φ, ξ = a sin φ, η = a cos φ, where a is given by (3.18). From (3.20) we have φ = φ 0 + φ 1 t for φ 0 and φ 1 constants. An example of this behavior 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.

Collapse
The above analysis addresses local spectral stability at a critical point. There are other solutions in which a → 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 1983, and references therein). The case of overlapping vorticity and wave action has been analyzed by McIntyre (2019).
The conditions for collapse are clear from equation (3.18). Collapse occurs at the time t * at which (3.32) As an example we suppose that H 0 = R 0 = 0. Then the system collapses as κ → 0. Under these assumptions, the governing equation for κ reduces tȯ  This can be integrated, and we arrive at t = t 0 ± 1 4 (12ϑ − 8 sin 2ϑ + sin 4ϑ) . Thus the collapse corresponds to the formation of a cusp in κ. Figure 3 confirms these results. The top right panel shows the convergence of the wave packet and the vortex. The bottom right panel compares the theoretical prediction of κ(t) (where we have taken the negative branch of the solution corresponding toκ < 0) to the numerical result. The two curves are indistinguishable.
In our model the wave action A is fixed. Therefore, the wave energy A ω(κ) vanishes as κ → 0 since ω(κ) ∝ √ κ 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.

Blow-up
There are also solutions in which κ, the wave number 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 (2001). In this paper we consider only ideal solutions, and we do not include wave breaking.
A necessary condition for κ(t) to grow without bound is that In the bottom panel, we display κ. 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.
for large κ. For large κ, hence blow-up requires (3.42) As an example, we take H 0 = −1. Then which implies that the blow-up solution takes the form We examine this numerically in figure 4, and find agreement with the theoretical prediction.

Two vortices and one wave packet
We now examine the system comprising a single wave packet with action A p and wavevector (k p , l p ) located at (x p , y p ); a point vortex of strength −Γ located at x 1 , y 1 ; and a second point vortex of strength +Γ located at x 2 , y 2 . Refer to figure 5. Initially, and, by symmetry, these conditions hold at all later times. The Lagrangian is (with We vary all the dependent variables, and then apply the symmetry conditions (4.1) to obtain a closed set of four equations. (It is illegal to apply the symmetry condition before taking the variations.) The reduced set of equations is where c g is the x-component of the group velocity, and is the squared distance between the wave packet and either vortex. Because of the symmetry conditions (4.1), we do not need the evolution equations for y p , l p , x 2 , and y 2 . The equations (4.3)-(4.6) conserve energy in the form and momentum in the form The angular momentum vanishes. Defining we rewrite (4.3)-(4.6) as three equationṡ in the three unknowns k p , y 1 and X, where now d 2 = X 2 + y 2 1 . 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) 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 = −∞ with k p > 0, approaching the vortex pair from the left, as shown in the top panel of figure 5. 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 ω r 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. k p reaches its maximum value at X = 0, where (4.15) 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 → ∞.
The numerical solution shown in figure 5 confirms this analysis.

Wave-packet induced drift
If Γ ≪ A, (4.9-4.13) implyẊ = c g (k p ) − A p k p 2πd 4 (X 2 − y 2 1 ), (4.16) k p = 0, (4.17) Hence k p and c g are constants. The governing equations (with A p = 2πa 2 and k p = g = 1) reduce toẊ (4.20) 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 (1870, see also Morton 1913Darwin 1953). In the reference frame moving with the wave packet, the stream function is an integral of motion. Hence is constant. DefineẊ =Ẋ − c g0 and θ = tan −1 (y 1 /X). Using (4.21) and following Maxwell (1870) and Darwin (1953) we find that X = a 2 cos 2θ Let cos θ = −sn(ν) with the suppressed modulus of the Jacobi elliptic function understood to be (4.23) It is tedious but straight-forward to show that Similar expressions may be found for y 1 (ν) and t(ν) but fall outside the scope of our discussion (Darwin 1953). From (4.24), the total drift in the x-direction is where E, K are the complete elliptical integrals of the first and second kind, respectively. The drift volume, D, defined as (4.26) has the same order of magnitude as the vertically integrated Stokes drift. Note, the connection between Stokes drift (specifically the motion in the vertical plane) and Darwin drift has been examined by Eames & McIntyre (1999).

1 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 on a ring of geostrophic vortices in Morikawa & Swenson 1971). The single vortex of strength Γ = 2π remains stationary at x = 0. The N wave packets at x p = (x p , y p ) have equal actions A p = A, and wave vectors of equal magnitude |k p | = κ. They lie symmetrically on the circle |x p | = χ. Both κ and χ are constants. We take the arbitrary depth parameter to be unity, H 0 = 1.
A vortex at x induces the velocity field We also need The wavepacket at x p induces the velocity at the vortex. The equations of motion reduce tȯ and p A(−2l p x p y p − k p (x 2 p − y 2 p ), −2k p x p y p + l p (x 2 p − y 2 p )) = 0. (5.8) The last equation is the condition that the vortex remains stationary. We look for solutions exhibiting simple harmonic motion with k p = κ(sin θ p , − cos θ p ) and x p = χ(cos θ p , sin θ p ). This leads to the constraints and A(cos θ p , sin θ p ) = 0. (5.11) Taking A = 1, we observe that solutions to this system are related to the N -th roots of unity. This sets the initial phase of θ p . We find that (5.12) Figure 6 shows a solution of this type with four wave packets.

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 (1870) and Darwin (1953). 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πn, y), where n is any integer. Again we take H 0 = 1.
The governing equations are (2.49) -(2.51) with U m and U d now calculated from the stream functions and (6.2)

Limit of weak point vortices
We begin by considering the wave-packet-induced motion of the vortices when |Γ i | ≪ 1. As in §4-see particularly §4.1-this motion takes a nontrivial form.
We consider a wave packet traveling along the x-axis so that l = y p = 0. To O(1), (2.49)-(2.51) reduce toẋ where U d is computed from (6.2). As we are considering one (weak) vortex per period, i = 1 and we take y 1 = y.
Although we have been unable to find closed form solutions, asymptotic analysis reveals interesting properties of this system. We expand the transcendental pieces of (6.13) and (6.14) as sinh y sin χ (cosh y − cos χ) 2 = 2 ∞ n=1 ne −ny sin nχ (6.15) and 1 − cosh y cos χ (cosh y − cos χ) 2 = −2 ∞ n=1 ne −ny cos nχ. (6.16) This expansion assumes y > 0 ; a similar formula holds for y < 0. The equations of motion may then be written aṡ ne −ny cos nχ, (6.17) ne −ny sin nχ, (6.18) or, more compactly, asż Let b > 0 to be the initial vertical coordinate of the particle, and (for reasons that become clear in §7) let π/2 be its horizontal coordinate. Let the initial location of the wave packet be (π, 0). As the wave packets and vortex are sufficiently far apart, a natural small parameter arises, namely (6.20) We expand z as z = z 0 + ǫz 1 + ǫ 2 z 2 + · · · (6.21) To O(ǫ 2 ), we finḋ We solve this system iteratively. To lowest order, where z 0 0 = −π/2 + ib. This implies where χ 0 = −π/2 − c g0 t and we define θ ≡ c g0 t. (6.25) 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) so that upon substitution of the result for z 1 , we find (6.27) At this order there is a mean drift in the direction of wave propagation. Unlike the Stokes drift which is O(A), this term is O(A 2 ). This has potentially important implications for the stability of vortex streets, as discussed in §7. Figure 7 compares our approximate analytical solutions to numerical integrations of the full equations for ǫ = 1/10 (left panel) and ǫ = 1/3 (right panel). We see that the asymptotic theory works relatively well for small values of ǫ.

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 = ±b, and spaced π apart in the x-direction. Refer to figure 8. This is the minimal system to illustrate the instability of a periodic vortex street (Domm 1956). The stability of vortex streets was first considered by von Kármán (1911), and discussed in detail by Lamb (1932, §156). Domm (1956) examined its 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 1956). The equations of motion are (2.49) -(2.51) with the stream functions given by (6.1) and (6.2). We take Γ i = ǫ 2 γ so that Γ 1 = Γ 2 = ǫ 2 γ and Γ 3 = Γ 4 = −ǫ 2 γ. We seek equations of motion valid to O(ǫ 2 ). We begin by expanding the variables describing the wave packet as x p = x p0 + c g0 t + ǫ 2 x p2 , y p = 0, (7.1) 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 iṡ Expanding hyperbolic functions, we find thaṫ where γ ≡ γ 1 . Eqn (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 1932, §156). Expandinġ 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 (1932, §156) finds that the self advection of the vortex street leads to a uniform translation of the vortices with speed U m and U d to O(ǫ 2 δ). At this order, the vortex induced flow is the same as it would be in the absence of the wave packet, hence, from Lamb (1932), we have the O(ǫ 2 δ) result (7.14) At this order U m (z 1 ) and U m (z 2 ) only depend on vortex 1 and 2, hence we omit U m (z 3 ) and 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(δ) contributions to U d (z 1,2 ) are found by substituting the expansion given by (7.8) into (6.19) so that To solve this system of equations, we must also expand the perturbations in ǫ, so that We assume the time dependence of the order zero terms goes like e Λt , where Λ = ǫ 2 λ (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 1932). 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 z δ11 = i µ 2c g0 z * δ10 (e iθ − 1), (7.17) and z δ21 = −i µ 2c g0 z * δ20 (e iθ − 1). (7.18) We now have sufficient information to solve for the stability of the vortex street configuration, which will be dictated by the evolution of {z δ10 , z δ20 }. Substituting the lower order expansions into U d and U m , we find that the dynamics of {z δ10 , z δ20 } are given by the phase averaged equations of motion. The evolution equations are given bẏ where σ = µ 2 2c g0 ; ν = γ 4π . (7.21) Define z δ10 = αe ǫ 2 λt , z δ20 = βe ǫ 2 λt . Then, (7.19-7.20), together with their complex conjugates, imply the following eigenvalue problem The eigenvalues λ are found to be λ 2 = ν(ν + σ). (7.23) In the limit that σ = 0, we recover the result of Lamb (1932, §156) that the symmetric vortex street is linearly unstable. It follows from (7.23) that the system is stable when We note that ν can be positive or negative depending on the sign of γ, so that the inequality may hold only if the sign of ν and σ 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 ν < 0, the system remains linearly stable.

Long-time behavior
In the previous subsection, we analyzed the linear stability of the vortex street to perturbations. However, as was shown by Domm (1956), 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×2π). We take ǫ = 1/4, k 0 = 1, A p = 1, Γ = ǫ 2 , and we integrate the system for 3 × 10 4 seconds. This is shown in the left panel of figure 9, where it is seen that the vortex street remains stable for long times in this configuration. In the absence of the wave packet (A p = 0), shown on the right panel of figure 9, the system is unstable, and the vortex street eventually dissolves, analogous to the behaviour found for two vortex pairs (Love 1893;Tophøj & Aref 2013).

Conclusion
It has long been recognized in the oceanographic community that surface waves may freely exchange momentum and energy with underlying currents. The dynamics are governed by wave action conservation (Longuet-Higgins & Stewart 1962;Whitham 1965), 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 1966), there are still many fundamental questions regarding wave -current interaction, including a need to better understand their two way coupling (McWilliams et al. 2004). This is thrown in to particularly stark relief by numerical models of climate, which are beginning to resolve the submesoscale (on the order of 1-10km), where these interactions may be especially pronounced (McWilliams 2016;Romero et al. 2017). 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 1999;Pizzo & Melville 2013). The interaction of this flow with the wave field is thought to be significant for establishing Langmuir circulations (Leibovich 1983), 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.  figure 8. The wave packet motion is shown in blue, while the vortices are shown in black (Γ i < 0) and red (Γ i > 0). The left panel shows that the vortex street is stable over this long time integration. The right panel shows 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. 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.
The 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 (1974) to see 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 possible implications for the advection of jetsam, flotsam and pollution at the ocean surface. This drift also has implications for the stability of a vortex street in the presence of a wave packet. Our analysis shows that the wave packet may stabilize the vortex street.
Numerical calculations confirm that this system is stable for long times, suggesting that this phenomenon might be observable 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.
Declaration of Interests: The authors report no conflicts of interest.
is the wave energy per unit area, and A = E/ω is the action. Independent variations of A and θ are equivalent to independent variations of A and θ, and either choice of variation yields the dispersion relation and the action conservation equation for surface waves. The Lagrangian yields the same two equations and is equivalent to (2.1).
Our contention is that ( 17) is equivalent to ∂ ∂φ f (φ, g(φ)) = 0 ( 18) Clearly ( 18) is equivalent to f 1 (φ, g(φ)) + f 2 (φ, g(φ))g ′ (φ) = 0 and the fact that ( 16) solves ( 15) means that f 2 (φ, g(φ)) = 0. Thus ( 18) is indeed equivalent to ( 14). To see that this proves our contention about the Lagrangian, 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 (φ, ψ) where now ψ stands for a vector whose components are all the gridded values of ψ, and φ stands for a vector whose components are all the gridded values of φ. If there are N spacetime gridpoints, then ( 14) and ( 15) 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 ( 16) to eliminate some but not all of the ψ-terms in f (φ, ψ) 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.