2. Set-up and governing equations

We begin by introducing Cartesian axes, $Oxyz'$, describing a frame rotating around the vertical ($z'$) axis with angular velocity $\varOmega$. We consider the evolution of a spherical vortex of (spherical) radius $a(t)$ moving with speed $U(t)$ in the $z'$ direction so we define a new vertical coordinate

(2.1)\begin{equation} z = z'-\int_0^t U(t') \textrm{d} t', \end{equation}

such that the coordinate origin moves with the vortex. This transformation is only applied to coordinates and not to the associated velocity fields. Our non-dimensional system is described by the rotating, incompressible Euler equations

(2.2a)$$\begin{gather} \left[\frac{\partial }{\partial t}-U\frac{\partial }{\partial z}\right] \boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} + \epsilon \hat{\boldsymbol{z}}\times\boldsymbol{u} ={-}\boldsymbol{\nabla} p, \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

where $\boldsymbol{u}$ and p describe the fluid velocity and pressure respectively. We note that all velocities are relative to the background fluid and we do not gain any fictitious forces beyond the Coriolis force (Johnson & Crowe Reference Johnson and Crowe2021). Here, we have non-dimensionalised velocities by $V$, lengths by $L$, time by $L/V$ and pressure by $\rho _0 V^2$ for density $\rho _0$. The length scale $L$ and velocity scale $V$ are set by the initial size and translational velocity of the vortex. The strength of the rotation is represented in our non-dimensional equations by the inverse Rossby number, defined as

(2.3)\begin{equation} \epsilon = \frac{1}{{Ro}} = \frac{2\varOmega L}{V}, \end{equation}

which will be taken to be a small parameter. A diagram of our system is given in figure 1.

Figure 1. Diagram of our set-up showing Cartesian coordinates, $(x,y,z)$, cylindrical radius, $\rho$, azimuthal angle, $\phi$, spherical radius, $r$ and polar angle, $\theta$. A spherical vortex of radius $a(t)$ moves vertically with velocity $U(t)$ and the system is rotating with angular velocity of $\varOmega$ corresponding to an inverse Rossby number of $\epsilon$.

We now take the coordinate origin to be the centre of the spherical vortex and introduce cylindrical coordinates $(\rho ,\phi ,z)$ where the cylindrical radius is given by $\rho ^2 = x^2 + y^2$ and the azimuthal angle $\tan \phi = y/x$. Under the assumption that the flow is independent of the azimuthal angle, $\phi$, the velocity field may be written as

(2.4)\begin{equation} \boldsymbol{u}(\rho,z) = u_\rho\hat{\boldsymbol{\rho}} + w \hat{\boldsymbol{\phi}} + u_z\hat{\boldsymbol{z}} ={-}\frac{1}{\rho}\frac{\partial \psi}{\partial z}\hat{\boldsymbol{\rho}} + w \hat{\boldsymbol{\phi}} + \frac{1}{\rho} \frac{\partial \psi}{\partial \rho}\hat{\boldsymbol{z}}, \end{equation}

for Stokes streamfunction $\psi$ and swirl velocity $w$. We will later use spherical polar coordinates $(r,\theta ,\phi )$ when solving for the vortex. Here the spherical radius $r$ is defined by $r^2 = \rho ^2 + z^2$ and the polar angle by $\tan \theta = \rho /z$.

As we will show later, there are two relevant time scales in the problem. The first is the transient time scale on which the wave field develops and the second is the slow evolution time scale on which the vortex loses energy. We therefore expand our time derivative using a multiple scales approach to write

(2.5)\begin{equation} \frac{\partial }{\partial t} \to \frac{\partial }{\partial t} + \delta \frac{\partial }{\partial T}, \end{equation}

where $T = \delta t$ is the slow evolution time scale and $\delta$ is a small parameter. It will be shown later that $\delta = \epsilon ^4$. Finally we will ignore the problem of the transient evolution and look for solutions with a fully developed wave field. Therefore, we take ${\partial }/{\partial t} = 0$ and consider only the slow evolution on time scale $T$. Since the vortex speed and radius may slowly evolve as energy is lost we take $U = U(T)$ and $a = a(T)$.

2.1. The vorticity equation

Taking the curl of (2.2a) and neglecting the transient evolution gives the vorticity equation

(2.6)\begin{equation} \left[\delta\frac{\partial }{\partial T} -U\frac{\partial }{\partial z}\right]\boldsymbol{\omega} + \boldsymbol{\nabla}\times(\boldsymbol{\omega}\times\boldsymbol{u}) - \epsilon \frac{\partial \boldsymbol{u}}{\partial z} = 0, \end{equation}

for vorticity

(2.7)\begin{equation} \boldsymbol{\omega} ={-}\frac{1}{\rho}\frac{\partial }{\partial z}(\rho w) \hat{\boldsymbol{\rho}} - \frac{1}{\rho} D^2\psi \hat{\boldsymbol{\phi}} + \frac{1}{\rho}\frac{\partial }{\partial \rho}(\rho w)\hat{\boldsymbol{z}}. \end{equation}

Here, $D^2$ is related to the Laplacian operator and is given by

(2.8)\begin{equation} D^2f = \nabla^2f - \frac{2}{\rho}\frac{\partial f}{\partial \rho} = \rho\frac{\partial }{\partial \rho}\left(\frac{1}{\rho}\frac{\partial f}{\partial \rho}\right) + \frac{\partial ^2f}{\partial z^2}, \end{equation}

for some field f. Substituting (2.4) and (2.7) into (2.6) gives two equations for $\psi$ and $w$ from the three vorticity components. Both the $\rho$ and $z$ components of the vorticity equation give

(2.9)\begin{equation} \left[\delta\frac{\partial }{\partial T} -U\frac{\partial }{\partial z}\right](\rho w) + \frac{1}{\rho} J[\psi,\rho w] - \epsilon \frac{\partial \psi}{\partial z} = 0, \end{equation}

while the $\phi$ component gives

