2.1. Finding stationary solutions using a filament approach

2.1.1. Filament framework

We use the same vortex filament approach as in Durán Venegas & Le Dizès (Reference Durán Venegas and Le Dizès2019). Vorticity is considered to be concentrated along slender vortex filaments moving as material lines in the fluid according to

(2.1)\begin{equation} \frac{{{\rm d}} \boldsymbol{X}_j}{{{\rm d}}t} = \boldsymbol{U}(\boldsymbol{X}_j) - {U}_z^{\infty}\boldsymbol{\hat{e}_z}, \end{equation}

where $\boldsymbol {X}_j=(r_j,\theta _j,z_j)$ is the position of the $j$th vortex filament, ${U}_z^{\infty }$ an external velocity field and $\boldsymbol {U}(\boldsymbol {X}_j)=(U_r , r_j \varOmega, U_z)$ is the velocity induced by the vortices at $\boldsymbol {X}_j$. In a frame with rotation rate $\varOmega _R$, it is convenient to parametrize the position of vortex filaments in terms of two angular coordinates, $\zeta =\varOmega _R (t-t_0)$ and $\psi = \varOmega _R t$, and transform (2.1) as

(2.2)\begin{equation} \frac{\partial \boldsymbol{X}_j}{\partial \psi} + \frac{\partial \boldsymbol{X}_j}{\partial \zeta} = \frac{1}{\varOmega_R}\left[ \boldsymbol{U}(\boldsymbol{X}_j) + \boldsymbol{U}^{\infty}(\boldsymbol{X}_j) \right], \end{equation}

where $\boldsymbol {U}^{\infty }(\boldsymbol {X}_j)=(0,-r_j \varOmega _R, -{U}_z^{\infty })$ contains the contributions from the rotating frame and a constant free-stream velocity. In this context, $\zeta$ is often referred to as the wake age and corresponds to a spatial coordinate, while $\psi$ is a proxy of time (Leishman, Bhagwat & Bagai Reference Leishman, Bhagwat and Bagai2002). Filaments are discretized in straight segments $[\boldsymbol {X}^{n}_j, \boldsymbol {X}^{n+1}_j]$ in order to compute the velocity field using the Biot–Savart law. The equation is de-singularized using a cutoff approach with a Gaussian vorticity profile. Local contributions to the velocity field at $\boldsymbol {X}_j^{n}$ are obtained by replacing the straight segments by an arc of circle passing through $[\boldsymbol {X}^{n-1}_j, \boldsymbol {X}^{n}_j, \boldsymbol {X}^{n+1}_j]$ and using the cutoff formula (see Durán Venegas & Le Dizès Reference Durán Venegas and Le Dizès2019).

For this work, we consider basic flow solutions that satisfy,

(2.3a–c)\begin{equation} \frac{{{\rm d}}r_j}{{{\rm d}}\zeta} = \frac{U_r(\boldsymbol{X}_j)}{\varOmega_R}, \quad \frac{{{\rm d}} \theta_j}{{{\rm d}}\zeta} = \frac{\varOmega (\boldsymbol{X}_j)-\varOmega_R}{\varOmega_R}, \quad \frac{{{\rm d}}z_j}{{{\rm d}}\zeta} = \frac{U_z(\boldsymbol{X}_j) - U^{\infty}_z}{\varOmega_R}. \end{equation}

These solutions are $\psi$-independent solutions to (2.2). They are stationary in the moving frame in the sense that flow remains tangent to the vortex structure at all times:

(2.4)\begin{equation} (\boldsymbol{U}(\boldsymbol{X}_j) + \boldsymbol{U}^{\infty}(\boldsymbol{X}_j)) \times \boldsymbol{T}_j = 0, \end{equation}

where $\boldsymbol {T}_j$ is the tangent unit vector. As shown in Castillo-Castellanos et al. (Reference Castillo-Castellanos, Le Dizès and Durán Venegas2021), the spatial evolution of the wake as we move away from the rotor plane is obtained by numerically solving (2.3a–c) with boundary conditions on the rotor, at $\zeta =0$, where the position of each vortex is prescribed and far-field boundary conditions at $\zeta \to \infty$. In the following, we focus on the wake geometry in the far field, where spatially periodic solutions are obtained.

2.1.2. Periodic solutions in the far field

The solution in the far field has been obtained in Castillo-Castellanos et al. (Reference Castillo-Castellanos, Le Dizès and Durán Venegas2021). It was shown that it is no longer uniform as for a single helix but it nevertheless exhibits a certain spatial periodicity. In addition to the azimuthal symmetry $\theta \to \theta + 2 {\rm \pi}/N$ for $N$ vortex pairs, solutions are invariant by the double operation $z \to z + h$ and $\theta \to \theta + 2{\rm \pi} /\beta$. The parameters $h$ and $\beta$ are chosen such that there is a single location in an axial period $h$ where the vortices of a given pair are at the same azimuth. This azimuth is taken to define the mean radius $R$ and the separation distance $d$ of each vortex pair. Each vortex pair tends to form a helical braid on a larger helical structure of radius $R$ and pitch $H$ as illustrated in figure 1. The pitch of this larger helical structure is directly related to $h$ and $\beta$ by

(2.5)\begin{equation} H= h \beta, \end{equation}

and related to the pitch of the vortex pair through

(2.6)\begin{equation} h_\tau = \sqrt{H^{2} + (2{\rm \pi} R)^{2}}/\beta. \end{equation}

