2.1. Equation of motion

Let ${\boldsymbol{x}}_c(t)$ denote a critical point of a smooth, two-dimensional unsteady vorticity field $\omega ({\boldsymbol{x}},t)$ . By definition, it fulfils

(2.1) \begin{equation} \boldsymbol{\nabla }\omega ({\boldsymbol{x}}_c(t),t) = \boldsymbol{0}. \end{equation}

If the Hessian $\unicode{x1D643} = \boldsymbol{\nabla }(\boldsymbol{\nabla }\omega )$ in rectangular coordinates $(x,y)$ given by

(2.2) \begin{equation} \unicode{x1D643} = \begin{pmatrix} \partial _{xx}\omega & \partial _{xy}\omega \\ \partial _{xy}\omega & \partial _{yy}\omega \end{pmatrix}, \end{equation}

evaluated at ${\boldsymbol{x}}_c(t)$ , is regular, then the type of the critical point is determined by the sign of the determinant, $\det \unicode{x1D643}$ . If $\det \unicode{x1D643}\gt 0$ , then ${\boldsymbol{x}}_c$ is a local extremum of vorticity; if $\det \unicode{x1D643}\lt 0$ , then it is a saddle point. Local maxima of positive vorticity and local minima of negative vorticity are considered centres of vortices. The physical interpretation of saddle points is less obvious, but as we will see, they play an important topological role in the creation and destruction of vortices.

The eigenvalues $\lambda _1, \lambda _2$ of $\unicode{x1D643}$ describe the shape of a vortex in a neighbourhood of an extremum. Close to such a point, the iso-vorticity curves are closed, and to second order in the distance, the curves are ellipses. If the half-axes are denoted $a,b$ , then the eccentricity is

(2.3) \begin{equation} e = \sqrt {1 - b^2/a^2} = \sqrt {1-\lambda _1/\lambda _2}, \end{equation}

where $|\lambda _1|\leqslant |\lambda _2|$ is assumed. The eccentricity is a measure of how much the shape of the vortex deviates from axisymmetry. The absolute values of the eigenvalues measure the concentration of the vortex. The larger they are, the more rapidly the magnitude of vorticity decays away from the extremum.

Implicit differentiation of (2.1) yields

(2.4) \begin{equation} \unicode{x1D643}\dot {\boldsymbol{x}}_c + \partial _t\boldsymbol{\nabla }\omega = \boldsymbol{0}, \end{equation}

from which the velocity of the critical point can be found:

(2.5) \begin{equation} \dot {\boldsymbol{x}}_c = -\unicode{x1D643}^{-1}\,\partial _t\boldsymbol{\nabla }\omega . \end{equation}

Further information can be obtained from the vorticity transport equation

(2.6) \begin{equation} \partial _t\omega = -{\boldsymbol{u}}{\boldsymbol \cdot} \boldsymbol{\nabla }\omega + \nu \Lambda, \end{equation}

where $\boldsymbol{u}$ is the underlying fluid velocity field, $\nu$ is the kinematic viscosity, and $\Lambda =\Delta \omega$ is the Laplacian of the vorticity. Taking the gradient on both sides of this equation, one obtains

(2.7) \begin{equation} \partial _t\boldsymbol{\nabla }\omega = -(\boldsymbol{\nabla }{\boldsymbol{u}})\,\boldsymbol{\nabla }\omega - \unicode{x1D643}{\boldsymbol{u}} + \nu\, \boldsymbol{\nabla }\Lambda . \end{equation}

Evaluating at a critical point of vorticity where $\boldsymbol{\nabla }\omega = \boldsymbol{0}$ , this reduces to

(2.8) \begin{equation} \partial _t\boldsymbol{\nabla }\omega = - \unicode{x1D643}{\boldsymbol{u}} + \nu\, \boldsymbol{\nabla }\Lambda . \end{equation}

Inserting (2.8) in (2.5) then gives

(2.9) \begin{equation} \dot {\boldsymbol{x}}_c = {\boldsymbol{u}} - \nu \unicode{x1D643}^{-1}\,\boldsymbol{\nabla }\Lambda . \end{equation}

This equation of motion for a critical point of vorticity was first derived in rectangular coordinates by Brøns & Bisgaard (Reference Brøns and Bisgaard2010), but not analysed in any detail. The form (2.9) is a general, coordinate-free vector equation. It expresses the velocity of the critical point $\dot {\boldsymbol{x}}_c$ as the sum of the advection by the fluid ${\boldsymbol{v}}_a={\boldsymbol{u}}$ and a viscous drift velocity ${\boldsymbol{v}}_d=-\nu \unicode{x1D643}^{-1}\,\boldsymbol{\nabla }\Lambda$ . The viscous drift velocity is proportional to viscosity but otherwise depends only on spatial derivatives of the vorticity field, evaluated at the critical point. However, the effect of viscosity is not isolated to the viscous drift velocity. The velocity field itself satisfies the Navier–Stokes equations, hence the advection velocity ${\boldsymbol{v}}_a$ implicitly depends on the viscosity too.

2.2. The viscous drift velocity

The viscous drift velocity obviously vanishes if the viscosity is zero. Thus in an inviscid fluid, the critical points of vorticity are material points. This also follows directly from the fact that in an inviscid fluid, vorticity is simply advected by the flow.

The viscous drift velocity is also zero if

(2.10) \begin{equation} \boldsymbol{\nabla }\Lambda =\boldsymbol{0} \end{equation}

at the critical point. To understand (2.10), consider a Cartesian coordinate system $(x,y)$ centred at a critical point of vorticity, with corresponding polar coordinates $(r,\theta )$ . In this coordinate system, (2.10) becomes

(2.11) \begin{equation} \partial _{xxx}\omega _0 + \partial _{xyy}\omega _0 = 0, \quad \partial _{xxy}\omega _0 + \partial _{yyy}\omega _0 = 0, \end{equation}

where the subscript $0$ denotes evaluation at the critical point.