(2.10)\begin{equation} \left[\delta\frac{\partial }{\partial T} -U\frac{\partial }{\partial z}\right]\left[ \frac{1}{\rho} D^2 \psi \right] + J\left[ \psi,\frac{1}{\rho^2}D^2 \psi \right] + \frac{1}{\rho^3}\frac{\partial }{\partial z}(\rho w)^2 + \epsilon \frac{\partial w}{\partial z} = 0. \end{equation}

The Jacobian operator $J$ is taken to be between $\rho$ and $z$ hence

(2.11)\begin{equation} J[f,g] = \frac{\partial (\,f,g)}{\partial (\rho,z)} = \frac{\partial f}{\partial \rho}\frac{\partial g}{\partial z} - \frac{\partial f}{\partial z}\frac{\partial g}{\partial \rho}, \end{equation}

for some fields f and g.

The rotation terms in (2.9) and (2.10) may now be combined with the swirl velocity and the $U\partial _z$ terms may be written using the Jacobian derivative to give

(2.12)\begin{equation} \delta\frac{\partial }{\partial T}\varGamma + \frac{1}{\rho} J\left[\psi-\frac{U}{2}\rho^2,\varGamma\right] = 0, \end{equation}

and

(2.13)\begin{equation} \delta\frac{\partial }{\partial T} \left[ \frac{1}{\rho} D^2 \psi \right] + J\left[ \psi-\frac{U}{2}\rho^2,\frac{1}{\rho^2}D^2 \psi \right] + \frac{1}{\rho^3}\frac{\partial }{\partial z}\varGamma^2 = 0, \end{equation}

where

(2.14)\begin{equation} \varGamma = \rho w+\frac{\epsilon\rho^2}{2}, \end{equation}

is the total azimuthal circulation divided by $2{\rm \pi}$. Equations (2.12) and (2.13) are our governing equations which will be used throughout this analysis. On the $z$ axis ($\rho = 0$) we impose the condition that $w = 0$ and $\psi = 0$ such that all components of velocity and vorticity are finite and continuous.

2.2. Asymptotic approach

To proceed we will assume that $\delta$ is small and expand $\psi$ and $w$ using a regular perturbation expansion of the form

(2.15)\begin{equation} (\psi,w) = (\psi_0,w_0) + \delta (\psi_1,w_1) +O(\delta^2), \end{equation}

where each term in the expansion is a function of $\epsilon$. We will later show that $\epsilon$ can be related to $\delta$ through $\epsilon ^4 = \delta$ and at this point we could expand each term in powers of $\epsilon = \delta ^{1/4}$ before reordering terms such that the expansion is asymptotically well-ordered in $\delta$. Anticipating the relation $\epsilon ^4 = \delta$, we therefore ignore any components which are $O(\epsilon ^4)$ or higher when calculating the leading-order solution in $\delta$.

To leading order in $\delta$ (2.12) and (2.13) give that

(2.16)\begin{equation} \frac{1}{\rho} J\left[\psi_0-\frac{U}{2}\rho^2,\rho w_0 + \frac{\epsilon\rho^2}{2}\right]= 0, \end{equation}

and

(2.17)\begin{equation} J\left[ \psi_0-\frac{U}{2}\rho^2,\frac{1}{\rho^2}D^2 \psi_0 \right] + \frac{1}{\rho^3}\frac{\partial }{\partial z}\left(\rho w_0+\frac{\epsilon\rho^2}{2}\right)^2 = 0. \end{equation}

We begin by using (2.16) to relate the swirl velocity, $w$, to $\psi$.

2.3. The swirl velocity and Kelvin's circulation theorem

Equation (2.12) states that the azimuthal circulation is conserved following a ring of particles of radius $\rho$, where the term $\epsilon \rho ^2/2$ describes the circulation of the rotating frame. This can be thought of as a consequence of Noether's theorem: axial angular momentum in conserved due to the invariance of the system to azimuthal rotation. Consider an incoming ring of particles that far from the vortex is arbitrarily close to the $z$ axis, i.e. $\rho \ll 1$, so $\varGamma$ is vanishingly small. The circulation $\varGamma$ is conserved as the ring of particles expands to pass around the vortex and so $\varGamma$ remains arbitrarily small arbitrarily close to the vortex with $\varGamma$ vanishing on the vortex boundary by continuity. Thus

(2.18)\begin{equation} w_0 ={-}\frac{\epsilon\rho}{2}\quad \text{at}\ r = a. \end{equation}

Equation (2.16) implies that

(2.19)\begin{equation} \rho w_0 + \frac{\epsilon\rho^2}{2} = K\left(\psi_0-\frac{U}{2}\rho^2\right), \end{equation}

for some function $K = K(\varphi )$. This result is simply Kelvin's circulation theorem for a ring of particles an axisymmetric flow. We now take

(2.20)\begin{equation} \psi_0 = \frac{U}{2}\rho^2\quad \text{at}\ r = a, \end{equation}

and use (2.18) to impose the condition that $K(0) = 0$. We note that (2.20) is equivalent to the requirement that the vortex boundary, $r = a$, is a streamsurface in a frame moving with the vortex.

For the fluid outside the spherical vortex, all streamlines must originate from the ultra far field where $\psi _0 \to 0$ and $w_0 \to 0$. Therefore (2.19) gives that

(2.21)\begin{equation} \frac{1}{2}\epsilon\rho^2 = K\left(-\frac{U}{2}\rho^2\right), \end{equation}

hence, outside the vortex, we can determine $K$ as the linear function

(2.22)\begin{equation} K(\varphi) ={-}\epsilon k \varphi, \end{equation}

for some field $\varphi$. Here, we have defined $k(T) = 1/U(T)$ to be an effective wavenumber.

Inside the vortex (2.19) remains valid, however, as streamlines do not necessarily originate from the far field, we are unable to determine $K$ using a far-field condition. The choice of this function $K$ is arbitrary and picking a general linear function, $K(\varphi ) = \alpha \varphi$, for some parameter $\alpha$ results in the family of vortices considered by Moffatt (Reference Moffatt1969). To ensure that our vortex reduces to Hill's vortex as $\epsilon \to 0$, we assume that (2.22) continues to hold within the vortex. This assumption corresponds to a vortex with a slow swirl and ensures that ${\partial w_0}/{\partial r}$ is continuous across $r = a$. The value of $\alpha = -\epsilon k$ can be thought of as an internal wavenumber describing the azimuthal rotation rate of the vortex.