Figure 1. (a) Geometric parameters of periodic solutions: separation distance $d$, radius $R$, axial pitches $H$ and $h$, core size $a$ and $\beta = H/h$. Representation of (b) $N=1$ pairs without hub vortex and (c) $N=2$ pairs with hub vortex for ($R^{*}=7$, $H^{*}/N=5.25$, $\beta =3/4$).

Depending on the value of $\beta$, the double-helix structure may describe (i) a leapfrog-type pattern, where vortices trade places every $1/\beta$ turns; (ii) a relatively sparse braid; or (iii) a dense ‘telephone cord’-type pattern. We are typically in situation (i) when $\beta < 1$, and in situation (iii) when $\beta > 10$. The parameter $\beta$ also characterizes the axial periodicity of the solution. It becomes axially periodic only if $\beta$ is a rational number $p/q$. In that case, it means that the double helix makes $q$ turns on itself as the large helix does $p$ turns. The axial period is thus $pH=qh$. In the following, we only consider rational values of $\beta$.

As soon as the vortex core size $a$ is fixed, the far field is then defined by 5 non-dimensional parameters which are

(2.7a–e)\begin{equation} R^{*} \equiv \frac{R}{d}, \quad H^{*} \equiv \frac{H}{d}, \quad \varepsilon^{*} \equiv \frac{a}{d},\quad \beta,\quad N. \end{equation}

The solution is obtained by solving the system (2.3a–c) with the prescribed symmetry. The frame velocities $(\varOmega _R,U^{\infty }_z)$ are unknown quantities. However, these quantities are proportional to the vortex circulation $\varGamma$ as it is the unique quantity of the vortex system involving time. For this reason, we can fix $\varGamma$ to 1.

As shown in Castillo-Castellanos et al. (Reference Castillo-Castellanos, Le Dizès and Durán Venegas2021), the problem can be treated as a nonlinear minimization problem using an iterative procedure. The convergence is rapid if we start for each pair from an initial guess given by an undeformed double helix on a larger helix and estimates for $(\varOmega _R,U^{\infty }_z)$ obtained from uniform helices. The converged solution is found to exhibit spatial variations but in most cases the initial guess turns out to be a good approximation of the solution.

In the present study, we consider two different configurations: one composed of a single helical pair ($N=1$) without a central hub vortex (figure 1b), and another composed of $N=2$ helical pairs with a central hub vortex (figure 1c). Also, we fix $\varepsilon ^{*}=0.1$ and vary the remaining parameters ($R^{*}$, $H^{*}$ and $\beta$). In the following, solutions that satisfy (2.3a–c) are denoted $\boldsymbol {X}_j^{B}(\zeta )$.

2.2. Inviscid global instability analysis

The stability of $\boldsymbol {X}_j^{B}$ is analysed by considering the evolution of infinitesimal perturbation displacements

(2.8)\begin{equation} \boldsymbol{X}_j(\zeta,\psi) = \boldsymbol{X}_j^{B}(\zeta) + \boldsymbol{X}'_j(\zeta,\psi). \end{equation}

Equation (2.2) is linearized around $\boldsymbol {X}_j^{B}(\zeta )$ to obtain a linear dynamical system