We now consider a Fourier expansion of the vorticity, $\omega = \sum _{n=-\infty }^{\infty }c_n(r)\,\textrm{e}^{\textrm{i} n\theta }$ . A Taylor expansion in $r$ then shows that the third-order derivatives of $\omega$ occur only in the terms

(2.12) \begin{align} &\,\,\,c_1 = \overline {c_{-1}}= \frac {1}{16}(\partial _{xxx}\omega _0 + \partial _{xyy}\omega _0 - \textrm{i}(\partial _{xxy}\omega _0 + \partial _{yyy}\omega _0))r^3 + \mathcal O(r^5), \end{align}

(2.13) \begin{align} &c_3 = \overline { c_{-3}}= \frac {1}{48}(\partial _{xxx}\omega _0 - 3\partial _{xyy}\omega _0 -\textrm{i}(3\partial _{xxy}\omega _0 - \partial _{yyy}\omega _0))r^3 + \mathcal O(r^5). \end{align}

Thus from (2.11) and (2.12), the critical point follows the flow if and only if the basic angular mode of vorticity, represented by $c_1$ , is zero to third order in $r$ . In that case, (2.13) becomes

(2.14) \begin{equation} c_3 =\overline { c_{-3}} = \frac {1}{12}(\partial _{xxx}\omega _0 + \textrm{i}\,\partial _{yyy}\omega _0)r^3 + \mathcal O(r^5), \end{equation}

which in general is non-zero.

We note that (2.10) is satisfied for an axisymmetric vorticity field since all $c_n$ with $n\neq 0$ then vanish. Interestingly, there is no connection to the eccentricity (2.3); vortices of any local shape as defined by the eccentricity may (or may not) follow the flow. We stress that (2.10) will, in general, hold at only one time instant. The acceleration of a critical point of vorticity will usually be different from the acceleration of the fluid particle on which it resides at the given time instant. In that case, (2.10) will be violated immediately. Only if (2.10) holds for an interval of time will a critical point of vorticity also be a material point.

If the eigenvalues of $\unicode{x1D643}$ are numerically large, the eigenvalues of $\unicode{x1D643}^{-1}$ are small, and the viscous velocity correction will also be small. That is, the more concentrated a vortex is, the closer it follows the fluid flow – for given viscosity and $\boldsymbol{\nabla }\Lambda$ .

The strength of a vortex measured by the value $\omega _0$ at the extremum plays no role for its motion, as (2.9) depends only on derivatives of $\omega$ . In particular, the analysis above also holds for local maxima of negative vorticity and local minima of positive vorticity. We call such points anti-vortices. They are reminiscent of the holes inside a $Q$ -vortex observed by Nielsen et al. (Reference Nielsen, Andersen, Hansen and Brøns2021). The rate of change of vorticity at a critical point is

(2.15) \begin{equation} \frac {\textrm{d}}{\textrm{d} t}\omega ({\boldsymbol{x}}_c(t),t) = \boldsymbol{\nabla }\omega {\boldsymbol \cdot} \dot {\boldsymbol{x}}_c + \partial _t\omega = \boldsymbol{\nabla }\omega {\boldsymbol \cdot} (\dot {\boldsymbol{x}}_c-{\boldsymbol{u}}) + \nu \Lambda = \nu \Lambda, \end{equation}

where we used the vorticity transport equation (2.6) and the definition (2.1) of a critical point. Since $\Lambda = {\textrm{tr}}(\unicode{x1D643})$ , it follows that vorticity decreases for a maximum of vorticity and increases for a minimum of vorticity. Hence an anti-vortex will preserve its type, whereas a vortex can turn into an anti-vortex as time progresses. For saddle points, $\Lambda$ can have either sign.

3.1. The generic case

As long as $\unicode{x1D643}$ is regular at a critical point of vorticity, the point will move according to (2.9), and the type of the point (extremum or saddle) will not change. If, however, at some time instant, $\unicode{x1D643}$ becomes singular, then a local bifurcation event is expected, as the velocity of the critical point is no longer well-defined from (2.4). Here, we consider the case where $\unicode{x1D643}$ attains a simple zero eigenvalue at $t=0$ . We choose the coordinate system such that the critical point is at the origin, and $\unicode{x1D643}$ is diagonal. The Hessian then has a simple zero eigenvalue exactly when one of $\partial _{xx}\omega$ and $\partial _{yy}\omega$ is zero. Here, we assume

(3.1) \begin{equation} \partial _{xy}\omega _0 = 0, \quad \partial _{xx}\omega _0 = 0, \quad \partial _{yy}\omega _0 \neq 0, \end{equation}

where the subscript $0$ here and in the following indicates evaluation at $(x,y,t)=(0,0,0)$ , so any quantities with a subscript $0$ are constants. Specifically, we have

(3.2) \begin{equation} \unicode{x1D643}_0 = \begin{pmatrix} 0 & 0 \\ 0 & \partial _{yy}\omega _0 \end{pmatrix}, \end{equation}

and since $\unicode{x1D643}_0$ is singular, it is not possible to solve $\boldsymbol{\nabla }\omega (x,y,t) = \boldsymbol{0}$ for $x,y$ with the implicit function theorem as in § 2.1. However, assuming for now the non-degeneracy condition

(3.3) \begin{equation} \nu\, \partial _x\Lambda _0\neq 0, \end{equation}

the equation can be solved for $y,t$ since the Jacobian matrix

(3.4) \begin{equation} \unicode{x1D645}_0 = \left . \frac {\partial (\partial _x\omega, \partial _y\omega )}{\partial (y,t)}\right |_0 = \begin{pmatrix} \partial _{xy}\omega _0 & \partial _{xt}\omega _0 \\ \partial _{yy}\omega _0 & \partial _{yt}\omega _{0} \end{pmatrix} = \begin{pmatrix} 0 & \nu\, \partial _x\Lambda _0 \\ \partial _{yy}\omega _0 & -\partial _{yy}\omega _0v_0 + \nu\, \partial _{y}\Lambda _0 \end{pmatrix} \end{equation}

