1 Introduction
Flows consisting of a discrete dispersed phase carried by a continuous fluid phase are tremendously common both in nature and in industrial applications (Balachandar & Eaton Reference Balachandar and Eaton2010). Such multiphase flows are well known to show a range of complex behaviours. When the dispersed phase is present at high volume fraction, for example, it can modify the rheology of the fluid (Stickel & Powell Reference Stickel and Powell2005). But even a low-volume-fraction dispersed phase can be difficult to describe. Recent decades have seen the development of a large body of work to characterize the motion of transported particles by homogeneous fluids (Guha Reference Guha2008). The case of spherical particles that behave inertially due to a difference in density from the carrier fluid has been well studied, and although details remain to be worked out, their essential phenomenology is by now fairly well characterized (Toschi & Bodenschatz Reference Toschi and Bodenschatz2009; Balachandar & Eaton Reference Balachandar and Eaton2010). There is less consensus, however, about particles whose effective inertia comes from their finite size relative to the smallest length scales in the flow (Qureshi et al. Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007; Ouellette, O’Malley & Gollub Reference Ouellette, O’Malley and Gollub2008; Brown, Warhaft & Voth Reference Brown, Warhaft and Voth2009; Homann & Bec Reference Homann and Bec2010) or particles that are aspherical so that their shape plays a role in their dynamics (Shin & Koch Reference Shin and Koch2005; Andersson & Soldati Reference Andersson and Soldati2013; Voth & Soldati Reference Voth and Soldati2017).
Because such anisotropic particles are relevant in applications ranging from paper making (Lundell, Söderberg & Alfredsson Reference Lundell, Söderberg and Alfredsson2011) to the dynamics of icy clouds (Pinsky & Khain Reference Pinsky and Khain1998) to the fate of microplastic pollution in the ocean (Ryan et al. Reference Ryan, Moore, van Franeker and Moloney2009; Chubarenko et al. Reference Chubarenko, Bagaev, Zobkov and Esiukova2016; Sherman & Van Sebille Reference Sherman and Van Sebille2016), there has been a significant amount of recent work on their behaviour in flows (Voth & Soldati Reference Voth and Soldati2017). Because anisotropic particles couple to the flow in complex ways, however, the detailed structure of the flow may influence their behaviour more than it would for spherical particles. Interest in paper making or similar fibre-laden industrial flows has led to a number of studies of the dynamics of anisotropic particles in wall-bounded shear flows (Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008a ,Reference Mortensen, Andersson, Gillissen and Boersma b ; Challabotla, Zhao & Andersson Reference Challabotla, Zhao and Andersson2015a ,Reference Challabotla, Zhao and Andersson b ). Cloud physics and related motivations have spurred work to understand how anisotropic particles move in isotropic turbulence (Parsa et al. Reference Parsa, Calzavarini, Toschi and Voth2012; Ni, Ouellette & Voth Reference Ni, Ouellette and Voth2014; Byron et al. Reference Byron, Einarsson, Gustavsson, Voth, Mehlig and Variano2015; Ni et al. Reference Ni, Kramel, Ouellette and Voth2015). Flows in the ocean with relevance to microplastic transport, however, are distinct from these two cases, since the near-shore environment has its own characteristic fluid mechanics. Of particular relevance to this application are surface gravity waves (Isobe et al. Reference Isobe, Kubo, Tamura, Kako, Nakashima and Fujii2014; Chubarenko & Stepanova Reference Chubarenko and Stepanova2017; Hinata et al. Reference Hinata, Mori, Ohno, Miyao and Kataoka2017), which produce flows that are qualitatively different from wall-bounded shear flows or isotropic turbulence.
The tumbling of anisotropic particles in deep-water waves has previously been studied in the context of developing models for the effective viscosity of grease ice (de Carolis, Olla & Pignagnoli Reference de Carolis, Olla and Pignagnoli2005; Olla Reference Olla2006), where it was argued that the strain rate produced by the wave field controls the tumbling rate of the particles. More recently, we conducted a numerical study of the behaviour of both non-inertial and weakly inertial anisotropic particles under linear surface gravity waves (DiBenedetto, Ouellette & Koseff Reference DiBenedetto, Ouellette and Koseff2018). We found that the shape of these spheroidal particles strongly couples to their transport mechanics, primarily because they adopt dynamically preferred orientations (with a superimposed oscillation) with respect to the wave field. These preferred orientations modify their lift and drag relative to a sphere of comparable volume and mass, leading to differences in horizontal and vertical transport. Here, we follow on this work to illuminate the physical mechanism responsible for this preferential orientation. Since the equations of motion for inertial spheroids are very complicated, we consider the somewhat simplified case of neutrally buoyant, non-inertial point particles moving in the flow field generated by linear surface gravity waves, which we can treat analytically. We show that the preferred orientation is a consequence of finite wave amplitude and arises in a way that is analogous to the appearance of Stokes drift in the translational dynamics of a particle. We show analytically that the preferred orientation of the particle is simply related to its aspect ratio, just as we observed previously (DiBenedetto et al. Reference DiBenedetto, Ouellette and Koseff2018). Since we find this preferred orientation even without particle inertia, we argue that it is not an inertial effect. Finally, we also determine the magnitude of the oscillation of the particle about this preferred orientation and its phase lag relative to the orbital motion.
We begin in § 2 by determining the equations of motion that govern the dynamics of the particle. Then, in § 3, we solve these equations in a series of approximations, following the derivation of classic Stokes drift. We show that to lowest order the particle orientation simply drifts with time, but that the first-order nonlinear correction introduces a fixed point into the wave-averaged orientation of the particles. Finally, in § 4, we summarize our results and discuss the extensions necessary to connect our results to the case of real microplastic particles in the environment.
2 Equations of motion
We consider here a flow driven by linear surface gravity waves. We define the unit vector
$\hat{\boldsymbol{e}}_{x}$
to point in the direction of wave propagation and the unit vector
$\hat{\boldsymbol{e}}_{z}$
to point antiparallel to gravity. We assume no flow in the
$y$
direction. The flow field is then given by
$$\begin{eqnarray}\displaystyle u_{x} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\cosh [k(z+H)]}{\sinh kH}\cos (kx-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle u_{z} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\sinh [k(z+H)]}{\sinh kH}\sin (kx-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$k$
is the wavenumber of the wave,
$\unicode[STIX]{x1D714}$
is the frequency,
$A$
is the wave amplitude and
$H$
is the total depth of the water column. The dispersion relation
$\unicode[STIX]{x1D714}^{2}=gk\tanh kH$
where
$g$
is acceleration due to gravity, constrains the relationship between
$k$
,
$\unicode[STIX]{x1D714}$
and
$H$
. These expressions are derived from Airy wave theory, which assumes irrotational flow and linearizes the free surface boundary conditions by assuming small-amplitude waves. The results are valid as long as the waves are not too steep, requiring
$kA\ll 1$
. The only non-vanishing components of the rate of strain are
$S_{xx}=-S_{zz}$
and
$S_{xz}=S_{zx}$
, which are given by 
$$\begin{eqnarray}\displaystyle S_{xx} & = & \displaystyle -kA\unicode[STIX]{x1D714}\frac{\cosh [k(z+H)]}{\sinh kH}\sin (kx-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{xz} & = & \displaystyle kA\unicode[STIX]{x1D714}\frac{\sinh [k(z+H)]}{\sinh kH}\cos (kx-\unicode[STIX]{x1D714}t).\end{eqnarray}$$
The particles are assumed to be infinitesimal point particles, although they retain a non-trivial shape. We assume spheroidal particles. The aspect ratio
$\unicode[STIX]{x1D706}$
of the particles is given by the ratio
$a/b$
where
$a$
is the length of the particle along its symmetry axis and
$b$
is the length of the particle along one of its orthogonal axes. It is often more convenient to work in terms of the eccentricity
$\unicode[STIX]{x1D700}$
, given by
$$\begin{eqnarray}\unicode[STIX]{x1D700}=\frac{\unicode[STIX]{x1D706}^{2}-1}{\unicode[STIX]{x1D706}^{2}+1}.\end{eqnarray}$$
One can also write
$\unicode[STIX]{x1D706}$
in terms of
$\unicode[STIX]{x1D700}$
as
$$\begin{eqnarray}\unicode[STIX]{x1D706}=\sqrt{\frac{1+\unicode[STIX]{x1D700}}{1-\unicode[STIX]{x1D700}}}.\end{eqnarray}$$
Spheroids with
$\unicode[STIX]{x1D700}>0$
are prolate, while those with
$\unicode[STIX]{x1D700}<0$
are oblate; spheres have
$\unicode[STIX]{x1D700}=0$
. Note that
$\unicode[STIX]{x1D700}$
is bounded and lies between
$-1$
and
$1$
by construction.
We assume that the particles have no inertia relative to the flow. Thus, their translational degrees of freedom will be identical to those of a fluid element, and their rotational degrees of freedom will evolve according to Jeffery’s (Reference Jeffery1922) equation. Defining
$\boldsymbol{p}$
as a unit vector that points along the particle’s axis of symmetry, Jeffery’s equation gives
$$\begin{eqnarray}{\dot{p}}_{i}=\unicode[STIX]{x1D6FA}_{ij}p_{j}+\unicode[STIX]{x1D700}(\unicode[STIX]{x1D61A}_{ij}p_{j}-p_{i}p_{j}\unicode[STIX]{x1D61A}_{jk}p_{k}),\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}_{ij}$
is the rate of rotation tensor (which is zero in this flow) and
$\unicode[STIX]{x1D61A}_{ij}$
is the rate of strain tensor.
In general, the particle can point in any direction, and so we should consider all three components of Jeffery’s equation. Here, these can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{x}=S_{xx}p_{x}(1-(p_{x}^{2}-p_{z}^{2}))+S_{xz}p_{z}(1-2p_{x}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{y}=-p_{y}[S_{xx}(p_{x}^{2}-p_{z}^{2})+2S_{xz}p_{x}p_{z}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{z}=S_{xz}p_{x}(1-2p_{z}^{2})-S_{xx}p_{z}(1+(p_{x}^{2}-p_{z}^{2})). & \displaystyle\end{eqnarray}$$
By inspection, the equation of motion for
$p_{y}$
(2.6b
) has two fixed points where its right-hand side vanishes: one when
$p_{y}=1$
(and therefore
$p_{x}=p_{z}=0$
, since
$\boldsymbol{p}$
is a unit vector), and one when
$p_{y}=0$
. The first of these fixed points, which corresponds to the particle being aligned orthogonal to the plane of the waves, is a global fixed point for the full dynamical system, since the right-hand sides of (2.6a
) and (2.6c
) also vanish for this case. A straightforward linear stability analysis shows that this fixed point is unstable (see appendix A).
The second fixed point for
$p_{y}$
, corresponding to the particle lying fully in the plane of the waves, is somewhat more complex to analyse since it is not a global fixed point; when
$p_{y}=0$
, the only obvious constraint on
$p_{x}$
and
$p_{z}$
is that
$p_{x}^{2}+p_{z}^{2}=1$
, and so the right-hand sides of (2.6a
) and (2.6c
) do not vanish. However, our previous numerical work showed that particles will (over some period of time) align into the plane of the waves regardless of their initial orientation (DiBenedetto et al.
Reference DiBenedetto, Ouellette and Koseff2018), suggesting that the set corresponding to
$p_{y}=0$
is stable. Thus, for now we will make the ansatz that we can neglect the
$p_{y}$
dynamics and assume that the particles will align into the wave plane. As we will show, this assumption allows us to derive a limit cycle in the
$p_{x}$
and
$p_{z}$
dynamics, which (as we show in appendix A) is consistent with
$p_{y}=0$
defining a stable attracting set.
Once the particles have aligned into the wave plane, this becomes a two-dimensional problem, and we can work more conveniently in polar coordinates. Measuring the polar angle
$\unicode[STIX]{x1D719}$
down from the
$z$
axis (so that
$\unicode[STIX]{x1D719}=0$
when
$\boldsymbol{p}=\hat{\boldsymbol{e}}_{z}$
and
$\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}/2$
when
$\boldsymbol{p}=\hat{\boldsymbol{e}}_{x}$
), we can write
$$\begin{eqnarray}\displaystyle p_{x} & = & \displaystyle \sin \unicode[STIX]{x1D719},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{z} & = & \displaystyle \cos \unicode[STIX]{x1D719},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{p}}_{x} & = & \displaystyle \dot{\unicode[STIX]{x1D719}}\cos \unicode[STIX]{x1D719},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{p}}_{z} & = & \displaystyle -\dot{\unicode[STIX]{x1D719}}\sin \unicode[STIX]{x1D719}.\end{eqnarray}$$
$p_{y}=0$
and transforming, both the
$p_{x}$
and
$p_{z}$
equations of motion ((2.6a
) and (2.6c
)) consistently reduce to 
$$\begin{eqnarray}\unicode[STIX]{x1D700}^{-1}\dot{\unicode[STIX]{x1D719}}=S_{xx}\sin 2\unicode[STIX]{x1D719}+S_{xz}\cos 2\unicode[STIX]{x1D719}.\end{eqnarray}$$
If the strain rate were independent of time, there would be a fixed point in the
$\unicode[STIX]{x1D719}$
dynamics at
$$\begin{eqnarray}\unicode[STIX]{x1D719}^{\ast }=-\frac{1}{2}\tan ^{-1}\left(\frac{S_{xz}}{S_{xx}}\right),\end{eqnarray}$$
where
$^{\ast }$
denotes a fixed point. This point would moreover be independent of particle shape and only related to the strain angle in the flow, as has been observed in previous work (Gay Reference Gay1968; Reed & Tryggvason Reference Reed and Tryggvason1974). In our problem, however, the strain rate is a function of time, and so the dynamical system is non-autonomous; thus, there is no simple fixed point for
$\unicode[STIX]{x1D719}$
in this flow.
Figure 1. The evolution of the polar angle
$\unicode[STIX]{x1D719}$
as found by numerically integrating the dynamical system in (2.11). In all cases,
$kA=0.2$
,
$H=10$
and
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$
, the initial orientation of the particle was set at
$\unicode[STIX]{x1D719}(t=0)=0$
, and the initial position was set to
$x_{c}(t=0)=0$
and
$z_{c}(t=0)$
such that the particle remained just below the free surface even at its maximum vertical excursion. Each curve corresponds to a different particle eccentricity; from top to bottom,
$\unicode[STIX]{x1D700}=0.5$
, 0.1,
$-0.1$
and
$-0.5$
. In each case,
$\unicode[STIX]{x1D719}$
asymptotically tends to a limit cycle of sinusoidal oscillations about a preferred value.
Combining (2.1a
), (2.1b
) and (2.9), we can write the full dynamical system (assuming the particle has already rotated into the plane of the waves) governing the evolution of the coordinates of the particle centre
$(x_{c},z_{c})$
and its orientation
$\unicode[STIX]{x1D719}$
as
$$\begin{eqnarray}\displaystyle {\dot{x}}_{c} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\cosh [k(z_{c}+H)]}{\sinh kH}\cos (kx_{c}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}}_{c} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\sinh [k(z_{c}+H)]}{\sinh kH}\sin (kx_{c}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D719}} & = & \displaystyle \frac{kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH} [\!-\cosh [k(z_{c}+H)]\sin (kx_{c}-\unicode[STIX]{x1D714}t)\sin 2\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle +\,\sinh [k(z_{c}+H)]\cos (kx_{c}-\unicode[STIX]{x1D714}t)\cos 2\unicode[STIX]{x1D719}\!].\end{eqnarray}$$
$\unicode[STIX]{x1D645}$
of this dynamical system is given by 
$$\begin{eqnarray}\displaystyle \text{Tr}\,\unicode[STIX]{x1D645} & = & \displaystyle \frac{\unicode[STIX]{x2202}{\dot{x}}_{c}}{\unicode[STIX]{x2202}x_{c}}+\frac{\unicode[STIX]{x2202}{\dot{z}}_{c}}{\unicode[STIX]{x2202}z_{c}}+\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH} [\!\cosh [k(z_{c}+H)]\sin (kx_{c}-\unicode[STIX]{x1D714}t)\cos 2\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle +\,\sinh [k(z_{c}+H)]\cos (kx_{c}-\unicode[STIX]{x1D714}t)\sin 2\unicode[STIX]{x1D719}\!].\end{eqnarray}$$
As long as
$\unicode[STIX]{x1D700}\neq 0$
(that is, a non-spherical particle),
$\text{Tr}\,\unicode[STIX]{x1D645}\neq 0$
, and so via Liouville’s theorem there is the possibility of attractors in the dynamics. Thus, while there may be no fixed point for the particle orientation, there can be an attracting limit cycle. And indeed, when we numerically integrate equation (2.11), this is exactly the behaviour we find, as shown in figure 1. We note that since the Stokes drift speed, with which the particle moves in a wave-period-averaged sense, is not the same as the wave propagation speed, there will be a small Doppler shift of the wave frequency as seen in the reference frame of the particle. This slight change in the apparent value of
$\unicode[STIX]{x1D714}$
does not, however, impact the results below, and so we do not explicitly account for it.
3 Analysis
The equations of motion given in (2.11) are very complex, in that they are coupled, nonlinear and non-autonomous. While we can integrate them numerically, making analytical progress requires some approximations.
3.1 Linear approximation
To analyse the flows and transport produced by surface gravity waves, a common first approximation is to linearize the equations of motion for the particle centre, given that the wave steepness is assumed to be small (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2016). As is well known, in this limit fluid particles will move on closed elliptical orbits. In that spirit, let us write the position of the spheroid
$(x_{c},z_{c})$
as the sum of the (constant) position of the centre of the orbit
$(x_{o},z_{o})$
and deviations from this centre (that is, coordinates on the orbit)
$(\unicode[STIX]{x1D709},\unicode[STIX]{x1D701})$
, so that
$$\begin{eqnarray}\displaystyle x_{c} & = & \displaystyle x_{o}+\unicode[STIX]{x1D709},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle z_{c} & = & \displaystyle z_{o}+\unicode[STIX]{x1D701}.\end{eqnarray}$$
For small-amplitude waves (
$kA\ll 1$
), we can linearize the equations of motion and approximate
$$\begin{eqnarray}\displaystyle {\dot{x}}_{c} & = & \displaystyle u_{x}(x_{c},z_{c},t)\approx u_{x}(x_{o},z_{o},t)=\dot{\unicode[STIX]{x1D709}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}}_{c} & = & \displaystyle u_{z}(x_{c},z_{c},t)\approx u_{z}(x_{o},z_{o},t)=\dot{\unicode[STIX]{x1D701}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D719}} & = & \displaystyle f(x_{c},z_{c},t)\approx f(x_{o},z_{o},t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D709}} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\cosh [k(z_{o}+H)]}{\sinh kH}\cos (kx_{o}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D701}} & = & \displaystyle A\unicode[STIX]{x1D714}\frac{\sinh [k(z_{o}+H)]}{\sinh kH}\sin (kx_{o}-\unicode[STIX]{x1D714}t).\end{eqnarray}$$
$x_{o}$
and
$z_{o}$
are constant values that refer to the centre of the wave orbital (Kundu et al.
Reference Kundu, Cohen and Dowling2016). Thus, if we average these equation over a wave period, we find that
$\dot{\bar{x}}_{c}=0$
and
$\dot{\bar{z}}_{c}=0$
, where we have defined the wave-period average by, for example, 
$$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D709}}=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D709}\,\text{d}(kx_{o}-\unicode[STIX]{x1D714}t). & & \displaystyle\end{eqnarray}$$
Thus, we recover the expected closed elliptical orbits for the particle motion, as these equations can be directly integrated to obtain
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}(t) & = & \displaystyle -A\frac{\cosh [k(z_{o}+H)]}{\sinh kH}\sin (kx_{o}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}(t) & = & \displaystyle A\frac{\sinh [k(z_{o}+H)]}{\sinh kH}\cos (kx_{o}-\unicode[STIX]{x1D714}t).\end{eqnarray}$$
The equation of motion for the polar angle at this order of approximation is given by
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D719}} & = & \displaystyle \frac{kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH} [\!-\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\sin 2\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle +\,\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\cos 2\unicode[STIX]{x1D719}\!].\end{eqnarray}$$
Unlike for the translational equations of motion,
$\dot{\bar{\unicode[STIX]{x1D719}}}\neq 0$
, meaning that the particle will reorient as it moves on the elliptical orbit. This equation cannot be simply integrated in closed form, as it remains nonlinear and non-autonomous, but numerically integrating the equation shows that
$\unicode[STIX]{x1D719}$
does not tend to a constant value with a superimposed oscillation as we found from integrating the full equations of motion (2.11). Instead, the particle tumbles, with
$\unicode[STIX]{x1D719}$
increasing without bound linearly with a superimposed oscillation as shown in figure 2. At this level of approximation, spheroids do not show a preferred orientation. Thus, we can hypothesize that just as Stokes drift of the orbital centre is a nonlinear, finite-amplitude effect, so too is the preferred orientation of the spheroid.
Figure 2. (a) Numerical solutions of the linearized equation of motion for
$\unicode[STIX]{x1D719}$
(3.6) for
$\unicode[STIX]{x1D700}=0.5$
, 0.1,
$-0.1$
and
$-0.5$
. The two larger values of
$|\unicode[STIX]{x1D700}|$
are the steeper curves. Positive values of
$\unicode[STIX]{x1D700}$
are shown with solid lines, and negative values with dashed lines; the difference is simply a
$180^{\circ }$
phase shift. In all cases,
$kA=0.2$
,
$H=10$
,
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$
,
$x_{o}=0$
and
$z_{o}$
was chosen such that the particle always remained just below the free surface. (b) The same data shown only for short times, to make the phase difference between the positive and negative values of
$\unicode[STIX]{x1D700}$
clearer.
3.2 Finite-amplitude correction
To check this hypothesis, we can account for the leading-order finite-wave-amplitude correction to the particle motion by including the next term in the Taylor expansion in the approximation we made in (3.2a ). That is, we now write
$$\begin{eqnarray}\displaystyle {\dot{x}}_{c} & = & \displaystyle u_{x}(x_{c},z_{c},t)\approx u_{x}(x_{o},z_{o},t)+\unicode[STIX]{x1D709}\left.\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}\right|_{(x_{o},z_{o})}+\unicode[STIX]{x1D701}\left.\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}\right|_{(x_{o},z_{o})},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{z}}_{c} & = & \displaystyle u_{z}(x_{c},z_{c},t)\approx u_{z}(x_{o},z_{o},t)+\unicode[STIX]{x1D709}\left.\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}\right|_{(x_{o},z_{o})}+\unicode[STIX]{x1D701}\left.\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}\right|_{(x_{o},z_{o})},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D719}} & = & \displaystyle f(x_{c},z_{c},t)\approx f(x_{o},z_{o},t)+\unicode[STIX]{x1D709}\left.\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}x}\right|_{(x_{o},z_{o})}+\unicode[STIX]{x1D701}\left.\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}z}\right|_{(x_{o},z_{o})}.\end{eqnarray}$$
At this level of approximation, we expect to recover Stokes drift for the translational motion of the particle. Inserting the forms of
$\unicode[STIX]{x1D709}$
and
$\unicode[STIX]{x1D701}$
given in (3.5) and simplifying, the translational equations of motion become
$$\begin{eqnarray}\displaystyle {\dot{x}}_{c} & {\approx} & \displaystyle A\unicode[STIX]{x1D714}\frac{\cosh [k(z_{o}+H)]}{\sinh kH}\cos (kx_{o}-\unicode[STIX]{x1D714}t)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{kA^{2}\unicode[STIX]{x1D714}}{\sinh ^{2}(kH)}[\cosh ^{2}[k(z_{o}+H)]\sin ^{2}(kx_{o}-\unicode[STIX]{x1D714}t)+\sinh ^{2}[k(z_{o}+H)]\cos ^{2}(kx_{o}-\unicode[STIX]{x1D714}t)]\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle {\dot{z}}_{c}\approx A\unicode[STIX]{x1D714}\frac{\sinh [k(z_{o}+H)]}{\sinh kH}\sin (kx_{o}-\unicode[STIX]{x1D714}t). & & \displaystyle\end{eqnarray}$$
When averaged over the wave period, these equations give the correct Stokes drift forms of
$$\begin{eqnarray}\displaystyle \dot{\bar{x}}_{c}=kA^{2}\unicode[STIX]{x1D714}\frac{\cosh [2k(z_{o}+H)]}{2\sinh ^{2}kH} & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle \dot{\bar{z}}_{c}=0, & & \displaystyle\end{eqnarray}$$
so that the spheroid will slowly drift in the direction of propagation of the waves but will not (on average) rise or sink in the water column.
Figure 3. The evolution of the polar angle
$\unicode[STIX]{x1D719}$
as found by numerically integrating the approximation given by (3.12). In all cases,
$kA=0.2$
,
$H=10$
,
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$
,
$x_{o}=0$
and
$z_{o}$
was chosen such that the particle always remained just below the free surface. Again, each curve corresponds to a different particle eccentricity; from top to bottom,
$\unicode[STIX]{x1D700}=0.5$
, 0.1,
$-0.1$
and
$-0.5$
. With the finite wave-amplitude correction in (3.12), we recover the oscillation of
$\unicode[STIX]{x1D719}$
about a preferred orientation.
Similar calculations give the new equation of motion for the polar angle
$\unicode[STIX]{x1D719}$
as
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D719}} & {\approx} & \displaystyle \frac{kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH} [\!-\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\sin 2\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle +\,\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\cos 2\unicode[STIX]{x1D719}\!]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{(kA)^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{2\sinh ^{2}kH}[\sin [2(kx_{o}-\unicode[STIX]{x1D714}t)]\sin 2\unicode[STIX]{x1D719}+\sinh [2k(z_{o}+H)]\cos 2\unicode[STIX]{x1D719}].\end{eqnarray}$$
At this level of approximation, an important new feature emerges in the equation of motion: the last term in the equation is now autonomous. This new piece modifies the dynamics significantly, changing the linear drift of
$\unicode[STIX]{x1D719}$
seen at the lowest order of approximation to a steady oscillation about a preferred value, as shown in figure 3. Thus, we find an intriguing symmetry in the spheroid dynamics. If we neglect the finite amplitude of the waves, the particle centre oscillates about a fixed location but its orientation drifts linearly. But when we account for the finite amplitude, the particle centre drifts linearly but its orientation oscillates about a fixed value.
Physically, this result suggests that the same mechanism that is responsible for Stokes drift is also responsible for fixing the orientation of the particle at a preferred value. However, concluding much more from (3.12) is difficult, since it is complicated, nonlinear and non-autonomous. Thus, making more analytical progress requires further approximations.
3.3 Change of variables
A large part of the difficulty in analysing (3.12) further is due to the appearance of
$\unicode[STIX]{x1D719}$
on its right-hand side only in transcendental functions. Thus, let us introduce a change of variables to remove this explicit transcendental dependence.
Let
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}=\sin 2\unicode[STIX]{x1D719} & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}=\cos 2\unicode[STIX]{x1D719}. & & \displaystyle\end{eqnarray}$$
The time derivatives of these new variables are related to
$\dot{\unicode[STIX]{x1D719}}$
by
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D712}} & = & \displaystyle 2\dot{\unicode[STIX]{x1D719}}\cos 2\unicode[STIX]{x1D719}=2\unicode[STIX]{x1D713}\dot{\unicode[STIX]{x1D719}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D713}} & = & \displaystyle -2\dot{\unicode[STIX]{x1D719}}\sin 2\unicode[STIX]{x1D719}=-2\unicode[STIX]{x1D712}\dot{\unicode[STIX]{x1D719}}.\end{eqnarray}$$
$|\unicode[STIX]{x1D712}|\leqslant 1$
,
$|\unicode[STIX]{x1D713}|\leqslant 1$
and
$\unicode[STIX]{x1D712}^{2}+\unicode[STIX]{x1D713}^{2}=1$
. With this change of variables, the linearized equation of motion for
$\unicode[STIX]{x1D719}$
(3.6) can be transformed into the pair of coupled equations
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D712}} & = & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}[-\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}+\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D713}^{2}],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D713}} & = & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}[\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}^{2}-\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D712}} & = & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}[-\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}+\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D713}^{2}]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{(kA)^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh ^{2}kH}[\sin [2(kx_{o}-\unicode[STIX]{x1D714}t)]\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}+\sinh [2k(z_{o}+H)]\unicode[STIX]{x1D713}^{2}],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D713}} & = & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}[\cosh [k(z_{o}+H)]\sin (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}^{2}-\sinh [k(z_{o}+H)]\cos (kx_{o}-\unicode[STIX]{x1D714}t)\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{(kA)^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh ^{2}kH}[\sin [2(kx_{o}-\unicode[STIX]{x1D714}t)]\unicode[STIX]{x1D712}^{2}+\sinh [2k(z_{o}+H)]\unicode[STIX]{x1D712}\unicode[STIX]{x1D713}].\end{eqnarray}$$
$(\unicode[STIX]{x1D712},\unicode[STIX]{x1D713})$
phase-plane dynamics of this 2-D system is relatively simple. The solutions lie on the unit circle in the
$(\unicode[STIX]{x1D712},\unicode[STIX]{x1D713})$
phase plane. In the linearized case, the solutions move at a constant rate around this circle with a superimposed oscillation at the wave frequency, while in the nonlinear case they oscillate on an arc on the unit circle around their stationary value (after an initial transient), as shown in figure 4.
Figure 4. The particle orientational dynamics in the
$(\unicode[STIX]{x1D712},\unicode[STIX]{x1D713})$
phase plane. At lowest order, the dynamics of
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
move along the unit circle at a constant rate with a superimposed oscillation (thin solid line). When accounting for finite wave amplitudes, they instead oscillate on arcs on this circle (thick lines). Data are shown here for the same cases plotted in figure 3. In all cases,
$kA=0.2$
,
$H=10$
,
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$
,
$x_{o}=0$
and
$z_{o}$
was chosen such that the particle always remained just below the free surface, and data are shown for eccentricities
$\unicode[STIX]{x1D700}=0.5$
(lower right arc), 0.1 (upper right arc),
$-0.1$
(lower left arc) and
$-0.5$
(upper left arc).
Empirically, the solutions to these coupled equations consist of a fast mode that oscillates at the wave frequency
$\unicode[STIX]{x1D714}$
and a slow mode that evolves on a much longer time scale. The fast oscillation for
$\unicode[STIX]{x1D712}$
is very close to
$180^{\circ }$
out of phase with the wave driving and for
$\unicode[STIX]{x1D713}$
is
$90^{\circ }$
out of phase. These observations suggest the decomposition
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}(t)=\bar{\unicode[STIX]{x1D712}}-B\sin (kx_{o}-\unicode[STIX]{x1D714}t) & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(t)=\bar{\unicode[STIX]{x1D713}}+B\cos (kx_{o}-\unicode[STIX]{x1D714}t), & & \displaystyle\end{eqnarray}$$
where
$\bar{\unicode[STIX]{x1D712}}$
and
$\bar{\unicode[STIX]{x1D713}}$
are still functions of time but with characteristic time dynamics that are not at the wave frequency. One can think of
$\bar{\unicode[STIX]{x1D712}}$
and
$\bar{\unicode[STIX]{x1D713}}$
as
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
with the Fourier mode at the wave frequency (with the appropriate phase) removed. We note that this decomposition neglects a small phase shift between
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
and the wave driving; as we show, however, our results agree well with the full dynamics, providing an a posteriori justification for this approximation.
We make the requirement that the amplitude
$B$
of oscillation at the wave frequency is the same for both
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
because the two variables are not independent. Since by definition
$\unicode[STIX]{x1D712}=\sin 2\unicode[STIX]{x1D719}$
and
$\unicode[STIX]{x1D713}=\cos 2\unicode[STIX]{x1D719}$
,
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
must have the same overall amplitude; thus, their response to the wave driving also has the same amplitude. At this stage, however,
$B$
remains an unknown non-dimensional function of the physical parameters of the system (that is,
$k$
,
$A$
,
$\unicode[STIX]{x1D714}$
,
$H$
,
$x_{o}$
,
$z_{o}$
and
$\unicode[STIX]{x1D700}$
); we will set its value later.
We can now substitute these forms for
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
into their equations of motion and take a wave-period average. We begin with the linearized dynamics (3.16a
) and (3.16b
). The nonlinear products of
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
that appear in these equations are given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}\unicode[STIX]{x1D713} & = & \displaystyle \bar{\unicode[STIX]{x1D712}}\bar{\unicode[STIX]{x1D713}}-B^{2}\sin (kx_{o}-\unicode[STIX]{x1D714}t)\cos (kx_{o}-\unicode[STIX]{x1D714}t)\nonumber\\ \displaystyle & & \displaystyle +\,B\bar{\unicode[STIX]{x1D712}}\cos (kx_{o}-\unicode[STIX]{x1D714}t)-B\bar{\unicode[STIX]{x1D713}}\sin (kx_{o}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}^{2} & = & \displaystyle \bar{\unicode[STIX]{x1D713}}^{2}+2B\bar{\unicode[STIX]{x1D713}}\cos (kx_{o}-\unicode[STIX]{x1D714}t)+B^{2}\cos ^{2}(kx_{o}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}^{2} & = & \displaystyle \bar{\unicode[STIX]{x1D712}}^{2}-2B\bar{\unicode[STIX]{x1D712}}\sin (kx_{o}-\unicode[STIX]{x1D714}t)+B^{2}\sin ^{2}(kx_{o}-\unicode[STIX]{x1D714}t).\end{eqnarray}$$
$\bar{\unicode[STIX]{x1D712}}$
and
$\bar{\unicode[STIX]{x1D713}}$
are still functions of time, and so they cannot be passed through a wave-period average. But since their time variation is slow compared to the wave-induced oscillation, we can approximate them as constant over each wave period. If we do that, many terms in the equations vanish due to orthogonality relations between sines and cosines, and we are left with 
$$\begin{eqnarray}\displaystyle \dot{\bar{\unicode[STIX]{x1D712}}} & {\approx} & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\left(\frac{1}{2}\cosh [k(z_{o}+H)]+\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D713}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\bar{\unicode[STIX]{x1D713}}} & {\approx} & \displaystyle -\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\left(\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D712}}.\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \ddot{\bar{\unicode[STIX]{x1D712}}} & = & \displaystyle -\left(\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\right)^{2}\left(\frac{1}{2}\cosh [k(z_{o}+H)]+\sinh [k(z_{o}+H)]\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\left(\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D712}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \ddot{\bar{\unicode[STIX]{x1D713}}} & = & \displaystyle -\left(\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\right)^{2}\left(\frac{1}{2}\cosh [k(z_{o}+H)]+\sinh [k(z_{o}+H)]\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\left(\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D713}}.\end{eqnarray}$$
$\bar{\unicode[STIX]{x1D712}}$
and
$\bar{\unicode[STIX]{x1D713}}$
are simple harmonic oscillators at the ‘slow’ frequency: 
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{s} & = & \displaystyle \left(\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\right)\left[\left(\frac{1}{2}\cosh [k(z_{o}+H)]+\sinh [k(z_{o}+H)]\right)\right.\nonumber\\ \displaystyle & & \displaystyle \times \left.\left(\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]\right)\right]^{1/2}.\end{eqnarray}$$
Note that this frequency depends on
$B$
, since the original equations were nonlinear. This result also consistently recovers the purely linear increase of
$\bar{\unicode[STIX]{x1D719}}$
with time that we expect at this level of approximation.
Let us now include the finite-amplitude corrections in our wave-period average, using (3.17a ) and (3.17b ). These give us
$$\begin{eqnarray}\displaystyle \dot{\bar{\unicode[STIX]{x1D712}}} & = & \displaystyle \frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\left(\frac{1}{2}\cosh [k(z_{o}+H)]+\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{(kA)^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh ^{2}kH}\left[-\frac{B^{2}}{4}+\sinh [2k(z_{o}+H)]\left(\bar{\unicode[STIX]{x1D713}}^{2}+\frac{B^{2}}{2}\right)\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{\bar{\unicode[STIX]{x1D713}}} & = & \displaystyle -\frac{2kA\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh kH}B\left(\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]\right)\bar{\unicode[STIX]{x1D712}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{(kA)^{2}\unicode[STIX]{x1D714}\unicode[STIX]{x1D700}}{\sinh ^{2}kH}\sinh [2k(z_{o}+H)]\bar{\unicode[STIX]{x1D712}}\bar{\unicode[STIX]{x1D713}}.\end{eqnarray}$$
$\bar{\unicode[STIX]{x1D712}}$
and
$\bar{\unicode[STIX]{x1D713}}$
are now quite different, which is part of what leads to the new behaviour we see when including the finite-amplitude correction. Importantly, they are no longer harmonic oscillators. But what is most interesting is that these autonomous differential equations now contain fixed points, which we will label
$\bar{\unicode[STIX]{x1D712}}^{\ast }\neq 0$
and
$\bar{\unicode[STIX]{x1D713}}^{\ast }\neq 0$
. Let us consider the
$\dot{\bar{\unicode[STIX]{x1D713}}}$
(3.24b
). Assuming that
$\bar{\unicode[STIX]{x1D712}}^{\ast }\neq 0$
(which, empirically, it is not), we can set
$\dot{\bar{\unicode[STIX]{x1D713}}}=0$
, divide through by
$\bar{\unicode[STIX]{x1D712}}$
, and find the fixed point
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D713}}^{\ast }=-B\frac{2\sinh kH}{kA}\left(\frac{\cosh [k(z_{o}+H)]+{\textstyle \frac{1}{2}}\sinh [k(z_{o}+H)]}{\sinh [2k(z_{o}+H)]}\right).\end{eqnarray}$$
Given the constraint that
$\unicode[STIX]{x1D712}^{2}+\unicode[STIX]{x1D713}^{2}=1$
, the fixed point for
$\bar{\unicode[STIX]{x1D712}}$
is then
$$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D712}}^{\ast }=\sqrt{1-\bar{\unicode[STIX]{x1D713}}^{\ast 2}}, & & \displaystyle\end{eqnarray}$$
and thus the preferred orientation of the spheroid about which it oscillates is
$$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D719}}^{\ast }={\textstyle \frac{1}{2}}\cos ^{-1}\bar{\unicode[STIX]{x1D713}}^{\ast }. & & \displaystyle\end{eqnarray}$$
However, we are not yet done, because (3.25) still depends on the unknown oscillation amplitude
$B$
.
To set the value of
$B$
, we can apply physical constraints to (3.25). As mentioned above, in general
$B$
can be a (non-dimensional) function of all the physical parameters in the problem:
$k$
,
$A$
,
$\unicode[STIX]{x1D714}$
,
$H$
,
$x_{o}$
,
$z_{o}$
and
$\unicode[STIX]{x1D700}$
. It is explicitly not, however, a function of time, by construction. Now, since the left-hand side of (3.25) is independent of time and the right-hand side has no explicit dependence on
$\unicode[STIX]{x1D714}$
,
$B$
also cannot be a function of
$\unicode[STIX]{x1D714}$
since the time units cannot be cancelled. Additionally, we know that by construction
$|\bar{\unicode[STIX]{x1D713}}^{\ast }|\leqslant 1$
, since it is defined as the cosine of an angle. The factors aside from
$B$
on the right-hand side of (3.25), however, can take any values (since they are not constrained at all by the motion of the particle but rather are properties of the wave) and can certainly lie outside this range. Thus, the only way to guarantee that the right-hand side of (3.25) lies between
$-1$
and 1 for all wave parameters is to require
$B$
to cancel out the rest of the right-hand side; that is,
$$\begin{eqnarray}\displaystyle B=F(\unicode[STIX]{x1D700})\frac{kA}{2\sinh kH}\left(\frac{\sinh [2k(z_{o}+H)]}{\cosh [k(z_{o}+H)]+{\textstyle \frac{1}{2}}\sinh [k(z_{o}+H)]}\right), & & \displaystyle\end{eqnarray}$$
where
$F(\unicode[STIX]{x1D700})$
is some unknown function of
$\unicode[STIX]{x1D700}$
alone. But since
$\unicode[STIX]{x1D700}$
itself lies between
$-1$
and 1 and is non-dimensional, the simplest choice for
$F$
that fits our constraints is
$F(\unicode[STIX]{x1D700})=\unicode[STIX]{x1D700}$
. That choice gives
$$\begin{eqnarray}\displaystyle B=\frac{kA\unicode[STIX]{x1D700}}{2\sinh kH}\left(\frac{\sinh [2k(z_{o}+H)]}{\cosh [k(z_{o}+H)]+\frac{1}{2}\sinh [k(z_{o}+H)]}\right) & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D713}}^{\ast }=-\unicode[STIX]{x1D700}. & & \displaystyle\end{eqnarray}$$
We can now set
$\bar{\unicode[STIX]{x1D719}}^{\ast }$
. We have
$$\begin{eqnarray}\displaystyle \cos 2\bar{\unicode[STIX]{x1D719}}^{\ast }=2\cos ^{2}\bar{\unicode[STIX]{x1D719}}^{\ast }-1=-\unicode[STIX]{x1D700}, & & \displaystyle\end{eqnarray}$$
so
$$\begin{eqnarray}\displaystyle \cos \bar{\unicode[STIX]{x1D719}}^{\ast }=\sqrt{{\textstyle \frac{1}{2}}(1-\unicode[STIX]{x1D700})}. & & \displaystyle\end{eqnarray}$$
Standard trigonometric identities then give
$$\begin{eqnarray}\displaystyle \tan \bar{\unicode[STIX]{x1D719}}^{\ast }=\frac{\sqrt{1-{\textstyle \frac{1}{2}}(1-\unicode[STIX]{x1D700})}}{\sqrt{{\textstyle \frac{1}{2}}(1-\unicode[STIX]{x1D700})}}=\sqrt{\frac{1+\unicode[STIX]{x1D700}}{1-\unicode[STIX]{x1D700}}}=\unicode[STIX]{x1D706}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}$
is the aspect ratio of the particle, just as we reported previously (DiBenedetto et al.
Reference DiBenedetto, Ouellette and Koseff2018). Finally, let us note that this value of
$B$
allows us to write the slow frequency
$\unicode[STIX]{x1D714}_{s}$
that appears in the linearized dynamics as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{s}=\unicode[STIX]{x1D714}\left(\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\right)^{2}\sinh [2k(z_{o}+H)]\left[\frac{2\tanh [k(z_{o}+H)]+1}{\tanh [k(z_{o}+H)]+2}\right]^{1/2}. & & \displaystyle\end{eqnarray}$$
3.4 Oscillation amplitude
Now, although
$B$
is the amplitude of the oscillation of
$\unicode[STIX]{x1D712}$
and
$\unicode[STIX]{x1D713}$
about their fixed points, it is not the amplitude of the oscillation of
$\unicode[STIX]{x1D719}$
. However, knowing
$\bar{\unicode[STIX]{x1D719}}^{\ast }$
, we can compute the amplitude of this oscillation by going back to (3.12). Empirically,
$\unicode[STIX]{x1D719}$
oscillates about
$\bar{\unicode[STIX]{x1D719}}^{\ast }$
at the wave frequency; however, there is a phase lag between the wave and the oscillation of
$\unicode[STIX]{x1D719}$
, and in this case it is not negligible. Thus, at steady state, we can write
$$\begin{eqnarray}\unicode[STIX]{x1D719}(t)=C\sin (kx_{o}-\unicode[STIX]{x1D714}t)+D\cos (kx_{o}-\unicode[STIX]{x1D714}t)+\bar{\unicode[STIX]{x1D719}}^{\ast }.\end{eqnarray}$$
$C$
and
$D$
will be different due to the unknown phase of the
$\unicode[STIX]{x1D719}$
oscillation. The amplitude
$A_{\unicode[STIX]{x1D719}}$
of the oscillation can be found by adding them in quadrature, since sines and cosines are orthogonal:
$$\begin{eqnarray}\displaystyle A_{\unicode[STIX]{x1D719}}=\sqrt{C^{2}+D^{2}}. & & \displaystyle\end{eqnarray}$$
With this form, we have
$$\begin{eqnarray}\dot{\unicode[STIX]{x1D719}}=-C\unicode[STIX]{x1D714}\cos (kx_{o}-\unicode[STIX]{x1D714}t)+D\unicode[STIX]{x1D714}\sin (kx_{o}-\unicode[STIX]{x1D714}t).\end{eqnarray}$$
We can now insert this form into (3.12) and evaluate it at two different values of the wave phase
$kx_{o}-\unicode[STIX]{x1D714}t$
: since it is now an algebraic equation, it must hold at all phases. If we choose
$(kx_{o}-\unicode[STIX]{x1D714}t)=0$
, we obtain
$$\begin{eqnarray}C=-\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\sinh [k(z_{o}+H)]\cos [2(D+\bar{\unicode[STIX]{x1D719}}^{\ast })]\left[1+\frac{kA}{\sinh kH}\cosh [k(z_{o}+H)]\right].\end{eqnarray}$$
If alternatively we choose
$(kx_{o}-\unicode[STIX]{x1D714}t)=\unicode[STIX]{x03C0}/2$
, we obtain
$$\begin{eqnarray}\displaystyle D & = & \displaystyle -\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\cosh [k(z_{o}+H)]\sin [2(C+\bar{\unicode[STIX]{x1D719}}^{\ast })]\nonumber\\ \displaystyle & & \displaystyle \times \,\left[1-\frac{kA}{\sinh kH}\sinh [k(z_{o}+H)]\cot [2(C+\bar{\unicode[STIX]{x1D719}}^{\ast })]\right].\end{eqnarray}$$
The second term in brackets in each of these equations is of order
$(kA)^{2}$
, and arises from the finite-wave-amplitude correction to the equation of motion for
$\unicode[STIX]{x1D719}$
. But while these corrections are very important for modifying the dynamics of
$\unicode[STIX]{x1D719}$
(in that they are the origin of the preferential alignment), they do not play a significant role in setting the magnitude of the oscillation. Thus, we can drop these higher-order terms in setting the oscillation amplitude
$A_{\unicode[STIX]{x1D719}}$
, and write
$$\begin{eqnarray}\displaystyle C & {\approx} & \displaystyle -\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\sinh [k(z_{o}+H)]\cos [2(D+\bar{\unicode[STIX]{x1D719}}^{\ast })],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle D & {\approx} & \displaystyle -\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\cosh [k(z_{o}+H)]\sin [2(C+\bar{\unicode[STIX]{x1D719}}^{\ast })].\end{eqnarray}$$
$C$
and
$D$
, we can compute the oscillation amplitude
$A_{\unicode[STIX]{x1D719}}$
via (3.36). In figure 5, we show the predictions for
$\bar{\unicode[STIX]{x1D719}}^{\ast }$
and
$A_{\unicode[STIX]{x1D719}}$
along with the same data as shown in figure 1. The agreement is excellent.
Figure 5. The same curves of
$\unicode[STIX]{x1D719}$
as shown in figure 1, now augmented with our predictions for the preferred alignment
$\bar{\unicode[STIX]{x1D719}}^{\ast }$
((3.33); solid horizontal lines) and the oscillation amplitude
$A_{\unicode[STIX]{x1D719}}$
((3.36); dashed horizontal lines). As above, from top to bottom the curves correspond to the eccentricities
$\unicode[STIX]{x1D700}=0.5$
, 0.1,
$-0.1$
and
$-0.5$
.
Figure 6. The phase lag
$\unicode[STIX]{x1D711}$
between
$\unicode[STIX]{x1D719}$
and the wave plotted as a function of the eccentricity
$\unicode[STIX]{x1D700}$
, as computed from (3.41) for
$kA=0.2$
,
$H=10$
,
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}$
,
$x_{o}=0$
and
$z_{o}$
, chosen such that the particle always remained just below the free surface.
Finally, we can also compute the phase lag
$\unicode[STIX]{x1D711}$
between
$\unicode[STIX]{x1D719}$
and the wave, which is given by
$$\begin{eqnarray}\displaystyle \tan \unicode[STIX]{x1D711}=\frac{D}{C} & & \displaystyle\end{eqnarray}$$
and plotted in figure 6.
These results – that the particles settle onto a stable limit cycle in the plane of the waves about a preferential angle
$\overline{\unicode[STIX]{x1D719}^{\ast }}$
– hold for all particle and wave parameters. However, because our analysis assumes that
$kA$
is small, our predictions for
$A_{\unicode[STIX]{x1D719}}$
and
$\unicode[STIX]{x1D711}$
will break down for large values of
$kA$
.
4 Summary and discussion
We have considered here the dynamics of non-inertial point spheroids moving in the flow field generated by linear surface gravity waves. Numerically solving the equations of motion shows that the particles move on ellipsoidal orbits with a superimposed Stokes drift, just as one would expect, but that they adopt preferred orientations about which they oscillate as they move through the wave field. By analysing the equations of motion, we have shown that this preferential orientation is a consequence of finite wave amplitude, and appears in the angular dynamics just as Stokes drift appears in the translational dynamics at the same level of approximation. We analytically confirmed our previous result that the preferred polar angle of the particles is simply related to the aspect ratio (DiBenedetto et al. Reference DiBenedetto, Ouellette and Koseff2018), and additionally were able to predict the amplitude and phase of the oscillation about this preferred angle.
Our main results are summarized here. Asymptotically, the particles settle onto a limit cycle consisting of an oscillation with amplitude
$A_{\unicode[STIX]{x1D719}}$
and phase lag
$\unicode[STIX]{x1D711}$
about a preferential angle
$\overline{\unicode[STIX]{x1D719}^{\ast }}$
. These quantities are all independent of the particle’s initial orientation, which only controls how long the particles take to approach this limit cycle. The preferential angle
$\overline{\unicode[STIX]{x1D719}^{\ast }}$
is entirely determined by the particle shape and is independent of the wave parameters. The amplitude and phase of the oscillation, however, are functions of the wave parameters, the particle’s initial position within the wave field, and the particle shape. For easy reference, the limit cycle is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}(t)=A_{\unicode[STIX]{x1D719}}\sin (kx_{o}-\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D711})+\bar{\unicode[STIX]{x1D719}}^{\ast }, & & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D719}}^{\ast }=\tan ^{-1}\unicode[STIX]{x1D706}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{\unicode[STIX]{x1D719}}=\sqrt{C^{2}+D^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}=\tan ^{-1}\frac{D}{C}, & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle C & = & \displaystyle -\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\sinh [k(z_{o}+H)]\cos [2(D+\bar{\unicode[STIX]{x1D719}}^{\ast })],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle D & = & \displaystyle -\frac{kA\unicode[STIX]{x1D700}}{\sinh kH}\cosh [k(z_{o}+H)]\sin [2(C+\bar{\unicode[STIX]{x1D719}}^{\ast })].\end{eqnarray}$$
In the Introduction, we motivated this study by considering the problem of the transport of microplastic particles in the ocean. We have, however, analysed only a relatively idealized case here. Real microplastics, for example, are not true point particles, are not perfect spheroids, and have (weak) inertia relative to the flow since they are not perfectly neutrally buoyant (Ryan et al.
Reference Ryan, Moore, van Franeker and Moloney2009; Chubarenko et al.
Reference Chubarenko, Bagaev, Zobkov and Esiukova2016). Nevertheless, we argue that our work still has relevance to the more applied question of microplastic transport. Although microplastics are not point particles (indeed, they are typically defined as being up to 5 mm in linear size, although they can be as small as
$1~\unicode[STIX]{x03BC}\text{m}$
), they tend to have small Stokes and particle Reynolds numbers so that a point-particle approximation is reasonable. They are often close to spheroidal in shape. And although they are not neutrally buoyant, the density of most plastics is close to that of seawater, so that microplastics are typically only slightly negatively or positively buoyant. Thus, we would expect that adding these complications will not qualitatively change the behaviour of the particles from the results we have found here, although there may be some quantitative changes. Even those, however, are likely to be small. Our previous numerical study, for example, included the effects of weak particle inertia and found very similar results to what we show here (DiBenedetto et al.
Reference DiBenedetto, Ouellette and Koseff2018). Our results here demonstrate that preferential orientation arises from the fundamental interaction between particle shape and wave motion, and is not a consequence of particle inertia. Additionally, it is known that Stokes drift is important for the translational dynamics of real microplastics (Isobe et al.
Reference Isobe, Kubo, Tamura, Kako, Nakashima and Fujii2014; van den Bremer & Breivik Reference van den Bremer and Breivik2017); and since we have shown that the preferred orientation of the particles arises from the same physics, it stands to reason that it also ought to persist.
Thus, we expect that using the results we have derived here to model the more realistic case will be a good approximation, and certainly superior to simply assuming that the particles are point-like non-inertial spheres (Kukulka et al. Reference Kukulka, Proskurowski, Morét-Ferguson, Meyer and Law2012; Maximenko, Hafner & Niiler Reference Maximenko, Hafner and Niiler2012). In particular, it would be interesting to use our results to estimate the lift and drag felt by weakly inertial spheroids, since the preferred orientation would imply an angle of attack relative to the flow that depends on particle shape in a way that is fundamentally different from the spherical case, as a first step to developing more accurate models of microplastic dispersion in the ocean.
Acknowledgements
We thank J. R. Koseff for helpful discussions over the course of this work. M.H.D. gratefully acknowledges the support of the Stanford Gerald J. Lieberman Fellowship. This work was supported by the US National Science Foundation under grant no. CBET-1706586.
Appendix A. Stability of fixed points
Here, we revisit the question of the stability of the
$p_{y}$
fixed points in the dynamical system given by (2.6a
)–(2.6c
). For convenience, these equations are
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{x} & = & \displaystyle S_{xx}p_{x}(1-(p_{x}^{2}-p_{z}^{2}))+S_{xz}p_{z}(1-2p_{x}^{2}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{y} & = & \displaystyle -p_{y}[S_{xx}(p_{x}^{2}-p_{z}^{2})+2S_{xz}p_{x}p_{z}],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}^{-1}{\dot{p}}_{z} & = & \displaystyle S_{xz}p_{x}(1-2p_{z}^{2})-S_{xx}p_{z}(1+(p_{x}^{2}-p_{z}^{2})).\end{eqnarray}$$
There are two fixed points for the
$p_{y}$
equation:
$p_{y}=1$
(which requires
$p_{x}=p_{z}=0$
), and
$p_{y}=0$
. For the first of these, we can use simple linear stability analysis to determine its stability. Evaluating the Jacobian
$\unicode[STIX]{x1D645}$
of the dynamical system for
$p_{y}=1$
and
$p_{x}=p_{z}=0$
, we find
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D645}=\unicode[STIX]{x1D716}\left(\begin{array}{@{}ccc@{}}S_{xx} & 0 & S_{xz}\\ 0 & 0 & 0\\ S_{xz} & 0 & -S_{xx}\end{array}\right). & & \displaystyle\end{eqnarray}$$
The eigenvalues of this matrix are 0 and
$\pm \sqrt{S_{xx}^{2}+S_{xz}^{2}}$
. Thus, since one of these eigenvalues is always real and positive, this fixed point is unstable.
For
$p_{y}=0$
, the situation is more complicated. It is not a fixed point for the full dynamical system, and the eigenvalues of the Jacobian evaluated at
$p_{y}=0$
are not sign definite (as one would expect, given that we have argued above that the particle dynamics in the wave plane settle onto a limit cycle rather than a fixed point). To make progress, it is convenient to switch to spherical coordinates. Using the polar angle
$\unicode[STIX]{x1D719}$
and the azimuthal angle
$\unicode[STIX]{x1D703}$
, we have
$$\begin{eqnarray}\displaystyle p_{x} & = & \displaystyle \sin \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D703},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{y} & = & \displaystyle \cos \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D703},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{z} & = & \displaystyle \cos \unicode[STIX]{x1D719}.\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}^{-1}\dot{\unicode[STIX]{x1D703}} & = & \displaystyle -{\textstyle \frac{1}{2}}S_{xx}\sin 2\unicode[STIX]{x1D703}-S_{xz}\cot \unicode[STIX]{x1D719}\sin \unicode[STIX]{x1D703},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D700}^{-1}\dot{\unicode[STIX]{x1D719}} & = & \displaystyle {\textstyle \frac{1}{2}}S_{xx}\sin 2\unicode[STIX]{x1D719}(1+\cos ^{2}\unicode[STIX]{x1D703})+S_{xz}\cos 2\unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$
$p_{y}=0$
corresponds to the family of points
$\unicode[STIX]{x1D703}=0+n\unicode[STIX]{x03C0}$
for integer
$n$
. To evaluate the stability of these points for any
$\unicode[STIX]{x1D719}$
, we first compute 
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}=-\unicode[STIX]{x1D700}\left[S_{xx}\cos 2\unicode[STIX]{x1D703}+S_{x}z\cot \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D703}-S_{xz}\csc ^{2}\unicode[STIX]{x1D719}\sin \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right].\end{eqnarray}$$
Evaluating this expression at
$\unicode[STIX]{x1D703}=0+n\unicode[STIX]{x03C0}$
, we have
$$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right|_{\unicode[STIX]{x1D703}=0+2\unicode[STIX]{x03C0}n} & = & \displaystyle -\unicode[STIX]{x1D700}(S_{xx}+S_{xz}\cot \unicode[STIX]{x1D719}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right|_{\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}+2\unicode[STIX]{x03C0}n} & = & \displaystyle -\unicode[STIX]{x1D700}(S_{xx}-S_{xz}\cot \unicode[STIX]{x1D719}).\end{eqnarray}$$
$S_{xx}$
and
$S_{xz}$
are purely oscillatory with fixed amplitude (with
$S_{xx}$
oscillating about
$0$
and
$S_{xz}$
oscillating about a non-zero positive value). As we have shown above, if we assume that the spheroid is aligned in the plane of the waves,
$\unicode[STIX]{x1D719}$
settles onto a limit cycle that oscillates at the same frequency as the waves with constant amplitude about a fixed value. Thus, to evaluate the stability of the in-plane fixed point, we need only consider the behaviour of the envelope about these oscillations, and can therefore average equations (A 6) over a wave period. Doing so leads to the new equations 
$$\begin{eqnarray}\displaystyle \left.\overline{\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}}\right|_{\unicode[STIX]{x1D703}=0+2\unicode[STIX]{x03C0}n} & = & \displaystyle -\unicode[STIX]{x1D700}\overline{S_{xz}\cot \unicode[STIX]{x1D719}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \left.\overline{\frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}}\right|_{\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}+2\unicode[STIX]{x03C0}n} & = & \displaystyle \unicode[STIX]{x1D700}\overline{S_{xz}\cot \unicode[STIX]{x1D719}}.\end{eqnarray}$$
$\unicode[STIX]{x1D719}$
tends asymptotically to a limit cycle,
$\overline{S_{xz}\cot \unicode[STIX]{x1D719}}$
will be constant in time. Thus, at least one of the two cases in (A 7) will be stable – and since they correspond to the same physical case (given the symmetries of the problem), the set defined by
$p_{y}=0$
is thus stable as well.