2.4. The streamfunction

Using (2.19), (2.17) becomes

(2.23)\begin{equation} J\left[ \psi_0-\frac{U}{2}\rho^2,\frac{1}{\rho^2}D^2 \psi_0 \right] + \frac{1}{\rho^3}\frac{\partial }{\partial z}\left[K\left(\psi_0-\frac{U}{2}\rho^2\right)\right]^2 = 0, \end{equation}

and since it can be shown that

(2.24)\begin{equation} J\left[ \varphi, \frac{1}{\rho^2} K'\left(\varphi\right) K\left(\varphi\right) \right] = \frac{1}{\rho^3} \frac{\partial }{\partial z}\left[K\left(\varphi\right)\right]^2, \end{equation}

for some field $\varphi$, we can write

(2.25)\begin{equation} J\left[ \psi_0-\frac{U}{2}\rho^2,\frac{1}{\rho^2}\left\{D^2 \psi_0 +K'\left(\psi_0-\frac{U}{2}\rho^2\right) K\left(\psi_0-\frac{U}{2}\rho^2\right) \right\}\right] = 0. \end{equation}

Therefore, we have

(2.26)\begin{equation} D^2 \psi_0 +K'\left(\psi_0-\frac{U}{2}\rho^2\right) K\left(\psi_0-\frac{U}{2}\rho^2\right) = \rho^2 Z\left(\psi_0-\frac{U}{2}\rho^2\right), \end{equation}

for some function $Z(\varphi )$. We note that the function $Z$ can be related to the pressure head $H(\psi _0-U\rho ^2/2) = p_0 + |\boldsymbol {u}_0|^2/2$ using Bernoulli's principle (Batchelor Reference Batchelor2000, § 7.5) by $Z = H'$.

We can now determine $Z$ outside the vortex using the same method as was used to determine $K$. In the far field $D^2\psi _0 \to 0$, hence we have

(2.27)\begin{equation} -\frac{1}{2} \epsilon^2 k^2 U = Z\left(-\frac{U}{2}\rho^2\right), \end{equation}

so $Z$ is the constant function

(2.28)\begin{equation} Z ={-}\tfrac{1}{2} U \epsilon^2 k^2, \end{equation}

where we recall that $k U = 1$. For the flow inside the vortex we follow Moffatt (Reference Moffatt1969) and take

(2.29)\begin{equation} H(\varphi) = H_0 + \lambda \varphi \quad \implies\quad Z(\varphi) = \lambda, \end{equation}

for some constant $\lambda$. Therefore, our leading-order streamfunction satisfies

(2.30)\begin{equation} \left[ D^2 + \epsilon^2 k^2 \right] \psi_0 = \begin{cases} 0, \quad & r > a,\\ \left(\lambda+\frac{1}{2}U\epsilon^2 k^2 \rho^2\right) \rho^2, & r < a, \end{cases} \end{equation}

such that $\psi _0 = U\rho ^2/2$ on $r = a$ and ${\partial \psi _0}/{\partial r}$ is continuous across $r = a$ so the velocity is continuous across the vortex boundary. These boundary conditions will allow us to uniquely determine $\lambda$. Once we have solved for $\psi _0$, we can find $w_0$ from $\psi _0$ using (2.19).

We observe from (2.30) that there are two natural length scales in the problem – the vortex radius, $a$, and the wavelength $(\epsilon k)^{-1}$. We therefore expect our solution to consist of an interior vortex region with $\epsilon k r \ll 1$ and an exterior wave field where $\epsilon k r \sim 1$. We begin by determining the interior vortex solution and then determine the exterior wave field, the wave amplitude is obtained by matching to the interior solution.

3. Interior vortex solution

Here, we solve for the leading-order streamfunction, $\psi _0$, in the interior region around the vortex where $\epsilon k r \ll 1$.

3.1. The flow inside the vortex

Inside the vortex, where $r\leq a$, we have

(3.1)\begin{equation} \left[D^2+\epsilon^2k^2\right] \psi_0 = \left(\lambda+\frac{1}{2}U\epsilon^2 k^2 \rho^2\right) \rho^2, \end{equation}

from (2.30), with general solution

(3.2)\begin{align} \psi_0 = \left[\frac{\lambda}{\epsilon^2k^2}+\frac{U}{2}\right]r^2 \sin^2 \theta + r \sin \theta \sum_{m=1}^\infty \textrm{P}^1_m(\cos \theta) \left[ A_m \textrm{j}_m(\epsilon k r) +B_m \textrm{y}_m(\epsilon k r)\right],\quad r \leq a. \end{align}

Here, the $\textrm {P}^1_n$ are associated Legendre polynomials and the $\textrm {j}_n$ and $\textrm {y}_n$ are spherical Bessel functions of the first and second kind respectively and we have imposed the conditions that our solution is $2{\rm \pi}$ periodic and $\theta = 0,{\rm \pi}$ are streamlines. Since we require $\psi _0$ to be non-singular at the origin we must have $B_m = 0$ for all $m$ and the condition that $\psi _0 = U\rho ^2/2$ at $r = a$ gives that $A_m = 0$ for $m \geq 2$ and

(3.3)\begin{equation} A_1 ={-}\frac{\lambda a}{\epsilon^2 k^2 \textrm{j}_1(\epsilon k a)}. \end{equation}

Therefore, our solution is

(3.4)\begin{equation} \psi_0 = \left(\frac{U}{2} + \frac{\lambda}{\epsilon^2 k^2}\left[1 -\frac{a \textrm{j}_1(\epsilon k r)}{r \textrm{j}_1(\epsilon k a)} \right]\right) r^2 \sin^2\theta,\quad r \leq a, \end{equation}

where the value of $\lambda$ is determined by the continuity of polar velocity, $u_\theta$, across $r = a$.

3.2. The flow around the vortex

Outside the vortex, for $r \geq a$, we have

(3.5)\begin{equation} \left[ D^2 + \epsilon^2k^2 \right] \psi_0 = 0, \end{equation}

from (2.30) with general solution

(3.6)\begin{equation} \psi_0 = r \sin \theta \sum_{m=1}^\infty \textrm{P}^1_m(\cos \theta) \left[ C_m \textrm{j}_m(\epsilon k r) +D_m \textrm{y}_n(\epsilon k r)\right],\quad r \geq a. \end{equation}

Requiring that $\psi _0 = U\rho ^2/2$ on $r = a$ and our solution decays away from the origin for $\epsilon k r \ll 1$ gives that

(3.7)\begin{equation} \psi_0 = \frac{Ua}{2}\frac{\textrm{y}_1(\epsilon k r)}{\textrm{y}_1(\epsilon k a)}r \sin^2\theta,\quad r \geq a. \end{equation}

3.3. The full interior solution

Our final condition on $\psi _0$ is that the polar velocity, $u_{\theta 0}$, is continuous across $r = a$ which gives that ${\partial \psi _0}/{\partial r}$ must be continuous. Using (3.7) and (3.4) we have

(3.8)\begin{equation} \frac{\lambda a}{\epsilon^2 k^2} \frac{\epsilon k a \textrm{j}_2(\epsilon k a)}{\textrm{j}_1(\epsilon k a)} ={-}\frac{U a}{2} \frac{\epsilon k a \textrm{y}_2(\epsilon k a)}{\textrm{y}_1(\epsilon k a)}, \end{equation}

so

(3.9)\begin{equation} \lambda ={-}\frac{15U}{2a^2}\left( 1-\frac{38}{105}(\epsilon k a)^2 + O(\epsilon^4)\right). \end{equation}

Noting that $\epsilon k r \ll 1$ for $r \leq a$ and neglecting terms of order $O(\epsilon ^4)$ and higher, we can expand (3.4) as

(3.10) \begin{align} \psi_0 &= \frac{U}{2}\left( 1+ \left[\frac{3}{2}-\frac{19}{35}(\epsilon k a)^2\right] \left[1-\frac{r^2}{a^2}\right] +\frac{3(\epsilon k a)^2}{280}\left[ 9-\frac{14r^2}{a^2}+\frac{5r^4}{a^4}\right]\right) r^2 \sin^2\theta,\notag\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad r\leq a. \end{align}

In the limit of $\epsilon k a \to 0$ this result reduces to the classical vortex solution of Hill (Reference Hill1894) given by

(3.11)\begin{equation} \psi_0 = \frac{3U}{4}\left[\frac{5}{3}-\frac{r^2}{a^2}\right] r^2 \sin^2\theta, \quad r\leq a. \end{equation}

Similarly, expanding (3.7) for $\epsilon k r \ll 1$ and neglecting terms of order $O(\epsilon ^4)$ or higher we obtain

(3.12)\begin{equation} \psi_0 = \frac{U a^3} {2}\left[ 1 + \frac{(\epsilon k a)^2}{2}\left( \frac{r^2}{a^2}-1\right) \right] \frac{\sin^2\theta}{r},\quad r \geq a, \end{equation}

which reduces to inviscid flow around a sphere in the case of $\epsilon ka \to 0$.

4. Exterior wave field

We now determine the wave field in the exterior region where $\epsilon k r \sim 1$. From (2.30) we have

(4.1)\begin{equation} \left[ D^2 + \epsilon^2 k^2 \right] \psi_0 = 0, \end{equation}

for $r > a$. This equation may alternatively be derived from the linearised Euler equations under the assumption that far from the vortex the wave wake is small in amplitude (Johnson & Crowe Reference Johnson and Crowe2021) showing that the linear wave field satisfies the nonlinear system. Since solutions to (4.1) will vary over a length scale of order $O(\epsilon ^{-1})$ we now balance the two terms in this equation by introducing the long radial scale $R = \epsilon r$. Switching to spherical coordinates gives

(4.2)\begin{equation} \frac{\partial ^2\psi_0}{\partial R^2} + \frac{\sin \theta}{R^2}\frac{\partial }{\partial \theta}\left[ \frac{1}{\sin \theta}\frac{\partial \psi_0}{\partial \theta}\right] + k^2 \psi_0 = 0, \end{equation}

which can be solved using separation of variables for solution

(4.3)\begin{equation} \psi_0 = kR \sum_{n=1}^\infty \sin \theta \textrm{P}^1_n(\cos \theta) \left[ A_n \textrm{j}_n(kR) +B_n \textrm{y}_n(kR)\right]. \end{equation}

Long (Reference Long1953) observed that waves are not seen ahead of objects moving through a rotating flow. These observations are supported by the theoretical work of McIntyre (Reference McIntyre1972) for the case of an unbounded rotating flow motivating us to use the radiation condition of no upstream waves. Fraenkel (Reference Fraenkel1956) showed that this condition can be satisfied by $B_n = 0$ for $n \geq 2$ and

(4.4)\begin{equation} A_n = \begin{cases} 0, & n\text{ odd},\\ -\dfrac{(2n+1)(n-2)!}{2^n\left(\dfrac{n}{2}+1\right)!\left(\dfrac{n}{2}-1\right)!}B_1, & n\text{ even}, \end{cases} \end{equation}

so our solution may be written as

(4.5)\begin{equation} \psi_0 = BR \sin\theta \left[\textrm{y}_1(kR)\sin\theta + \sum_{n=1}^\infty \frac{(4n+1)(2n-2)!}{2^{2n}\left(n+1\right)!\left(n-1\right)!} \textrm{j}_{2n}(kR) \textrm{P}^1_{2n}(\cos \theta) \right], \end{equation}

where $B = -kB_1$ can be found by matching to the interior solution. We note that spherical Bessel functions take on a trigonometric form for $k R = O(1)$ justifying the definition of $k = 1/U$ as a wavenumber. Substituting $r = R/\epsilon$, we can see that wavelength of these waves is $O(1/\epsilon )$ assuming that $U$ remains $O(1)$.

In order to find the form of $B$ we match (4.5) to the inner solution in (3.7) by substituting $R = \epsilon r$. Noting that

(4.6a,b)\begin{equation} \textrm{j}_n(x) \sim \frac{x^n}{(2n+1)!!} \quad \text{and} \quad \textrm{y}_1(x) \sim{-}\frac{1}{x^{2}}, \end{equation}

for small argument $x$ we have

(4.7)\begin{equation} \psi_0 \sim \epsilon B \textrm{y}_1(\epsilon k r) r \sin^2\theta + O(\epsilon^4). \end{equation}

Hence, by comparison with (3.7), we have

(4.8)\begin{equation} B = \frac{Ua}{2\epsilon \textrm{y}_1(\epsilon k a)} ={-}\frac{\epsilon U a^3 k^2}{2}+O(\epsilon^3), \end{equation}

to leading order in $\delta$. We note that close to the vortex ($R\ll 1$), $\psi _0$ is dominated by the $\textrm {y}_1(kR)$ term. This behaviour is consistent with the expected dipole structure near the vortex. Additionally, since $B = O(\epsilon )$, the wave amplitude is $O(\epsilon )$ and hence small for $R = O(1)$.

5. Summary of leading-order solution

Before moving on to the solvability conditions required to determine the slow time evolution of the vortex, we present the full leading-order solution. The interior Hill's vortex solution is given by

(5.1) \begin{align} \psi_0 = \begin{cases} \dfrac{U}{2}\left( 1+ \left[\dfrac{3}{2}-\dfrac{19}{35}(\epsilon k a)^2\right] \left[1-\dfrac{r^2}{a^2}\right] +\dfrac{3(\epsilon k a)^2}{280}\left[ 9-\dfrac{14r^2}{a^2}+\dfrac{5r^4}{a^4}\right]\right) r^2 \sin^2\theta,\\ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad r \leq a,\\ \dfrac{U a^3} {2}\left( 1 + \dfrac{(\epsilon k a)^2}{2}\left[ \dfrac{r^2}{a^2}-1\right] \right) \dfrac{\sin^2\theta}{r}, \quad r \geq a, \end{cases} \end{align}

for $\epsilon k r \ll 1$. The terms of order $O((\epsilon k a)^2)$ in this result correspond to both the near-field inertial waves and the modification of internal vortex structure by the rotation. Far from the vortex, the exterior wave field streamfunction is

(5.2)\begin{equation} \psi_0 = \frac{a U R \sin\theta}{2\epsilon \textrm{y}_1(\epsilon k a)} \left[\textrm{y}_1(kR)\sin\theta + \sum_{n=1}^\infty \frac{(4n+1)(2n-2)!}{2^{2n}\left(n+1\right)!\left(n-1\right)!} \textrm{j}_{2n}(kR) \textrm{P}^1_{2n}(\cos \theta) \right], \end{equation}

for long radial scale $R = \epsilon r$. We note that the amplitude of the wave field is $O(\epsilon )$ since $\textrm {y}_1(\epsilon k a) = O(\epsilon ^{-2})$.

Figure 2(a) shows the streamfunction, $\psi _0$, for the interior vortex solution from (5.1) for $U = 1$, $a = 1$ and $\epsilon = 0.2$. The black lines denote streamlines in a frame moving with the vortex – corresponding to streamlines of $\psi '_0 = \psi _0 - U\rho ^2/2$. The blue line denotes the streamline $\psi '_0 = 0$ so corresponds to the vortex boundary, $r = a$. The streamlines of $\psi '_0$ describe the particle paths of incoming fluid over the short time scale, $t$, and will vary over the slow scale, $T$, as the vortex speed and radius change. Figure 2(b) shows the streamfunction, $\psi _0$, for the exterior wave field given by (5.2). We plot $\psi _0/|B|$ (for $B$ defined in (4.8)) as a function of the long spatial scales $(\epsilon k\rho ,\epsilon k z) = kR(\sin \theta ,\cos \theta )$ so the plotted result is unchanged for all values of $\epsilon$, $U$ and $a$. We observe a wave wake in the region $z < 0$ with decaying upstream perturbation. The wave amplitude is maximal along the line $\theta =2{\rm \pi} /3$ as discussed in Appendix B.

Figure 2. (a) Plot of $\psi _0$ from (5.1) for $U = 1$, $a = 1$ and $\epsilon = 02$. The black lines denote streamlines of $\psi '_0 = \psi _0 - U\rho ^2/2$ and the blue line denotes the vortex boundary, $r = a$. (b) Plot of $\psi _0/|B|$ from (5.2) as a function of the long spatial scales $(\epsilon k\rho ,\epsilon k z) = kR(\sin \theta ,\cos \theta )$.

6. Solvability conditions

We now consider the energy balance and azimuthal vorticity to determine solvability conditions which will allow us to examine the long-time decay of the vortex.

6.1. Conservation of energy

The leading-order kinetic energy of the vortex is given by

(6.1)\begin{equation} E_0 = \frac{1}{2} \int \left(u_{\rho 0}^2+u_{z 0}^2+w_0^2\right) \textrm{d} V = \int \frac{1}{2\rho^2}\left(|\boldsymbol{\nabla} \psi_0|^2 + \epsilon^2 k^2 \psi_0^2\right) \textrm{d} V, \end{equation}

where the integral is taken over the interior region. The generation of the wave field removes energy from the vortex leading to vortex decay with the loss of vortex energy balancing the outward flux of wave energy.

We can show (see Appendix A) that the total energy flux is given by

(6.2)\begin{equation} \delta\frac{\textrm{d}E_0}{\textrm{d}T} = U \int \left\{\frac{1}{2\rho^2}\left[\left( \frac{\partial \psi_0}{\partial z} \right)^2-\left( \frac{\partial \psi_0}{\partial \rho} \right)^2 \right] \hat{\boldsymbol{z}}\boldsymbol{\cdot}\hat{\boldsymbol{n}} + \frac{1}{\rho^2}\frac{\partial \psi_0}{\partial \rho}\frac{\partial \psi_0}{\partial z}\hat{\boldsymbol{\rho}}\boldsymbol{\cdot}\hat{\boldsymbol{n}} \right\}\textrm{d} S, \end{equation}

where the right-hand integral is evaluated over the surface of a sphere of radius $r = R/\epsilon$ and is independent of $R$. The required expression for $\psi _0$ on the sphere is given by the far-field streamfunction, (5.2), however, as the wave energy flux is $R$ independent we can use the asymptotic form of (5.2) for small $R$ given by

(6.3)\begin{equation} \psi_0 = \frac{\epsilon U a^3}{2}\left[\frac{\sin^2\theta}{R}+\frac{k^4 R^3 \sin^2 \theta \cos\theta}{8}\right]+O(R^4), \end{equation}

to simplify calculation of (6.2). Evaluation of (6.2) gives that

(6.4)\begin{equation} \delta\frac{\textrm{d}E_0}{\textrm{d}T} ={-}\frac{{\rm \pi}\epsilon^4 k^4 U^3 a^6}{4} + O(R^4), \end{equation}

where the $O(R^4)$ term may be set to zero as the result is independent of $R$. Noting that $k = 1/U$, defining $\delta = \epsilon ^4$ and calculating $E_0$ as the leading-order kinetic energy of a Hill's vortex,

(6.5)\begin{equation} E_0 = \frac{10{\rm \pi} U^2a^3}{7}, \end{equation}

our final result for the loss of vortex energy to the wave field is

(6.6)\begin{equation} \frac{\textrm{d}}{\textrm{d}T}\left[\frac{10{\rm \pi} U^2a^3}{7}\right] ={-}\frac{{\rm \pi} a^6}{4 U}. \end{equation}

6.2. Conservation of centre vorticity

As the vortex decays, we expect vortical fluid to escape and hence we do not necessarily expect all vortex streamlines to form closed surfaces correct to $O(\delta )$. We proceed by following the peak vorticity approach (Flierl & Haines Reference Flierl and Haines1994; Johnson & Crowe Reference Johnson and Crowe2021) and assuming that the innermost streamsurface of the streamfunction in a frame moving with the vortex

(6.7)\begin{equation} \psi' ={-}\frac{U}{2}\rho^2 + \psi, \end{equation}

remains closed to $O(\delta )$. The vorticity

(6.8)\begin{equation} \omega_\phi={-}\frac{1}{\rho}D^2 \psi, \end{equation}

contained within this streamsurface is assumed to be conserved. Here this innermost streamsurface corresponds to a torus of vanishing minor axis about the maximum of $\psi '$ so all quantities within this streamsurface can be evaluated at the position of this maximum, referred to from now on as the ‘vortex centre’.

To leading order we have

(6.9)\begin{equation} \omega_\phi = \frac{15U\rho}{2a^2}, \end{equation}

and as the maximum of $\psi '$ occurs for $r = a/\sqrt {2}$, $\theta = {\rm \pi}/2$ we have that the vorticity within the innermost streamsurface is

(6.10)\begin{equation} \omega_{\phi c} = \frac{15U}{2\sqrt{2}\,a}. \end{equation}

Taking this quantity to be constant to $O(\delta )$ gives that

(6.11)\begin{equation} \frac{\textrm{d}}{\textrm{d}T}\left[\frac{U}{a}\right] = 0, \end{equation}

hence $U/a$ is equal to its initial value, $U(T_0)/a(T_0)$, at $T = T_0$. We note that this approach may also be thought of as assuming that the value of $\omega _\phi$ is conserved following a particle within the innermost streamsurface.

7. Vortex decay

By combining the loss of vortex energy to the wave field, (6.6), with the conservation of centre vorticity, (6.11), we have

(7.1)\begin{equation} a(T) = a_0 \exp\left[ -\frac{7 a_0^3}{200 U_0^3} (T-T_0)\right], \end{equation}

and

(7.2)\begin{equation} U(T) = U_0 \exp\left[ -\frac{7 a_0^3}{200 U_0^3}(T-T_0)\right], \end{equation}

where $(U,a) = (U_0,a_0)$ at some initial time $T=T_0$ corresponding to $t = t_0 = \epsilon ^{-4} T_0$.

To leading order in $\epsilon$, the maximum of $\psi '$ within the vortex occurs for $r = a/\sqrt {2}$ and $\theta = {\rm \pi}/2$ and is given by

(7.3)\begin{equation} \psi_c = \frac{3Ua^2}{16}, \end{equation}

where we use a subscript ‘$c$’ to denote the value of a quantity at the vortex centre. Using (7.1) and (7.2) we can write

(7.4)\begin{equation} \psi_c = \frac{3U_0a_0^2}{16} \exp\left[ -\frac{21 a_0^3}{200 U_0^3} (T-T_0)\right], \end{equation}

which may be tested against numerical simulations. Additionally, the centre vorticity may be written using (6.10) as

(7.5)\begin{equation} \omega_{\phi c} = \frac{15 U_0}{2\sqrt{2} a_0}, \end{equation}

which clearly remains independent of time by the assumption of conservation of centre vorticity. Finally the leading order (in $\epsilon$) vortex energy should decay as

(7.6)\begin{equation} E_0 = \frac{10{\rm \pi} U_0^2 a_0^3}{7}\exp\left[ -\frac{7 a_0^3}{40 U_0^3}(T-T_0)\right]. \end{equation}

8. Numerical simulations

To test the model presented here we run numerical simulations of the non-dimensional, axisymmetric governing equations (see (2.2)) using the Dedalus package (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). Simulations are run for rotation parameters $\epsilon \in \{0,0.1,0.2,0.3,0.4,0.5\}$. We solve the axisymmetric Navier–Stokes equations on the domain $(\rho ,z) \in [0,51.2]\times [-51.2,51.2]$ with 1024 gridpoints in the $\rho$ direction and 2048 gridpoints in the $z$ direction. Fields are expanded in terms of Fourier modes in $z$ and Chebyshev polynomials in $\rho$ and integrated for $t \in [0,50]$ using a second-order semi-implicit backwards differentiation formula scheme with a timestep of $2\times 10^{-3}$. Due to the choice of basis functions, our fields satisfy periodic boundary conditions in the $z$ direction and we impose rigid wall conditions ($u_\rho = 0$) on $\rho = 51.2$. The boundary conditions at $\rho = 0$ are set by multiplying all equations through by the highest power of $\rho$ and substituting $\rho = 0$. This is dealt with internally by Dedalus and efficiently deals with the coordinate singularity.

Since we do not know $U(t)$ in advance, numerical simulations are performed in coordinates moving at speed $1$ in the $z$ direction by setting $U = 1$ in (2.2a). Therefore, a vortex moving with speed $U \approx 1$ should remain near the centre of the domain and a vortex found to be moving in the negative $z$ direction has $U < 1$. Viscous diffusion with a viscosity of $\nu = 10^{-4}$ is included for numerical stability. This diffusion acts to damp instabilities which would otherwise form near the rear stagnation point (Moffatt & Moore Reference Moffatt and Moore1978; Pozrikidis Reference Pozrikidis1986; Protas & Elcrat Reference Protas and Elcrat2016) leading to a breakdown of the vortex.

The $\rho$ and $z$ velocity components, $(u_\rho ,u_z)$, are initialised using the classical Hill's vortex solution (corresponding to the $\epsilon = 0$ limit of $\psi _0$ in (5.1)) with $U = 1$ and $a = 1$ while the swirl velocity $w$ is set by

(8.1)\begin{equation} w ={-}\frac{\epsilon k}{2\rho}\psi_0. \end{equation}

At early times we do not expect to see the vortex decay predicted in § 7 since this requires an established wave field. Instead we expect fast transient behaviour where the wave field is formed and the vortex adjusts to the effects of rotation. Once the transient adjustment near the vortex has abated, comparison with our theory can be made. We therefore begin all comparison from $t_0 = 10$ which we observe to be sufficient for the wave field to develop and the vortex to adjust. It should be noted that the vortex speed and radius may change due this transient adjustment phase and hence the values of $U(t_0) = U_0$ and $a(t_0) = a_0$ may differ from the values of $U = a = 1$ imposed at $t = 0$.

Figure 3 shows the streamfunction $\psi$ as a function of $(\rho ,z)$ for $\epsilon = 0.4$ and $t \in \{0,25,50\}$. The formation of the wave field can be clearly seen with waves forming behind the vortex of wavelength $\lambda = 2{\rm \pi} U/\epsilon \approx 16$. Small additional disturbances are seen in the wake of the vortex resulting from the initial transient evolution and from fluid escaping the vortex core.

Figure 3. The streamfunction, $\psi$, from numerical simulations with $\epsilon = 0.4$. Results are shown for (a) $t = 0$, (b) $t = 25$ and (c) $t = 50$.

A supplementary movie available at https://doi.org/10.1017/jfm.2021.386 shows the temporal evolution of the velocity fields for $\epsilon = 0.4$ – we opt to plot the velocity fields since they show the wave propagation more clearly than the streamfunction. Due to the fast decay of $u_\rho$ and $u_z$ with $\rho$, we plot $\rho u_\rho$ and $\rho u_z$ instead so that the wave field is visible for large $\rho$. We observe a wave field of similar structure to the prediction shown in figure 2(b) with no disturbances propagating upstream. At later times the vortex appears to move downwards, corresponding to a velocity, $U$, of less than $1$. Additionally, a small disturbance left by the initial adjustment of the vortex can be seen travelling downwards. Due to the periodicity of the domain in the $z$ direction and the boundary at $\rho = 51.2$ acting to reflect outgoing waves, we might expect waves to loop around the domain and interact with the vortex. While some disturbances are seen re-entering the domain, the simulation stop time, $t = 50$, is seen to be sufficient to prevent these waves from reaching the vortex.

To leading order in $\delta$ we expect that the swirl velocity, $w$, both inside and outside the vortex, satisfies

(8.2)\begin{equation} \rho w = K(\psi) ={-}\epsilon k \psi, \end{equation}

from (2.19) and (2.22). Similarly the vorticity and streamfunction are related by

(8.3)\begin{equation} -\rho \omega_\phi + \epsilon^2 k^2 \psi = \begin{cases} 0, & r > a,\\ -\dfrac{15U}{2a^2}\left[1-\dfrac{3}{7}(\epsilon k a)^2\right] \rho^2, & r < a, \end{cases} \end{equation}

using (2.30), (3.9) and $-\rho \omega _\phi = D^2\psi$. If these relations are satisfied throughout the vortex evolution then the vortex will remain close to a Hill's vortex and our theory should be valid. These relations tested are for $\epsilon = 0.4$ and $t = 30$ in figure 4 using scatter plots. These plots show the numerical values of two quantities at some point in $(\rho ,z)$ space; each point represents a different numerical gridpoint and all gridpoints are included. Points are coloured – red points represent points within the vortex ($r < 0.9a$), green points represent points near the vortex edge ($0.9a\leq r \leq 1.1a$) and blue points represent points outside the vortex ($r>1.1a$). Theoretical predictions using (8.2) and (8.3) are included as dashed lines.

Figure 4. Scatter plots for $\epsilon = 0.4$ and $t = 30$ showing (a) $\rho w$ against $\psi$, (b) $-\rho \omega _\phi + \epsilon ^2k^2\psi$ against $\rho ^2$ and (c) $\rho \omega _\phi$ against $\psi$. The dashed lines give our theoretical predictions.

Figure 4(a) shows the angular momentum, $\rho w$, against the streamfunction, $\psi$. We can see that all points lie very close to a line of constant gradient so (8.2) remains accurate late in the evolution. The value of $-\epsilon k$ can be determined by fitting the points to a straight line. Figure 4(b) shows $-\rho \omega _\phi + \epsilon ^2k^2\psi$ against $\rho ^2$ with a focus on the region close to the vortex, $\rho \leq 1$. Points close to the centre of the vortex lie close to our prediction (see (8.3)), however, we note that the predicted relation does not hold near the edge of the vortex where viscosity acts to smooth the discontinuity in vorticity. Our prediction is determined using the fitted value of $\epsilon k$ to determine $U = 1/k$ and the method described below to estimate $a$. Finally, figure 4(c) shows $\rho \omega _\phi$ against $\psi$ which, by (8.3), we expect to be linearly related outside the vortex. As expected we observe that the blue points (corresponding to points outside the vortex) lie close to a line of gradient $\epsilon ^2 k^2$.

While viscosity acts to maintain the vortex integrity by damping small scale instabilities, it also leads to a loss of around 10 % of the domain integrated energy. These viscous effects are particularly pronounced near the vortex edge as seen in figure 4 since Hill's vortex has a vorticity discontinuity across $r = a$. The smoothing of this discontinuity makes it particularly difficult to determine an accurate value for the vortex radius, $a$, while the drag associated with the boundary layers can modify the vortex speed, $U$.

To account for viscous effects, we assume that a diagnosed quantity, $f$, is a function of the linear wave drag with coefficient $\epsilon ^4$ and the viscous damping with coefficient $\nu$. If both $\epsilon ^4$ and $\nu$ are small, we may expand

(8.4)\begin{equation} f(\epsilon^4,\nu,t) = f(0,0,t)+\epsilon^4 \frac{\partial f}{\partial \epsilon^4}(0,0,t)+\nu \frac{\partial f}{\partial \nu}(0,0,t)+O(\epsilon^8,\epsilon^4\nu,\nu^2), \end{equation}

and hence we may determine the effect of the wave drag by running a simulation with $\epsilon = 0$ and subtracting the results from a non-zero $\epsilon$ run to get

(8.5)\begin{align} f(\epsilon^4,\nu,t) - [f(0,\nu,t) -f(0,0,t)] = f(0,0,t)+\epsilon^4\frac{\partial f}{\partial \epsilon^4}(0,0,t) +O(\epsilon^8,\epsilon^4\nu,\nu^2), \end{align}

where $f(0,0,t) = f(0,0,0)$ is constant in time and known from the analytical Hill's vortex solution or equivalently the initial conditions. We now use this method to estimate the maximum streamfunction, $\psi _c$, and vortex energy, $E_0$, adjusting for viscosity. The unadjusted value of $E_0$ is determined by integrating $\psi \,\omega _\phi$ over the region $r< a$ (see (A10)) where the value of $a$ is estimated as the radial distance from the centre of the vortex to the point where the vorticity is reduced to 10 % of its maximum value.

Figure 5 shows the centre streamfunction, $\psi _c$, centre vorticity, $\omega _{\phi c}$, and vortex energy, $E_0$. Here, the centre streamfunction, $\psi _c$, is the maximum value of $\psi '$ (defined in (6.7)) inside the vortex and the centre vorticity, $\omega _{\phi c}$, is the value of $\omega _\phi$ at the position where $\psi ' = \psi _c$. As stated above, we perform all comparisons from $t = t_0 = 10$ onwards; we also normalise all quantities by their value at $t_0$. Solid lines denote data from our numerical simulations and dashed lines denote our predictions from § 7. Due to the difficulty in determining the radius and speed of the adjusted vortex to the accuracy required to calculate $(a_0/U_0)^3$ at $t = t_0$, we instead take $a_0/U_0 = 1.15$ by fitting the predictions for $\psi _c$ to the numerical data in $t \in [10 ,20]$. We expect the value of $a_0/U_0$ to depend on $\epsilon$ so perform a fitting for all value of $\epsilon$ and find that all values satisfy $a_0/U_0\approx 1.15$. It should be noted that this value of $a_0/U_0$ lies within the range of possible values calculated using different methods for determining the vortex speed and radius so is consistent with our numerical data.

Figure 5. (a) The centre streamfunction, $\psi _c$, (b) the centre vorticity, $\omega _{\phi c}$, (c) the vortex energy, $E_0$. All quantities are normalised by their value at $t = 10$ and shown from $t = 10$ for all non-zero $\epsilon$ values simulated.

Figure 5(a) shows the comparison between the numerical and theoretical values of the centre streamfunction, $\psi _c$. We observe close agreement between the prediction and the theory though, as noted previously, some fitting was involved to account for the difficulties in determining $a_0/U_0$. As the results are consistent over the full range of $\epsilon$, we conclude that the time scaling $t\sim O(1/\epsilon ^4)$ is likely correct. Figure 5(b) shows the centre vorticity, $\omega _{\phi c}$, from our numerical simulations. We observe that $\omega _{\phi c}$ is conserved to within 3 % throughout the evolution for $\epsilon \leq 0.4$ which is consistent with our assumption that the centre vorticity is conserved throughout the evolution. We do observe larger, oscillatory deviations for $\epsilon = 0.5$ which may suggest the presence of wavelike oscillations inside the vortex due to the stronger rotation – these phenomena are likely to be strongly nonlinear and may act to modify the vortex decay.

Figure 5(c) shows a comparison between the vortex energy, $E_0$, from simulations and theory. We note that our predictions are less accurate for the larger $\epsilon$ values and this discrepancy may be due to several factors. Firstly the difficulties in determining the vortex radius are particularly important here due to the integration domain depending on $a$. Secondly, the expression for $E_0$ does not include the rotational component of the energy, $w^2$, as this is an $O(\epsilon ^2)$ contribution. However, our numerical simulations suggest that this component may be significant in the finite $\epsilon$ case and may even increase with time due to finite $\epsilon$ effects which could be related to conservation of angular momentum. Including this component would increase the predicted value back towards the numerical value. Finally, the neglected higher-order components of the translational kinetic energy $T = |\boldsymbol {\nabla }\psi |^2/(2\rho ^2)$, and the wave drag will modify the energy balance in the finite $\epsilon$ case.