is then regular. In § 3.2, we consider a symmetric flow where (3.3) is violated. The expressions in the last column of the last matrix follow from (2.8), using that $\unicode{x1D643}_0{\boldsymbol{u}}_0 = (0,\partial _{yy}\omega _0v_0)^{\textrm{T}}$ , where ${\boldsymbol{u}}_0=(u_0,v_0)^{\textrm{T}}$ is the velocity at the origin. Thus it follows from the implicit function theorem that there exist functions $y=y^{*}(x),\ t=t^{*}(x)$ defined for $x$ close to $0$ with $y^{*}(0)=0,\ t^{*}(0)=0$ such that

(3.5) \begin{equation} \partial _x\omega (x,y^{*}(x),t^{*}(x)) = 0, \quad \partial _y\omega (x,y^{*}(x),t^{*}(x)) = 0. \end{equation}

Derivatives of $y^{*}, t^{*}$ at $x=0$ can be found by repeated differentiation of (3.5), the first two orders being

(3.6) \begin{equation} \begin{pmatrix} y^{*\prime }(0) \\ t^{*\prime }(0) \end{pmatrix} = -\unicode{x1D645}_0^{-1} \begin{pmatrix} \partial _{xx}\omega _0 \\ \partial _{xy}\omega _0 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \end{equation}

and

(3.7) \begin{equation} \begin{pmatrix} y^{*\prime \prime }(0) \\ t^{*\prime \prime }(0) \end{pmatrix} = -\unicode{x1D645}_0^{-1} \begin{pmatrix} \partial _{xxx}\omega _0 \\ \partial _{xxy}\omega _0 \end{pmatrix} = \begin{pmatrix} -\dfrac {\partial _{xxx}\omega _0v_0}{\nu\, \partial _x\Lambda _0} + \dfrac {\partial _{xxx}\omega _0\, \partial _{y}\Lambda _0 - \partial _{xxy}\omega _0\,\partial _{x}\Lambda _0}{\partial _{yy}\omega _0\,\partial _x\Lambda _0}\\[10pt] -\dfrac {\partial _{xxx}\omega _0}{\nu\, \partial _x\Lambda _0} \end{pmatrix}, \end{equation}

where (3.1) has been used. A Taylor expansion yields

(3.8) \begin{equation} y^{*}(x) = \frac {y^{*\prime \prime }(0)}{2}x^2 + \mathcal O(x^3),\quad t^{*}(x) = \frac {t^{*\prime \prime }(0)}{2}x^2 + \mathcal O(x^3), \end{equation}

and further assuming the non-degeneracy condition

(3.9) \begin{equation} \partial _{xxx}\omega _0 \neq 0, \end{equation}

it follows that $t^{*\prime \prime }(0) \neq 0$ such that (3.8) can be solved to yield trajectories of critical points:

(3.10) \begin{equation} x_c(t) = \pm \sqrt {\frac {2t}{t^{*\prime \prime }(0)}} + \mathcal O(t), \quad y_c(t) = \frac {y^{*\prime \prime }(0)}{t^{*\prime \prime }(0)}t + \mathcal O(t^{3/2}). \end{equation}

The dominating terms can be found from (3.7). For $x_c(t)$ , it is interesting to resolve the $\mathcal O(t)$ term, which requires finding $t^{*\prime \prime \prime }(0)$ . We omit the tedious computations and only state the final result:

(3.11a) \begin{align} &x_c(t) = \pm \sqrt {-2\nu \frac {\partial _x\Lambda _0}{\partial _{xxx}\omega _0}t} \nonumber \\ &{}+ \left (u_0 + \nu \left [\frac {\partial _{xxx}\omega _0\,\partial _{xxy}\omega _0\,\partial _y\Lambda _0 - \partial _{xxy}\omega _0^2\,\partial _x\Lambda _0}{\partial _{xxx}\omega _0^2\,\partial _{yy}\omega _0} + \frac {\partial _{xxxx}\omega _0\,\partial _x\Lambda _0 - 3\partial _{xxx}\omega _0\,\partial _{xx}\Lambda _0}{3\partial _{xxx}\omega _0^2}\right]\right)t\nonumber \\ &{}+ \mathcal O(t^{3/2}), \end{align}

(3.11b) \begin{align} &\qquad\qquad y_c(t) = \left ( v_0 + \nu \frac {\partial _{xxy}\omega _0\,\partial _x\Lambda _0 - \partial _{xxx}\omega _0\,\partial _y\Lambda _0}{\partial _{yy}\omega _0\,\partial _{xxx}\omega _0}\right )t + \mathcal O(t^{3/2}). \end{align}

If $\partial _x\Lambda _0/\partial _{xxx}\omega _0\gt 0$ , then two critical points exist for $t\lt 0$ , one for each sign of $x_c(t)$ ; these merge and disappear at the origin at $t=0$ . If $\partial _x\Lambda _0/\partial _{xxx}\omega _0\lt 0$ , then two critical points are created at $t=0$ and exist for $t\gt 0$ . The types of the critical points are found from the eigenvalues of the Hessian

(3.12) \begin{equation} \unicode{x1D643}(x,y^{*}(x),t^{*}(x)) = \begin{pmatrix} \partial _{xxx}\omega _0x + \mathcal O(x^2) & \mathcal O(x)\\ \mathcal O(x) & \partial _{yy}\omega _0 + \mathcal O(x) \end{pmatrix}. \end{equation}

The eigenvalues are $\lambda _1 = \partial _{xxx}\omega _0x + {O}(x^2)$ and $\lambda _2 = \partial _{yy}\omega _0 + {O}(x)$ . Since $x\gt 0$ for one critical point and $x\lt 0$ for the other, one point is an extremum, and the other is a saddle.