(2.9a)$$\begin{gather} \frac{\partial r_i'}{\partial\psi} ={-}\frac{\partial r_i'}{\partial\zeta}+ \frac{1}{\varOmega_R}\sum_{j=1}^{N}\left[\frac{\partial U_r(\vec{X}_i^{B})}{\partial r_j} r_j' + \frac{\partial U_r(\vec{X}_i^{B})}{\partial\theta_j} \theta'_j + \frac{\partial U_r(\vec{X}_i^{B}) }{\partial z_j} z'_j \right] \end{gather}$$

(2.9b)$$\begin{gather}\frac{\partial \theta_i'}{\partial\psi} ={-}\frac{\partial \theta_i'}{\partial\zeta} + \frac{1}{\varOmega_R}\sum_{j=1}^{N}\left[\frac{\partial \varOmega(\vec{X}_i^{B}) }{\partial r_j} r_j' + \frac{\partial \varOmega(\vec{X}_i^{B}) }{\partial \theta_j} \theta'_j + \frac{\partial \varOmega(\vec{X}_i^{B}) }{\partial z_j} z'_j\right] \end{gather}$$

(2.9c)$$\begin{gather}\frac{\partial z_i'}{\partial\psi} ={-}\frac{\partial z_i'}{\partial\zeta} + \frac{1}{\varOmega_R}\sum_{j=1}^{N} \left[\frac{\partial U_z(\vec{X}_i^{B})}{\partial r_j} r_j' + \frac{\partial U_z(\vec{X}_i^{B})}{\partial \theta_j} \theta'_j + \frac{\partial U_z(\vec{X}_i^{B}) }{\partial z_j} z'_j\right], \end{gather}$$

which has the following general form:

(2.10)\begin{equation} \frac{\partial\boldsymbol{q}'}{\partial\psi} = \boldsymbol{\mathsf{L}}(\zeta)\boldsymbol{q}'(\zeta, \psi), \end{equation}

where $\boldsymbol {q}'=(\boldsymbol {X}_1, \boldsymbol {X}_2, \ldots )$ is the total displacement vector, and $\boldsymbol{\mathsf{L}}$ is an $N$ by $N$ block matrix containing the Jacobian terms. The spatial derivative (on $\zeta$) is evaluated using the same finite differences scheme as the base solution, while the velocity gradient is derived from the discretized Biot–Savart equations.

Each sub-matrix of $\boldsymbol{\mathsf{L}}$ can be written as

(2.11)\begin{equation} \boldsymbol{\mathsf{L}}_{ij} ={-}\delta_{ij}\frac{\partial}{\partial \zeta} + \frac{1}{\varOmega_R}\frac{\partial(\boldsymbol{U}(\boldsymbol{X}_i^{B}))}{\partial(\boldsymbol{X}_j)} \end{equation}

and contains $3n_s$ by $3n_s$ elements, where $n_s$ is the total number of discretization points of each vortex. The base flow is periodic with respect to $\zeta$ with a period $\zeta _B$ that satisfies $z_j(\zeta _B) = h$ and $\theta _j(\zeta _B)=2{\rm \pi} /\beta$. The difference between $\zeta _B$ and $\theta _j(\zeta _B)$ comes from the induced angular velocity $\varOmega (\boldsymbol {X}_j)$ in (2.3b). For small values of $\beta$, the ratio $C_\theta \equiv \zeta _B/\theta _j(\zeta _B)$ is close to 1, but decreases as $\beta$ increases (figure 2a,b). Each vortex contains 96 points per period $\zeta _B$. For the perturbations we consider a longer domain with a period $\zeta _p = n_p \zeta _B$ with $n_p = 32\beta$ such that $z_j(\zeta _p) = n_p h = 32 H$ and $\theta _j(\zeta _p)= n_p 2{\rm \pi} /\beta = 64{\rm \pi}$. This gives $n_p=24$ and $n_p=256$ for the values $\beta =3/4$ and $\beta =8$. For densely braided wakes ($\beta >20$), we use a shorter domain with $n_p = 4\beta$, such that $z_j(\zeta _p) = n_p h = 4 H$ and $\theta _j(\zeta _p)= n_p 2{\rm \pi} /\beta = 8{\rm \pi}$. This gives $n_p=84$ and $n_p=128$ for the values $\beta =21$ and $\beta =32$. Our objective is to apply periodic boundary conditions to the perturbations. We therefore discretize the operators $\boldsymbol{\mathsf{L}}_{ij}$ such that they satisfy this property.

Figure 2. Evolution of the ratio $C_\theta \equiv \zeta _B/\theta _j(\zeta _B)$ for ($N=1$, $R^{*}=7$): (a) $C_\theta$ as a function $H^{*}/R^{*}$ for $\beta =3/4$, and (b) as a function of $\beta$ for different $H^{*}/R^{*}$.

The final problem reduces to a linear system with constant coefficients that admits the formal solution

(2.12)\begin{equation} \boldsymbol{q}'(\zeta,\psi) = \exp(\boldsymbol{\mathsf{L}}\psi)\boldsymbol{q}_0'(\zeta) \end{equation}

for any initial perturbation $\boldsymbol {q}'_0(\zeta )=\boldsymbol {q}'(\zeta,0)$. In standard fashion, treating (2.12) as an eigenvalue problem can be used to describe the asymptotic limit $\psi \to \infty$. Here, because the domain is long but periodic, it also means that we consider perturbations that span the whole calculation domain.

We may decompose the linear operator as

(2.13)\begin{equation} \boldsymbol{\mathsf{L}} \boldsymbol{\varPhi} = \boldsymbol{\varPhi}\boldsymbol{\varLambda}, \end{equation}

where $\boldsymbol {\varPhi }$ is a matrix whose columns $(\boldsymbol \phi _1, \boldsymbol \phi _2, \ldots )$ are the eigenvectors of $\boldsymbol{\mathsf{L}}$, and $\boldsymbol {\varLambda }$ is a matrix whose diagonal entries are the corresponding (complex) eigenvalues $(\lambda _1, \lambda _2, \ldots )$. Here, $\sigma _m \equiv \mbox {Re}(\lambda _m)$ characterizes the temporal growth (or decay) of the perturbations, while $\omega _m\equiv \mbox {Im}(\lambda _m)$ characterizes the temporal oscillations. We may also introduce a dimensionless growth rate $\sigma _m^{*}$ and frequency $\omega _m^{*}$, based on the characteristic advection time scale of helical pairs $t_{adv} = \varGamma /(H^{2}/N^{2})$. Since $\boldsymbol{\mathsf{L}}$ is real valued, modes are either real or come in conjugate pairs. Equation (2.12) can be written as

(2.14)\begin{equation} \boldsymbol{q}'(\zeta,\psi) = \sum_m \boldsymbol{\phi}_m (\zeta) \exp{\lbrace\sigma_m\psi + {\rm i} \omega_m \psi\rbrace} b_m + c.c., \end{equation}

where the coefficients $b_m$ correspond to the coordinates of $\boldsymbol {q}'_0(\zeta ) = \boldsymbol {q}'(\zeta,0)$ expressed in the eigenvector basis, and $c.c.$ indicates the complex conjugate.

By construction, eigenvectors have the same dimensionality as $\boldsymbol {q'}$

(2.15)\begin{equation} \boldsymbol{\phi}_m = (\boldsymbol\phi_m^{(1)}, \boldsymbol\phi_{m}^{(2)}, \ldots), \end{equation}

where $\boldsymbol \phi _{m}^{(j)}$ is the $m$th eigenvector of the $j$th vortex filament. These eigenvectors represent a displacement vector and can be expressed in any coordinate system, e.g.

(2.16)\begin{equation} \boldsymbol\phi_{m}^{(j)} = (\check{r}_{m}^{(j)},\check{\theta}_{m}^{(j)},\check{z}_{m}^{(j)}) \end{equation}

indicates the components of $\boldsymbol \phi _{m}^{(j)}$ in global cylindrical coordinates. Conversely, displacement perturbations may also be expressed in terms of the local radial, azimuthal and axial coordinates associated with a uniform helix $\mathscr {H}$ of radius $R$ and pitch $H$,

(2.17)\begin{equation} \boldsymbol\phi_{m}^{(j)} = (\check{\rho}_{m}^{(j)}, \check{\phi}_{m}^{(j)}, \check{s}_{m}^{(j)}). \end{equation}

Due to the spatial periodicity, the eigenvectors can be expanded on a discrete Fourier basis

(2.18)\begin{equation} \boldsymbol{\phi}_m= \sum_{n=0}^{n_s-1} \boldsymbol{\hat{q}}_{mn} \exp{\left\lbrace \frac{\mathrm{i} 2{\rm \pi} n\zeta}{\zeta_p} \right\rbrace}, \end{equation}

where the azimuthal wavenumber $(k = n {\rm \Delta} k, {\rm \Delta} k = \beta /n_p)$ is normalized to ensure that $k=1$ corresponds to one helix turn, and $n_s$ the number of points used to discretize each vortex.

4.1. Typical displacement modes for leapfrogging wakes

For each pair, we introduce the following decomposition:

(4.1a,b)\begin{equation} \boldsymbol{X}_+\equiv (\boldsymbol{X}_1+\boldsymbol{X}_2)/2,\quad \boldsymbol{X}_-\equiv (\boldsymbol{X}_1-\boldsymbol{X}_2)/2, \end{equation}

where $\boldsymbol {X}_+$ characterizes the large-scale pattern traced by the vorticity barycentre, while $\boldsymbol {X}_-$ represents the rotation of the vortex pair relative to $\boldsymbol {X}_+$. Figure 1(b) depicts a leapfrogging wake with $\beta =p/q=3/4$, where $\boldsymbol {X}_+$ is represented as a tube enclosing the vortex pair. Over a single period, the vortex pair completes $p=3$ rotations around $\mathscr {H}$, while $\mathscr {H}$ completes $q=4$ rotations around the $z$ axis. Note that $\boldsymbol {X}_+$ and $\mathscr {H}$ are close but not exactly equal due to the effect of self-induction. In a similar vein, we introduce the following decomposition:

(4.2a,b)\begin{equation} \boldsymbol\phi_m^{+} \equiv (\boldsymbol\phi_m^{(1)} + \boldsymbol\phi_m^{(2)})/2,\quad \boldsymbol\phi_m^{-} \equiv (\boldsymbol\phi_m^{(1)} - \boldsymbol\phi_m^{(2)})/2, \end{equation}

where $\boldsymbol \phi _m^{+}$ and $\boldsymbol \phi _m^{-}$ indicate the displacement modes of $\boldsymbol {X}_+$ and $\boldsymbol {X}_-$, respectively.

Figure 5 presents the frequency spectrum, which is characterized by a set of modes distributed over three contiguous lobes at low frequencies, and a second set of modes over an additional lobe at higher frequencies. The maximum growth rate is observed near $\omega ^{*}=4.8$, with additional local maxima (in descending order) near $\omega ^{*}=51.9$, $\omega ^{*}=14.5$ and $\omega ^{*}=24.0$. These values, respectively, correspond to dominant wavenumbers $k_0=0.5$, $k_0=5.4$, $k_0=1.5$ and $k_0=2.5$. Dimensionless growth rates are slightly larger than the equivalent helical vortex (figure 5b). This can be explained by the change in the effective distance separating neighbouring loops.

Figure 5. Stability curves of a leapfrogging wake ($R^{*}=7$, $H^{*}=5.25$, $\beta =3/4$): dimensionless growth rate $\sigma ^{*}$ as a function of (a) $\omega ^{*}$ and (b) $k_0$. For reference, (b) shows the prediction for a uniform helix ($H/R=0.75$, $a_e/R=0.1$) in red lines. A schematic representation of the two groups of modes is shown in (c), where the symmetric (respectively anti-symmetric) part is shown in red (respectively green) arrows.

Low-frequency modes are clearly reminiscent of the unstable modes for the equivalent uniform helices, while those at higher frequencies are specific to this geometry. The two groups differ in the relative alignment between the displacements of the pair: predominantly aligned displacements (or symmetric with respect to the vorticity barycentre, figure 5(c) top) for low-frequency modes and predominantly opposed displacements (or anti-symmetric with respect to the vorticity barycentre, figure 5(c) bottom) for those at higher frequencies. A more quantitative way to illustrate this difference is through the spatial cross-correlation

(4.3)\begin{equation} (z_i' \star z_j')(\varDelta\theta) \equiv \int_\infty^{\infty} z_i'(\zeta/C_\theta) z_j'(\zeta/C_\theta + \varDelta\theta)\,{{\rm d}}\zeta, \end{equation}

where $\varDelta \theta$ is the delay in angular position. For instance, for $S_2$ the displacements between neighbouring turns, i.e. $\varDelta \theta =2{\rm \pi}$, are well anti-correlated since perturbations are in opposition to the phase (figure 6a,c). Conversely, for mode $A_1$ the auto-correlation between consecutive turns is negative, while the cross-correlation is positive (figure 6b,d), meaning that vortices move in opposite directions, but one of them is aligned with one of the vortices in the neighbouring loop.

Figure 6. Functions $(z_1'\star z_1')$ (a,b) and $(z_1'\star z_2')$ (c,d) for the case in figure 1(b) and modes: (a,c) $S_2$, and (b,d) $A_1$.

Figure 7 shows the eigenvectors $\check {s}_{m}^{+}$ and $\check {s}_{m}^{-}$ corresponding to the most unstable modes $S_1$ and $A_1$. Here, $\check {s}_{m}^{+}$ (respectively $\check {s}_{m}^{-}$) corresponds to the component of $\boldsymbol \phi _m^{+}$ (respectively $\boldsymbol \phi _m^{-}$) along the local axial direction. For the symmetric mode, $\check {s}_{m}^{+}$ is dominant and has a nearly constant envelope, while $\check {s}_{m}^{-}$ has a sinusoidal envelope, whereas the anti-symmetric mode displays the opposite behaviour.

Figure 7. Complex eigenvectors for the case in figure 5. Perturbations in the local axial direction (a,c) $\check {s}_{m}^{+}$, and (b,d) $\check {s}_{m}^{-}$ shown for modes (a,b) $S_1$, and (c,d) $A_1$.

Each mode displays a dominant wavenumber $k_0$ with additional peaks (typically in descending order of magnitude) at wavenumbers $k_n \approx k_0 \pm n \beta$ for integer values of $n$ (see Fourier spectra in figure 8(a–c) and table 2). This pattern is also observed in the spectra of the axial and radial components, $\check {z}_{m}$ and $\check {r}_{m}$. From these observations, we infer that perturbations propagate along the structure as

(4.4)\begin{align} s'_j(\zeta,\psi) \approx \sum_m \left[\left[\hat{s}_{m0} + \sum_{n=1} \hat{s}_{mn} \cos{(n\beta\zeta/C_\theta + \widehat{\varphi}_{mn})}\right] \cos(\omega_m \psi + k_0\zeta/C_\theta)\,{\rm e}^{\sigma_m \psi} \right], \end{align}

where $\hat {s}_{mn}$ and $\widehat {\varphi }_{mn}$ are the amplitudes and phase differences measured with respect to $k_0$. Consider the leading terms in (4.4). For symmetric modes $s_+^{\prime }$ roughly corresponds to a travelling wave, whereas $s_-^{\prime }$ approximates a wave that is modulated in amplitude by a cosine function of period $\zeta _B$, i.e. the periodicity of the base flow (figure 8a,b). Anti-symmetric modes display the opposite behaviour: $s_-^{\prime }$ approximates a travelling wave, whereas $s_+^{\prime }$ is modulated in amplitude by a function of period $\zeta _B$ (figure 8c). In both cases, the contributions from higher-order terms also correspond to waves modulated in amplitude by multiples of $\zeta _B$ in decreasing order of magnitude.

Figure 8. Fourier spectrum of modes: (a) $S_1$, (b) $S_2$ and (c) $A_1$, for the case in figure 5.

Table 2. Leading wavenumbers of the modes shown in figure 8. Here, ${\rm \Delta} k = 0.03125$.

Figure 9(a) presents the deformation due to the symmetric mode $S_1$ by plotting the perturbed geometry $\boldsymbol {X}_j=\boldsymbol {X}_j^{B} + \boldsymbol {X}_j'$ for some arbitrary amplitude. A developed plan view illustrates $p=1$ localized pairing events for every $q=2$ neighbouring turns of the large-scale pattern (seen as dashed lines in figure 9a,b). An additional example corresponding to mode $S_2$ is shown in figure 9(c,d), where $p=3$ localized pairing events are observed every $q=2$ neighbouring turns. This behaviour was expected since the predominantly aligned displacements result in a block displacement of the vortex pair. As a result, the large-scale pattern behaves like a uniform helix where perturbations with wavenumber $k_0=p/q$ repeat after $p$ cycles and display local pairing events at $q$ azimuthal locations (Widnall Reference Widnall1972).

Figure 9. Base state from figure 1(b) perturbed by modes (a,b) $S_1$ and (c,d) $S_2$. (a,c) Three-dimensional and (b,d) developed plan views illustrating the local pairing modes.

Anti-symmetric modes behave in a different manner. Here, the two vortices move towards (or away from) one another such that $\boldsymbol {X}_+$ deforms much less and only at specific positions (figure 10a). In other words, the pair predominantly displays an anti-symmetric motion with respect to the helical structure, hence the name. Displacements are localized and not necessarily aligned with the rotation of the pair. For instance, at the azimuths, where displacements are perpendicular to the line connecting the two vortices, the structure is twisted back and forth, whereas when displacements are parallel, the separation distance expands and contracts. The latter could potentially trigger the merging of the vortex pair. This localization can be deduced from the envelopes of the corresponding eigenvectors, illustrated as dashed lines enclosing each vortex in figure 10(b). If we consider a longitudinal cut, we can see that displacements of a given vortex are paired with one of the vortices in the neighbouring turn but not with its companion (see arrows at the bottom part of figure 10a). However, the choice of the characteristic wavenumber $k_0$ is not obvious. For instance, doubling $H^{*}$ from $5.25$ to $10.5$ shifts mode $A_1$ from $k_0=5.406$ to $k_0=5.938$, while increasing $\beta$ from $3/4$ to $7/2$ shifts the wavenumber to $k_0=5.552$. However, perturbations between consecutive pairs remain anti-correlated, suggesting the dominant wavenumber is selected by the geometry as the one that amplifies the local pairing. We shall explore this relation in the following section.

Figure 10. Base state from figure 1(b) perturbed by mode $A_1$. (a) Three-dimensional and (b) developed plan view. Arrows indicate the displacements at the centre plane in (a).

In Appendix A, we also analyse the long-wave instability of our solutions using the linear impulse response approach developed in Durán Venegas et al. (Reference Durán Venegas, Rieu and Le Dizès2021). The idea is to solve the linear perturbation equation (2.12) using a Dirac function as initial condition and analyse the spatio-temporal growth of the resulting wavepacket. A sufficiently long domain is considered so that the wavepacket does not reach the boundaries during the length of the simulations. As shown in Appendix A, the temporal modes can be recovered but their study is more difficult to perform with the linear impulse response as all the instability modes are simultaneously excited. However, the linear impulse response provides additional information by telling us how the instability spreads in space. In particular, we are able to show that, although the most unstable anti-symmetric perturbations propagate at a similar speed as the symmetric perturbations, their spreading in space is much less important. The transition from convective to absolute instability is therefore expected to be associated with the symmetric perturbations. Moreover, as for the linear spectrum, we also demonstrate that the spatio-temporal evolution of these symmetric perturbations can be well described by the spatio-temporal evolution of the linear impulse response on a uniform helical vortex of large core size.

4.2. Influence of $H^{*}$ and $R^{*}$ on the modes of leapfrogging wakes

Since the two kinds of pairing modes described in § 4.1 seem to involve neighbouring turns of $\mathscr {H}$, growth rates are expected to display a strong dependency on the relative pitch $H^{*}/R^{*}$. In general, the maximum growth rate is larger than the maximum growth rate obtained for uniform helices with the total circulation of the vortex pair (figure 11a), which is slightly larger than $\sigma =({\rm \pi} /2)(\varGamma /H^{2})$ obtained for an array of point vortices of circulation $2\varGamma$ and separated by a distance $H$. For the range of values considered, symmetric and anti-symmetric modes have comparable growth rates (figure 11a,b). For symmetric modes, $\sigma ^{*}$ initially increases for mode $S_1$ before vanishing, while modes $S_2$, $S_3$ and so on, gradually vanish as the pitch increases. For anti-symmetric modes, $\sigma ^{*}$ initially increases at a faster rate, but also vanishes more quickly, while the dominant wavenumber is approximated by $k_0=0.45H^{*}/R^{*}$ for $H^{*}/R^{*}>2$, see figure 11(c).

Figure 11. Growth rates and dominant wavenumbers for $\beta =3/4$. (a–c) shown as a function of $H^{*}/R^{*}$ for $R^{*}=7$; and (d–f) shown as a function of $R^{*}$ for $H^{*}/R^{*}=0.75$.

It is also interesting to fix the relative pitch $H^{*}/R^{*}$ and change the separation distance through $R^{*}\equiv R/d$ (figure 11d). Stability curves are essentially unchanged for modes $S_1$ and $S_2$, while the growth rate of $S_3$ is observed to decrease as the separation distance becomes smaller (figure 11c). This is also reminiscent of uniform helical vortices, where the maximum growth rates also decrease as the effective core size becomes smaller (see, for instance figure 4). For mode $A_1$, the growth rate remains constant (figure 11e), while $k_0$ seems to evolve linearly with $k_0 \sim 0.89R^{*}$ and overlaps the symmetric modes for $R^{*}<6$ (figure 11 f).

From this dependency on $H^{*}$ and $R^{*}$, we may deduce the following: (i) $\sigma$ scales with $\varGamma /H^{2}$ for modes $A_1$ and $S_1$ whatever $H^{*}$ and $R^{*}$. This is consistent with a local pairing acting over a distance comparable to $H$; and (ii) the wavenumber most amplified by $A_1$ increases linearly with $H^{*}/R^{*}$ (for constant $R^{*}$) and with $R^{*}$ (for constant $H^{*}/R^{*}$). In other words, mode $A_1$ deviates from the classical pairing of uniform helices, and amplifies a linear wavelength $l_0=2{\rm \pi} R k_0 \sim d^{-1}$ instead. This is reminiscent of a four-vortex system involving two co-rotating pairs separated by a distance $b=O(H)$ (see, for instance Crouch Reference Crouch1997; Fabre & Jacquin Reference Fabre and Jacquin2000). We shall revisit this matter in § 5 with the case of $N=2$ helical pairs.

4.3. Typical displacement modes for sparsely braided wakes

Figure 12(a) presents the stability curves as we move from leapfrogging to sparsely braided wakes ($1<\beta <10$). Modes $S_1$, $S_2$ and $A_1$, behave as leapfrogging wakes. Some modes may act on the large-scale pattern and on the distance separating the vortex pair (see, for instance mode $S_1$ in figure 13a,c,e). For these modes, small oscillations in $\sigma ^{*}(\beta )$ can be explained by a change in the effective distance between neighbouring vortices, which varies by a factor $d\cos (\phi )$ due to relative orientation of the vortex pairs, where $\phi = 2{\rm \pi} /q$ for $\beta =p/q$. For $S_1$ and $S_2$, the unstable frequencies and leading wavenumbers are unchanged, while the wavenumber most amplified by $A_1$ displays some dependency on $\beta$.

Figure 12. Stability of sparsely braided wakes for $R^{*}=10$ and $H^{*}/R^{*}=1$. (a) Growth rates and dominant wavenumbers; and (b) maximum growth rates as function of $\beta$.

Figure 13. Base state from figure 12 and $\beta =7$ perturbed by modes (a) $S_1$ and (b) $A_0$. (b) Inset includes the base state (in dashed lines) for reference.

We observe an additional set of low-frequency modes (figure 12a). Of this set of modes, the maximum growth corresponds to the zero-frequency mode $A_0$, which is almost a linear function of $\beta$ (figure 12b). The displacement produced by $A_0$ is characterized by a radial expansion of the vortex pair and a translation along the helical coordinate (see developed plan view in figure 13(b,d, f) and corresponding insets). The resulting perturbed state would correspond to a similar braided wake but one with slightly larger $d$ and $R$. Other modes in the same branch display a similar displacement, but one where radial expansion and the translation along $\mathscr {H}$ are modulated in amplitude by a wavenumber $k_0$. For the case of two uniform helices, an equivalent displacement would yield two uniform helices of slightly larger radius. In such a case, the initial displacement is not amplified and the system remains in neutral equilibrium. However, for helical braids, the perturbed state is not necessarily a solution of (2.4), explaining the positive growth rate. This branch is not always present and, given the size of the parameter space ($H^{*}$, $R^{*}$ and $\beta$), it is unclear how $\sigma ^{*}$ varies. For instance, for the case presented in figure 12(b) and $\beta =4$, $\sigma ^{*}=\sigma (H^{2}/\varGamma )$ is close to $1.5$. For the same $(R^{*}, \beta )$ and $H^{*}/R^{*}=2$, $\sigma ^{*}$ is nearly four times larger, suggesting these modes no longer scale with $\varGamma /H^{2}$. This behaviour extends to the case of densely braided wakes.

4.4. Typical displacement modes for densely braided wakes

Densely braided wakes correspond to the case $\beta > 10$, when the pitch $h_\tau$ becomes of the same order as the separation distance $d$ and the pairing modes of the vortex pair become dominant. The stability curves in figure 14(a,b), display nearly all of the modes introduced so far. Modes $S_1$, $S_2$ and $A_1$ remain unchanged but are dwarfed by the other modes. For $\beta =21$, the branch containing mode $A_0$ is still present, and now contains a new local maximum around $k_0/\beta =2$, denoted $\tilde {A}_2$ in figure 14(b). For $\beta =24$, the branch containing both modes splits in two. Modes $A_0$ and $\tilde {A}_2$ display a similar scaling and seem to vanish for large $\beta$ (figure 14c). The spatial structure of $A_0$ is unchanged (see, figure 13(b,d, f) insets and 15(a)), while mode $\tilde {A}_2$ displays a radial expansion modulated in amplitude with some spatial frequency (figure 15b).

Figure 14. Stability of densely braided wakes for $R^{*}=10$ and $H^{*}/R^{*}=1$: growth rate as a function of (a) $\omega$ and (b) $k_0$ for different $\beta$; and (c) growth rate of modes identified in (a,b) as a function of $\beta$.

Figure 15. Stability of densely braided wakes. Schematic representation of the displacement modes identified in figure 14: (a) $A_0$, (b) $\tilde {A}_2$, (c) $B_0$ and (d) $B_1$. For reference, each figure also includes the base state (in dashed lines).

Finally, we have the pairing modes of the vortex pair. For these modes, the scaling of the growth rate is different since the pairing acts over a distance comparable to $h_\tau /2$ instead of $H$. The maximum growth rate is comparable to the predictions for two interlaced helices and $\sigma = ({\rm \pi} /2) (2\varGamma /h_\tau ^{2})$ obtained for an array of point vortices of circulation $\varGamma$ and separated by $h_\tau /2$ (figure 14c). Here, mode $B_0$ corresponds to a special case with $\omega =0$ and $k_0=0$. As shown in figure 15(c), displacements are characterized by a radial expansion of the vortex pair and a translation along the helical coordinate for one vortex, and a radial contraction and a translation in the opposite direction for the other one. This results in a form of uniform pairing along the helical coordinate, analogous to the global pairing mode of two uniform helices. For $\beta >25$, a second maximum, denoted $B_1$, is observed near $k_0/\beta = 1$. As shown in figure 15(d), displacements approach the vortices in neighbouring turns at specific intervals, analogous to the local pairing mode of two uniform helices.

5.1. Typical displacement modes for leapfrogging wakes

We consider similar geometric parameters as in § 4.1, but increase the pitch so as to preserve the mean axial distance between neighbouring turns. The results are summarized in figures 16–18. The frequency spectra presented in figure 16(a,c) are characterized by a first set of modes distributed over three overlapping lobes at low frequencies, and a second set distributed over two additional lobes at higher frequencies and a small lobe containing $A_0$. As in the previous case, the former are symmetric modes, the latter are anti-symmetric. The structure of the eigenvectors is the same as before with dominant wavenumbers at $k_0$ and additional harmonic terms at $k_n=k_0 \pm n \beta$, see table 3. In general, the maximum growth rate is larger than the maximum growth rate obtained for two uniform helices with the total circulation of each vortex pair, which is slightly larger than $\sigma = ({\rm \pi} /2)(2\varGamma /H^{2})$ obtained for an array of point vortices of circulation $2\varGamma$ and separated by a distance $H/2$.

Figure 16. Stability curves for the case of two helical pairs (a) without and (b) with a central hub vortex and ($R^{*}=7$, $H^{*}/R^{*}=1.5$, $\beta =3/4$). Colour indicates the phase difference between the two helical pairs: $\varphi =0$ (in black), $\varphi ={\rm \pi}$ (in red). For reference, (a) include the values for two uniform helices ($H/R=1.5$, $a_e/R=0.1$) in dashed lines.

Table 3. Leading wavenumbers of the modes shown in figure 16. Here, ${\rm \Delta} k = 0.03125$.

As expected, symmetric modes are reminiscent of the unstable modes obtained for two equivalent helical vortices. Displacements between neighbouring turns are out of phase for modes in the branch containing $S_0$ and $S_2$ (in black), and phase aligned for modes in the branches containing $S_1$ and $S_3$ (in red). These displacements result in a local pairing at $2k$ azimuthal positions per turn of the large-scale helix. Mode $S_0$ corresponds to a special case. Displacements between neighbouring turns are out of phase, where one vortex pair expands in the radial direction, while the other one contracts, resulting in a uniform pairing of the large-scale pattern along the azimuthal direction (figure 17a,b). This is analogous to the global pairing mode $k=0$ of the interlaced helices (see, for instance Okulov & Sørensen Reference Okulov and Sørensen2009; Quaranta et al. Reference Quaranta, Brynjell-Rahkola, Leweke and Henningson2019).

Figure 17. (a) Three-dimensional view of the base state from figure 1(c) perturbed by mode $S_0$, and (b) developed plan view illustrating the global pairing.

Additionally, we note that anti-symmetric modes are observed over a similar range of values as in the case of one vortex pair. One notable difference is that two modes are now obtained for the same $k_0$ (each using a different colour in figure 16a,c), where each pair has a similar structure but shifted in phase (not shown).

The introduction a central hub vortex ensures the total circulation is zero and the angular velocity vanishes as $r \to \infty$. This has a small effect on the frame velocity and the base states with and without a central hub are qualitatively similar. The presence of a central hub modifies the stability properties to a small degree (figure 16b,d). For instance, the case with a central hub vortex has generally larger growth rates than the case without by 2 %–3 %. Additionally, the branches containing $A_0$ are suppressed by the hub vortex, while the out-of-phase perturbations are no longer neutrally stable near $k_0=1$. In the following, we consider only the case with a central hub.

5.2. Geometric dependency of the most unstable modes

As expected, symmetric and anti-symmetric modes display a different behaviour as we vary the geometric parameters. Varying the separation distance has a small influence on the symmetric modes. As in the case of one vortex pair, the change in growth rates is reminiscent of that of varying the effective core size in uniform helices: small for $S_0$ and $S_1$, but more important for higher wavenumbers (figure 18a). For anti-symmetric modes, the change in growth rates is also small (figure 18b), while the dominant wavenumber $k_0$ increases linearly with $R^{*}$ (figure 18c). For $H^{*}/R^{*}=1.5$, $k_0$ was found to be roughly the same as in the case of one vortex pair with equal effective pitch, suggesting the modes are selected by the same pairing mechanism described in § 4.

Figure 18. (a,d) Stability curves, (b,e) growth rates and (c, f) dominant wavenumbers as functions of $R^{*}$ for $(H^{*}/R^{*}=1.5,\beta =3/4)$ (a–c), and as functions of $H^{*}/R^{*}$ for $(R^{*}=7,\beta =3/4)$ (d–f).

Varying the relative pitch shows the transition between two different regimes. For small $H^{*}/R^{*}$, stability curves have a similar structure as before: symmetric modes distributed over two or more branches with small $k_0$ and anti-symmetric modes distributed over two overlapping branches with larger $k_0$. Growth rates are larger than predicted rates for the equivalent uniform helices, but still close to the point vortex prediction $\sigma \sim \varGamma N^{2}/H^{2}$. For large $H^{*}/R^{*}$, the stability curves are characterized by two branches at small $k_0$ containing modes $S_0$ and $S_1$, and a single branch containing $A_1$ at larger $k_0$. Modes $S_1$ and $A_1$ are shifted towards larger wavenumbers. Instead of vanishing, the dimensionless growth rate $\sigma ^{*}$ proceeds to increase, pointing to a different scaling law with $\sigma \sim \varGamma /R^{2}$ for large $H^{*}/R^{*}$ (figure 18d,e). A similar transition is observed for the most unstable wavenumber $k_0$: $k_0 \sim H^{*}/NR^{*}$ for small $H^{*}/R^{*}$, and $k_0 \sim H^{*}/R^{*}$ for large $H^{*}/R^{*}$ (figure 18 f).

This change of regime can be understood as follows. For the case of a single pair ($N=1$), the limit of large $H^{*}/R^{*}$ leads to a pair of parallel co-rotating vortices which are known to be stable with respect to long-wavelength perturbations (Jimenez Reference Jimenez1975). For the case of $N=2$ vortex pairs with a central hub, the same limit provides a system composed of two co-rotating pairs of vortices of circulation $\varGamma$ and one counter-rotating vortex of circulation $-4\varGamma$ at the centre. For this configuration, the instability is necessarily controlled by the distance between the vorticity centroids ($b=2R$ in our current notation), and the separation distance $d$. This configuration is similar to four-vortex systems (Crouch Reference Crouch1997; Fabre et al. Reference Fabre, Jacquin and Loof2002), which are known to be unstable, with a maximum growth rate scaling with $\varGamma /b^{2}$, and a most unstable wavenumber varying with $b$ and $d$.