The bifurcation is therefore a cusp or a saddle–centre bifurcation, and it is illustrated in figure 1. The sequence (a) $\rightarrow$ (b) $\rightarrow$ (c) shows the merging of an extremum and a saddle. The critical points follow each of the two branches of a locally parabolic trajectory for $t\lt 0$ (figure 1 a), coalesce into a degenerate critical point at $t=0$ (figure 1 b), and then disappear for $t\gt 0$ (figure 1 c).

The bifurcation occurs in viscous flows only, as the second condition in (3.3) requires $\nu \neq 0$ . In inviscid flow, critical points of vorticity never merge, and no new ones are created.

Figure 1.

Illustration of annihilation of critical points in a cusp bifurcation. The grey curves are iso-lines of vorticity. The red (black) point is the extremum (saddle point) of vorticity, and the red and black curves are the corresponding trajectories: (a) before bifurcation; (b) at the bifurcation point; (c) after bifurcation. For critical point creation, the temporal order is reversed.

The terms linear in $t$ in (3.11) have the same structure as in the regular case (2.9); the velocity of a critical point is the sum of an advection with the fluid and a drift velocity proportional to viscosity. For the motion in the $y$ -direction, the linear term is dominating near the bifurcation, but in the $x$ -direction, the leading-order term of the position is of order $t^{1/2}$ , which leads to a velocity that goes to infinity as $|t|^{-1/2}$ when the bifurcation is approached. Hence the viscous drift dominates the motion of the critical points completely, no matter the fluid velocity. That viscous drift must be present near the bifurcation is already clear from (3.3). It implies that $\boldsymbol{\nabla }\Lambda _0 \neq \boldsymbol{0}$ , hence violating (2.10) for zero viscous drift.

The eccentricity (2.3) of the level curves of vorticity close to the extremum of vorticity when the bifurcation is approached is

(3.13) \begin{equation} e = \sqrt {1 - \frac {\partial _{xxx}\omega _0}{\partial _{yy}\omega _0}x + \mathcal O(x^2)} \approx 1 - \frac {\partial _{xxx}\omega _0}{2\partial _{yy}\omega _0}x \approx 1 - \sqrt {\frac {\nu }{2}\frac {|\partial _{xxx}\omega _0\,\partial _x\Lambda _0|}{\left (\partial _{yy}\omega _0\right )^2}}\sqrt {|t|}, \end{equation}

where we used (3.11a ) in the last step. This is consistent with the plot of the vortex annihilation shown in figure 1. The almost elliptic iso-vorticity curves close to the extremum are squeezed together in the $y$ -direction faster than they shrink in the $x$ -direction, as measured by the eccentricity, which tends to the largest possible value $1$ as $t \rightarrow 0$ .

3.2. Symmetric bifurcation

If two vortices with identical distribution of vorticity are placed symmetrically around the origin in a Cartesian coordinate system, then the symmetry

(3.14) \begin{equation} \omega (-x,-y) = \omega (x,y) \end{equation}

will be fulfilled for the total vorticity of the two distributions. As the flow evolves, the symmetry will be preserved such that (3.14) holds for all times. This implies that all odd-order spatial derivatives at the origin are zero, and it follows that the origin is always a critical point of vorticity. If the Hessian at the origin is singular at some time instant, then a bifurcation event is expected. However, the analysis of § 3.1 does not apply; the non-degeneracy condition (3.3) is violated since $\partial _x\Lambda _0$ is a sum of third-order derivatives of $\omega$ . Here, we therefore derive the bifurcation structure under the symmetry (3.14), which we use to study the merging of two identical vortices in § 4. We omit some details as the analysis to a large extent follows that of § 3.1.

We consider a coordinate system such that the Hessian at the origin is singular and diagonal at $t=0$ under the assumptions (3.1), yielding the Hessian

(3.15) \begin{equation} \unicode{x1D643}(0,0,t) = \begin{pmatrix} \partial _{xxt}\omega _0t + \mathcal O(t^2) & \partial _{xyt}\omega _0t + \mathcal O(t^2) \\ \partial _{xyt}\omega _0t + \mathcal O(t^2) & \partial _{yy}\omega _0 + \mathcal O(t) \end{pmatrix}. \end{equation}

Furthermore, we assume the non-degeneracy condition

(3.16) \begin{equation} \partial _{xxt}\omega _0 = \nu\, \partial _{xx}\Lambda _0 \neq 0, \end{equation}

where the identity is obtained by differentiation of (2.8). As before, the subscript indicates evaluation at the bifurcation point $(x,y,t)=(0,0,0)$ . It follows that

(3.17) \begin{equation} \det \unicode{x1D643}(0,0,t) = \nu\, \partial _{xx}\Lambda _0\,\partial _{yy}\omega _0 t + \mathcal O(t^2), \end{equation}

such that the critical point at the origin changes between an extremum and a saddle at $t=0$ .

We introduce a new variable $\eta$ by $y=\eta x$ , and define

(3.18) \begin{equation} \tilde {f}(x,\eta, t) = \partial _x\omega (x,\eta x, t), \quad \tilde {g}(x,\eta, t) = \partial _y\omega (x,\eta x, t), \end{equation}

and note that $\tilde {f}(0,\eta, t) = \tilde {g}(0,\eta, t) = 0$ for all $\eta, t$ . It follows that $x$ is a factor in both $\tilde {f}$ and $\tilde {g}$ , and that critical points of vorticity other than the origin can be found from

(3.19) \begin{equation} f(x,\eta, t) = \tilde {f}(x,\eta, t)/x = 0, \quad g(x,\eta, t) = \tilde {g}(x,\eta, t)/x = 0. \end{equation}

Consider the Jacobian

(3.20) \begin{equation} \unicode{x1D645}_0 = \left . \frac {\partial (f,g)}{\partial (\eta, t)}\right |_0 = \begin{pmatrix} 0 & \partial _{xxt}\omega _0 \\ \partial _{yy}\omega _0 & \partial _{xyt}\omega _{0} \end{pmatrix} = \begin{pmatrix} 0 & \nu\, \partial _{xx}\Lambda _{0} \\ \partial _{yy}\omega _0 & -\partial _{yy}\omega _0\,\partial _xv_0 + \nu\, \partial _{xy}\Lambda _0 \end{pmatrix}, \end{equation}

where $\partial _{xyt}\omega _0$ is obtained by differentiation of (2.8). Since $\unicode{x1D645}_0$ is regular, it follows from the implicit function theorem that there exist functions $\eta =\eta ^{*}(x),\ t=t^{*}(x)$ defined for $x$ close to 0, with $\eta ^{*}(0)=0,\ t^{*}(0)=0$ , such that

(3.21) \begin{equation} f(x,\eta ^{*}(x),t^{*}(x)) = 0, \quad g(x,\eta ^{*}(x),t^{*}(x)) = 0. \end{equation}

By implicit differentiation of (3.21), one finds

(3.22) \begin{equation} \begin{aligned} \eta ^{*\prime }(0) = 0,\quad \eta ^{*\prime \prime }(0) & = -\frac {\partial _{xxxx}\omega _0\,\partial _xv_0}{3\nu \partial _{xx}\Lambda _0} + \frac {\partial _{xxxx}\omega _0\,\partial _{xy}\Lambda _0 - \partial _{xxxy}\omega _0\,\partial _{xx}\Lambda _0}{3\partial _{yy}\omega _0\,\partial _{xx}\Lambda _0}, \\ t^{*\prime }(0) = 0,\quad t^{*\prime \prime }(0) & = -\frac {\partial _{xxxx}\omega _0}{3\nu \partial _{xx}\Lambda _0}. \end{aligned} \end{equation}

Adding the non-degeneracy condition

(3.23) \begin{equation} \partial _{xxxx}\omega _0 \neq 0, \end{equation}

the Taylor expansion $t=t^{*\prime \prime }(0)\,x^2/2 + \mathcal O(x^3)$ can be solved for $x$ , and inserting the solution into $y = \eta x = \eta ^{*\prime \prime }(0)\,x^3/2 + \mathcal O(x^4)$ yields trajectories of two critical points:

(3.24) \begin{equation} x_c(t) = \pm \sqrt {-6\nu \frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}t} + \mathcal O(t), \quad y_c(t) = \pm \frac {1}{2}\eta ^{*\prime \prime }(0)\left (-6\nu \frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}t\right )^{3/2} + \mathcal O(t^2). \end{equation}

The two non-degeneracy conditions (3.16) and (3.23) will be fulfilled generically if no further symmetries or similar constraints are imposed. In the case that $\partial _{yy}\omega _0$ , $\partial _{xx}\Lambda _0$ and $\partial _{xxxx}\omega _0$ all have the same sign, it follows from (3.17) that the origin changes from a saddle to an extremum at $t=0$ . A simple computation shows that the critical points (3.24) that exist for $t\lt 0$ are both extrema. Hence the bifurcation is a supercritical pitchfork that describes the merging of two vortices as illustrated in figure 2.

If $(r,\theta )$ are polar coordinates for one of the extrema, then (3.24) yields

(3.25) \begin{equation} r^2 = \left |6\nu \frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}t\right | + \mathcal O(t^{3/2}), \quad \theta = \left |3\eta ^{*\prime \prime }(0)\,\nu \frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}t\right | + \mathcal O(t^{3/2}). \end{equation}

In analogy with the generic case of § 3.1, it follows that the radial velocity $\dot r$ of the extrema scales as $|t|^{-1/2}$ . The angular velocity of the line connecting the two extrema scales as a constant as the bifurcation is approached.

Figure 2.

Illustration of the final phase of merging of identical vortices in a pitchfork bifurcation. The grey curves are iso-lines of vorticity. The red (black) points are the extrema (saddle points) of vorticity, and the red curve is the trajectory of the two extrema. (a) Before bifurcation: the heavy grey curve is the separatrix loop associated with the saddle point. (b) At the bifurcation point: the critical points coalesce at the origin. (c) After bifurcation: a single extremum remains at the origin.

4.1. Numerical method

We solved the vorticity transport equation (4.2) numerically using a forward time centred space algorithm (E & Liu Reference E and Liu1996a , Reference E and Liub ; Andersen et al. Reference Andersen, Schreck, Hansen and Brøns2019), keeping the time step well below the von Neumann stability criterion. The Poisson equation was solved using a sparse linear equation solver (Hansen Reference Hansen2011; Eaton et al. Reference Eaton, Bateman, Hauberg and Wehbring2021). We validated the accuracy of the scheme by repeating the simulation for the largest Reynolds number ( $Re=1500$ ) against a central difference scheme using a fourth-order Runge–Kutta integrator. Computations were performed on a square grid with $500\times500$ grid points, using domain size $L/d=5$ . When imposing the initial condition (4.3), we used $\sigma _0 = 0.015$ , to ensure that the initial vorticity field resembles that associated with two point vortices. Tests on a finer grid show that the initial condition is sufficiently resolved to ensure grid-independent results. The initial time $t_0$ was set using the approach of Josserand & Rossi (Reference Josserand and Rossi2007), who determined $t_0$ such that $t=0$ corresponds to the time when the isolated Gaussian vortices in (4.3) would have zero core radius. On the overall time scale of the vortex merging process, this is so small that the difference from choosing $t_0=0$ is negligible.

Given the results of the simulations, we determined the critical points of vorticity (where $\boldsymbol{\nabla } \omega = {\textbf {0}}$ ) as the intersections of the iso-contour lines $\partial _x \omega =0$ and $\partial _y \omega =0$ . All spatial derivative calculations were performed using a fourth-order centred finite difference operator. The vorticity, the Hessian and $\boldsymbol{\nabla }\Lambda$ were found at the critical points by a two-dimensional cubic interpolation scheme. The velocity of the critical points was predicted at each time step using these quantities and (2.9).

We captured the motion of the critical points by tracking their pathlines, and calculated the velocity of a critical point by using a fourth-order centred finite difference operator. Figure 3(a) shows the trajectories of the vortex centres for the case $Re = 1500$ . We obtained perfect agreement between the velocity calculated from the trajectories and that predicted by (2.9), thus validating both the equation of motion and the numerical discretisation schemes used.

4.2. Numerical results

Figure 3.

(a) Trajectories of vortex centres for $Re = 1500$ with initial positions $(1/2,0)$ (black) and $(-1/2,0)$ (dotted green). The initial positions are marked with filled circles. The open circle marks the origin, where the merging occurs. The motion is anticlockwise. (b) Instantaneous positions of vortex centres of identical vortices $1,2$ (red) with definitions of distance $d(t)$ , angle $\theta (t)$ and polar decomposition $(v_r,v_{\theta })$ of the velocity of vortex centre 1.

Figure 3(b) shows the quantities that we will consider in the following: the distance $d(t)$ between the vortex centres, the angle $\theta (t)$ between the line connecting the vortex centres and the $x$ -axis, and the velocity $\dot {\boldsymbol{x}}_c$ of a critical point decomposed into its radial and azimuthal components, $\dot {\boldsymbol{x}}_c = v_r{\boldsymbol{e}}_r + v_{\theta }{\boldsymbol{e}}_{\theta } = \dot r{\boldsymbol{e}}_r + r\dot \theta {\boldsymbol{e}}_{\theta }$ .

Figure 4.

Distance $d(t)$ between the vortex centres during vortex merging. Scaling from (4.4) close to the bifurcation is shown as dashed black lines. Merging times are shown as inverted triangles. The upright triangles mark the inflection points.

Figure 4 shows the distance $d(t)$ as a function of time $t$ for a range of Reynolds numbers $Re$ . A common feature for all $Re$ is that initially $d(t)$ is almost constant. This can be understood from a visual inspection of the vorticity field, which shows that the vortices are highly concentrated and almost axisymmetric. The axisymmetry implies that the anti-symmetric vorticity field in the co-rotating frame is vanishing, hence there is almost no radial advection. As discussed in § 2.2, the axisymmetry and the high concentration of vorticity both imply that the viscous drift is also small, and hence that the total radial velocity of the critical points is initially small. Another common feature is the final behaviour as merging is approached. Here, the motion is completely dominated by viscous drift according to the bifurcation analysis in § 3.2. The theory predicts, from (3.25), a scaling

(4.4) \begin{equation} d(t)^2 = -a_dt + b_d \end{equation}

close to merging, where

(4.5) \begin{equation} a_d = \frac {12}{Re}\left |\frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}\right |,\quad b_d = \frac {12}{Re}\left |\frac {\partial _{xx}\Lambda _0}{\partial _{xxxx}\omega _0}\right | t_m, \end{equation}

with $t_m$ denoting the time of merging (bifurcation). Note that the coordinate system used here is that of § 3.2, which is aligned with the eigenvectors of $\unicode{x1D643}_0 = \unicode{x1D643}(0,0,0)$ as shown in figure 2, and not the computational frame of figure 3; the two systems differ by a rotation around the origin.

We determine $t_m$ by interpolation as the time where the Hessian $\unicode{x1D643}(0,0,t)$ has a zero eigenvalue. The merging times are marked with inverted triangles in figure 4. Computing $a_d,b_d$ in the rotated coordinate system yields for (4.4) the dashed lines in figure 4. The very good agreement between the scaling and the simulation data confirms the theory.

Figure 5.

Scaling coefficient $a_d$ in (4.4) as a function of $Re$ . The dashed line is a least squares linear fit $\ln (a) = -0.989 \ln (Re) + 2.584$ .

Figure 5 shows that the dependence of $a_d$ on $Re$ is well described by the relation

(4.6) \begin{equation} a_d = 13.25\, Re^{-0.989}, \end{equation}

i.e. $a_d$ is essentially inversely proportional to $Re$ . Comparing with (4.5) suggests that, at least for this flow, $\partial _{xx}\Lambda _0/\partial _{xxxx}\omega _0$ is independent of $Re$ .

Figure 6.

Radial advection velocity $v_{r,a}$ (solid lines) and radial viscous drift velocity $v_{r,d}$ (dashed lines) of vortex centre 1. Time of maximal inward advection $t_a$ is marked with circles. Time of maximal outward drift $t_d$ is marked with squares. Time is shown relative to the merging time $t_m$ : (a) low $Re$ , (b) high $Re$ , (c) zoom of (b) near $t_a$ and $t_d$ .

Figure 6 shows the radial velocity $v_r$ of vortex centre 1, decomposed into its radial advection velocity $v_{r,a}$ and its radial viscous drift velocity $v_{r,d}$ according to (2.9). Again, there are common features for all $Re$ . The radial advective velocity $v_{r,a}$ decreases monotonically from zero (with some minor initial oscillations for higher $Re$ ) until a minimum is reached. Subsequently, $v_{r,a}$ approaches zero again, monotonically for lower $Re$ , and with damped oscillations for higher $Re$ . We term the time where the global minimum of $v_{r,a}$ is reached the time of maximal inward advection $t_a$ . The radial drift velocity $v_{r,d}$ initially increases from zero to reach a small positive maximum. The viscous drift subsequently decreases to negative values, again possibly with some oscillations, and finally goes monotonically to $-\infty$ as the merging time $t_m$ is approached. We term the time when the global maximum of $v_{r,d}$ is reached the time of maximal outward drift $t_d$ .

Figure 7.

Relative time of maximal inward advection ${t_a}/t_m$ , and relative time of maximal outward drift ${t_d}/t_m$ , versus the Reynolds number.

The dependence of the key time instants $t_a$ and $t_d$ on $Re$ relative to the merging time $t_m$ is shown in figure 7. In general, ${t_d} \lt {t_a}$ . At low $Re$ , ${t_d}/t_m$ is much smaller than ${t_a}/t_m$ , but as $Re$ increases, the interval between the times become shorter. The fact that in the high- $Re$ range, $t_a$ and $t_d$ approach each other, is closely related to the development of an inflection point for $d(t)$ that occurs at time instant $t_i$ , where $\ddot d(t_i) = 0$ . Such inflection points appear in figure 4 for $Re=500$ and above. To see the connection, consider the second-order Taylor expansions of $v_{r,a}$ and $v_{r,d}$ based at $t_a$ and $t_d$ , respectively:

(4.7) \begin{align} v_{r,a}(t) &\approx v_{r,a}({t_a}) + \frac {1}{2}\ddot {v}_{r,a}({t_a})\,(t-{t_a})^2, \end{align}

(4.8) \begin{align} v_{r,d}(t) &\approx v_{r,d}({t_d}) + \frac {1}{2}\ddot {v}_{r,d}({t_d})\,(t-{t_d})^2. \end{align}

Both are valid simultaneously when $t_a$ and $t_d$ are close. Since $\ddot {d}(t) = 2(\dot {v}_{r,a}(t)+\dot {v}_{r,d}(t))$ , one finds from (4.7) and (4.8) that $\ddot {d}(t)=0$ is solved by a weighted average of $t_a$ and $t_d$ :

(4.9) \begin{equation} t_i \approx \frac {\ddot {v}_{r,a}({t_a})\,{t_a} + \ddot {v}_{r,d}({t_d})\,{t_d}}{\ddot {v}_{r,a}({t_a}) + \ddot {v}_{r,d}({t_d})}. \end{equation}

It appears from figure 6(c) that $|\ddot {v}_{r,d}({t_d})| \ll |\ddot {v}_{r,a}({t_a})|$ , and it follows that $t_i$ is located between $t_d$ and $t_a$ .

Figure 8.

Angle $\theta$ of the line connecting the two vortices. The black dashed line shows the evolution for point vortices according to (4.10). The red dashed lines are expressions (4.11). The triangles mark the inflection times $t_i$ at $Re = 500, 100, 1500$ .

We now turn to the evolution of the angle $\theta$ of the line connecting the two vortex centres. The results are shown in figure 8. Two point vortices with unit distance, each of circulation $1/2$ , give rise to a constant angular velocity such that

(4.10) \begin{equation} \theta = \frac {1}{2\pi }t, \end{equation}

independent of $Re$ . It appears from figure 8 that this describes the initial dynamics well for all $Re$ .

As the bifurcation is approached, the scaling (3.25) also implies a constant angular velocity, but now depending on $Re$ . Inserting (3.22) into (3.25), an expression of the form

(4.11) \begin{equation} \theta = a_{\theta }t + b_{\theta }, \end{equation}

analogous to (4.4), can be derived. Similar to (4.5), $a_{\theta }$ and $b_{\theta }$ depend on derivatives of the flow field evaluated at the bifurcation point in a complicated way. The relation (4.11) is shown as red dashed lines in figure 8 for the higher $Re$ , where they deviate significantly from the point vortex scaling (4.10); it appears that the theoretical result agrees very well with the simulations. It is interesting that at high $Re$ , the transition from the initial scaling (4.10) to the bifurcation scaling (4.11) takes place over a very narrow time interval around $t_i$ that is defined solely by the properties of the radial velocity.

A decomposition of the angular velocity into advection and viscous drift, $\dot \theta = \dot \theta _a + \dot \theta _d$ , shows that $\dot \theta _a$ is generally two orders of magnitude larger than $\dot \theta _d$ . Hence advection completely dominates the rotation of the vortex pair.

4.3. Phases of the merging process

We summarise the results of § 4.2 by describing the merging process as consisting of a series of phases. Each phase has a characteristic balance between advection $v_{r,a}$ and viscous drift $v_{r,d}$ . We use the term ‘phases’ to avoid confusion with the ‘stages’ defined by Cerretelli & Williamson (Reference Cerretelli and Williamson2003) on the basis of the dynamics of $d(t)$ .

The initial constant phase is characterised by both advection and viscous drift being small, which results in $d(t)$ being almost constant. The final bifurcation phase is characterised by viscous drift being dominating as $v_{r,d}\rightarrow -\infty$ for $t \rightarrow t_m$ with $d(t)$ following the scaling (4.4).

For $Re$ up to 100, where $t_d$ is much smaller than $t_a$ (see figure 7), the constant phase blends smoothly into the bifurcation phase before any significant advection has occurred; see figure 6(a).

For $Re = 500$ and above, where $t_d$ is close to $t_a$ and an inflection time $t_i$ exists, there is sufficient time for inward advection to grow while viscous drift is still negligible. This gives rise to an advective phase where the vortex centres move inwards due to advection. The constant phase now blends into the advective phase, and while we do not define a specific time instant where it starts, we formally define that the phase ends at $t_a$ , where $v_{r,a}$ starts damped oscillations; see figures 6(b, c). The oscillations of $v_{r,a}$ persist until merging, and as long as the viscous drift is not too large, oscillations of $d(t)$ ensue. We term the phase after $t_a$ where $d(t)$ oscillates the oscillatory phase, which continues until the monotonicallly decreasing viscous drift ultimately dominates, and the bifurcation phase is entered.

An overview of the phases and the resulting dynamics is shown in table 1. The last column shows the angular velocity $\dot \theta$ during each phase, according to the results shown in figure 8.

Table 1.

Overview of merging phases, classified by the variation of advection and viscous drift.

We note that in their experiments at $Re = 530$ , Cerretelli & Williamson (Reference Cerretelli and Williamson2003) also observed an early constant phase, termed the first diffusive stage, which after a well-defined transition was followed by a regime termed the convective stage, where $d(t)$ varies linearly with $t$ . This is in contrast to our results, where $d(t)$ varies smoothly as advection builds up. Similar smooth development of $d(t)$ is observed in the simulations by Josserand & Rossi (Reference Josserand and Rossi2007) and Sreejith & Anil (Reference Sreejith and Anil2021). However, it is interesting to note that the occurrence of an inflection point in $d(t)$ at the beginning of what we defined as the oscillatory phase implies the existence of a regime where $d(t) = d(t_i) + \dot d(t_i)\,(t-t_i) + o((t-t_i)^2)$ . Mathematically, this implies that in the vicinity of $t_i$ , $d(t)$ varies almost linearly, and while clearly visible in figure 4, this regime is much shorter than the convective stage found experimentally. This suggests that the experimentally observed behaviour must arise via mechanisms not included in our model, e.g. three-dimensional effects.

As mentioned in the Introduction, Cerretelli & Williamson (Reference Cerretelli and Williamson2003) argue that the anti-symmetric vorticity in the co-rotating frame is responsible for the radial advection of the vortices. It is important to note that this gives a complete description of the motion of the vortex centres only if the viscous drift is negligible. Our analysis shows that this holds only in the advective phase, which exists only for sufficiently high $Re$ . In the oscillatory phase and the bifurcation phase, viscous drift cannot be ignored, and for low $Re$ , the vortex motion is completely dominated by viscous drift.

4.4. Dynamics of anti-vortices

Figure 9.

Trajectories of critical points of vorticity at $Re=1000$ . The black and dotted green curves are the trajectories of the vortices (local maxima), with the initial positions marked with filled circles. Red curves are anti-vortices (local minima), and blue curves are saddle points. Open circles mark bifurcations: purple for creation, and black for annihilation or vortex merging.

The extrema of vorticity that we have described hitherto are positive maxima that we interpreted as vortex centres. At low $Re$ , no further vorticity extrema occur, but at higher $Re$ , the vorticity field deforms sufficiently for other critical points to be created through bifurcations. We show here that at $Re=1000$ , local minima of positive vorticity, termed anti-vortices in § 2.2, can develop through cusp bifurcations as described in § 3.1. The generic case is relevant here, as there is no symmetry around the points where the bifurcations occurs. However, as the global vorticity field satisfies the symmetry (3.14), a bifurcation event at some $(x_0,y_0)$ must be accompanied by an identical event at $(-x_0,-y_0)$ at the same time.

Figure 9 shows the trajectories of all critical points of vorticity at $Re=1000$ . In addition to the original vortex centres and the saddle point at the origin, two extremum/saddle pairs are created at $t=35.8$ in cusp bifurcations (c) $\rightarrow$ (b) $\rightarrow$ (a) in figure 1, and they are annihilated again at $t=60.9$ in the reverse bifurcation. It follows from the analysis in § 3.1 that close to the bifurcation, the extremum and the saddle essentially move along the eigendirection $\boldsymbol{e}$ of the zero eigenvalue of the Hessian that defines the $x$ -coordinate used in the analysis. The two points move in opposite directions, and their velocities go to infinity as the bifurcation is approached. Under the generic assumption that the azimuthal direction ${\boldsymbol{e}}_{\theta }$ (see figure 3 a) at the bifurcation point is not collinear with $\boldsymbol{e}$ , it follows that one of the critical points will rotate with the fluid flow in the anticlockwise direction, and the other will rotate against the flow, when sufficiently close to the bifurcation. Figure 9 shows that the saddle points rotate against the flow after the creation, but only for a very short time. As annihilation is approached, the anti-vortices exhibit the retrograde motion.

The annihilation process is illustrated in figure 10. Note that the bifurcations take place quite far from the origin. At the time of annihilation, the vorticity has been stretched out into filaments, and figure 10(a) shows that the anti-vortices and their corresponding saddles are present in the troughs between the filaments.

The anti-vortices continue to exist after the merging at $t=55.6$ of the original vortices. Thus the axisymmetrisation of the single vortex left after merging is a topologically complex process, with an intermediate phase with two anti-vortices and two saddles.

Equations (3.11) predict that the distance $d_{es}$ between an extremum and a saddle close to bifurcation follows a scaling of the form

(4.12) \begin{equation} d_{es}^2 = a_{es}t + b_{es}, \end{equation}

similar to (4.4), valid for the distance between the extrema in vortex merging. The fits shown in figure 11 confirm this theoretical prediction.

Figure 10.

Snapshots of the vorticity field at $Re=1000$ : (a) just before anti-vortex annihilation at $t=60.9$ ; (b) just after annihilation.The green point marks the merged vortex at the origin, cyan points are anti-vortices, and red points are saddles. Typical iso-vorticity curves are shown in white. In (a), the iso-vorticity curve through the saddle points is also shown. As $t \rightarrow 60.9$ , the saddles and the corresponding anti-vortices approach each other.

Figure 11.

The distance $d_{es}$ between the critical points in an anti-vortex/saddle pair at $Re=1000$ . Quadratic fits of the form (4.12) near the bifurcations are shown as dashed lines.

Recall that (2.15) showed that the strength of a positive vortex, measured by the value of vorticity at its centre, will decrease over the course of time. The same equation implies that, somewhat counter-intuitively, the opposite is the case for a positive anti-vortex. This is confirmed in figure 12, which shows that an anti-vortex increases its strength by six orders of magnitude between the creation at $t=35.8$ and the destruction at $t=60.9$ . Nevertheless, the strength of the anti-vortex at the time of annihilation is $7.9\times 10^{-3}$ times the strength of the merged vortex at the origin, so it is still rather weak.

While the anti-vortices in the merging process considered here do not interact with the vortices, they may do so in flows that are more complex. In fact, it may be the saddles that are created together with anti-vortices in cusp bifurcations that have a physical significance. These saddles may merge with a true vortex in a subsequent cusp bifurcation, and hence remove that vortex from the flow. Bifurcations to that effect have been observed in a preliminary study of the vortex dynamics in the flow around a cylinder in a channel (Ozdemir Reference Ozdemir2023).

Figure 12.

Development of the vorticity at an anti-vortex over its lifetime, at $Re=1000$ .