# Hybrid patterns and solitonic frequency combs in non-Hermitian Kerr Cavities

Salim B. Ivars  Departament de Física, Universitat Politècnica de Catalunya
(UPC), Rambla Sant Nebridi 22, 08222, Terrassa, Barcelona, Catalonia, Spain
Carles Milián  Institut Universitari de Matemàtica Pura i Aplicada,
Universitat Politècnica de València, 46022 València, Spain  Muriel Botey
Departament de Física, Universitat Politècnica de Catalunya (UPC), Rambla Sant
Nebridi 22, 08222, Terrassa, Barcelona, Catalonia, Spain  Ramon Herrero
Departament de Física, Universitat Politècnica de Catalunya (UPC), Rambla Sant
Nebridi 22, 08222, Terrassa, Barcelona, Catalonia, Spain  Kestutis Staliunas
Departament de Física, Universitat Politècnica de Catalunya (UPC), Rambla Sant
Nebridi 22, 08222, Terrassa, Barcelona, Catalonia, Spain  Institució Catalana
de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, E-08010,
Barcelona, Spain  Vilnius University, Faculty of Physics, Laser Research
Center, Sauletekio Ave. 10, Vilnius, Lithuania

###### Abstract

We unveil a new scenario for the formation of dissipative localised structures
in nonlinear systems. Commonly, the formation of such structures arises from
the connection of a homogeneous steady state with either another homogeneous
solution or a pattern. Both scenarios, typically found in cavities with normal
and anomalous dispersion, respectively, exhibit unique fingerprints and
particular features that characterise their behaviour. However, we show that
the introduction of a periodic non-Hermitian modulation in Kerr cavities
hybridises the two established soliton formation mechanisms, embodying the
particular fingerprints of both. In the resulting novel scenario, the
stationary states acquire a dual behaviour, playing the role that was
unambiguously attributed to either homogeneous states or patterns. These
fundamental findings have profound practical implications for frequency comb
generation, introducing unprecedented reversible mechanisms for real-time
manipulation.

Pattern formation in extended systems is a far-from-equilibrium phenomenon
that rules the dynamics of nonlinear physical, chemical or biological systems
Cross and Hohenberg (1993); Arecchi et al. (1999); Staliunas and Sanchez-
Morcillo (2003); Lugiato et al. (2015); Turing (1990); Knobloch (2015). The
existence and stability of pattern states (PS), together with that of the so-
called flat or homogeneous states (HS), lies at the core of the dissipative
localised structure (LS) formation Schäpers et al. (2000); Meron (1992);
Pomeau (1986), including optical solitons and frequency combs. The emergence
of solitons in different regimes of nonlinear systems is generally associated
with the existence of two stable solutions, either a PS and HS for anomalous
group velocity dispersion (GVD) or two stable HSs for normal GVD. Some non-
linear systems have additionally shown to hold tristability: three different
stable HS Dumeige and Féron (2011); Anderson et al. (2017); Milián et al.
(2018a), two PS and one HS Tlidi et al. (1993); Spinelli et al. (1998) or two
HS and one PS Zelnik et al. (2018); Al Saadi and Champneys (2021); Parra-Rivas
et al. (2022). However, the two paradigms of formation of LS do not coexist
and the character of the basal states is unique and well-defined.

In this letter, we unveil a novel hybrid scenario in the formation of solitons
and patterns opening a new avenue in the comprehension of nonlinear systems.
Specifically, we demonstrate that the introduction of a periodic modulation in
Kerr cavities with normal GVD is qualitatively similar to the effect of Turing
instability of the anomalous GVD, in the sense that the modulation introduces
low amplitude Turing-like (periodic) states without destabilising the system.
This is the key for the dual and hybridised mechanism that we demonstrate
below, which is understood as the blending of the two traditionally separated
LS forming regimes. As a result, new families of stable solitons, molecules,
and patterns emerge in the normal GVD regime.

In the last two decades, the damped-driven Nonlinear Schrödinger Equation
(NSE) has attracted significant attention Kaup and Newell (1978); Barashenkov
and Smirnov (1996); Chembo and Menyuk (2013); Lugiato and Lefever (1987);
Haelterman et al. (1992), for the mastering of microcavities and frequency
comb generationBao et al. (2020); Brasch et al. (2016); Kippenberg et al.
(2011). This model has been demonstrated to accurately describe the evolution
of the light field in Kerr resonators. It admits LSs in the two canonical
formation regimes associated with normal and anomalous GVD Matsko et al.
(2012); Barashenkov and Smirnov (1996).

The introduction of inhomogeneities in the damped-driven NSE has been studied
through the modulation of external injection or intracavity modulation Copie
et al. (2016). In the first case, the soliton dynamics has been studied for
temporal modulations of the driving field Todd et al. (2023); Hendry et al.
(2018), and the introduction of a spatial profile in the injected field has
shown to stabilise LSs Tabbert et al. (2019). In the second case, phase
modulations by Electro-Optical Modulators within the ring cavity have shown to
support synthetic dimensions Tusnin et al. (2020); Englebert et al. (2023);
Dutt et al. (2020) and stabilisation of 3D solitons Sun et al. (2023).

In parallel, the introduction of non-Hermitian (complex) potentials in non-
linear systems has demonstrated the ability to induce a wide range of
intriguing properties: unidirectional couplings arising from potential
asymmetries Ivars et al. (2022); Akhter et al. (2023), stabilisation of new
solutions Ivars et al. (2023a); Ahmed et al. (2018), the support of constant
intensity waves Makris et al. (2015), the occurrence of exceptional points and
jamming anomaly Barashenkov et al. (2016) or selective single-mode lasing Feng
et al. (2014), among others (see, e.g., Ashida et al. (2020); Konotop et al.
(2016); Suchkov et al. (2016) for reviews). In optics, non-Hermitian
potentials have been successfully implemented in recent experiments, Özdemir
et al. (2019); Cheng et al. (2023). For instance, EOM have shown to
effectively modulate the complex refractive index with integrated technology
(see Fig. 1(a) and Hu et al. (2024) for a recent review, and references
therein). In thin-film Fabry-Perot resonator Lugiato and Prati (2010); Hachair
et al. (2004); Brambilla et al. (1996); Peschel et al. (2003); Brambilla et
al. (1997), the refractive index and loss modulations can be introduced with
spatial phase and loss masks (see Fig. 1(b)). These cavities have emerged as
attractive options for generating frequency combs Cole et al. (2018); Wildi et
al. (2023). Despite obvious differences with rings, the nonlinear stationary
solutions of both systems share close similarities Cole et al. (2018).

The damped-driven NSE modelling passive Kerr dispersive cavities Chembo and
Menyuk (2013); Lugiato and Lefever (1987); Haelterman et al. (1992) can be
easily generalised to include a complex potential:

$\partial_{t}\psi=-i\partial_{x}^{2}\psi-(1+\delta
i)\psi+2i\left|\psi\right|^{2}\psi+V{\left(x\right)}\psi+h,\\\
V{\left(x\right)}=me^{-i\phi}\cos{\left(\frac{2\pi px}{L}\right)},$ (1)

see Tusnin et al. (2020); Englebert et al. (2023); Sun et al. (2023); Dutt et
al. (2020) for the derivation of the normalised model and estimation of
physical parameters. The cavity-laser detuning $\delta$ is our control
parameter and the equation includes a third order Kerr nonlinearity, an
external energy injection $h$ and a complex potential $V(x)$. Here, $m$
represents the depth of modulation, $p$ is the number of periods of length $L$
along the cavity path and $\phi$ represents the ratio between real and
imaginary parts. In the current study, we fix the number of oscillations to
$7$ without any loss of generality. For details on numerical methods see, e.g.
Hult (2007); Milián et al. (2018b).

Figure 1: Schematic representation of possible experimental systems holding
hybridised LSs formation. (a) Kerr optical ring cavity with integrated
electro-optical modulators for phase and amplitude modulation. (b) Thin Fabry-
Perot cavity with $\chi^{(3)}$ medium with loss and phase masks. Figure 2:
Bifurcation diagrams of HSs and PSs and corresponding spatial profiles.
(a)-(d) Norm versus detuning, for increasing values of $m$, showing the HSs
(black) and PSs (blue and turquoise).The thick/thin curves indicate
stable/unstable solutions, grey shaded area indicates HS bistability, and the
HS (grey) for $m=0$ is depicted as a reference. Insets $i)$ and $ii)$ zoom in
the overlap of the unmodulated PS branch (violet), the closed loop (turquoise)
evolving from the middle branch, and the stabilised PS (blue); $ii)$ also
shows the SN1, and the PS bifurcation point for the unmodulated system. Inset
$iii)$ shows the detachment of solutions from SN2. Inset $iv)$ highlights the
different families of stable PS, labelled as $P_{n}$, shown on the right
panels along with solitons $S_{n}$. The bottom right panel illustrates even
solitons, $S^{0-VI}$, with widths surpassing a single period $L$; periods are
shaded in grey and white. $h=1.7$.

Regarding the formation of LSs, for the homogeneous bistability case (normal
GVD), it is well-established that fronts or connections between two HS
typically exhibit one exponential decaying oscillatory tail, and one
monotonous connection related to the complex and real eigenvalues associated
with each stable state Ma et al. (2010); Parra-Rivas et al. (2016). Such an
oscillatory tail serves as a potential enabling two fronts to pin to each
other forming a LS. When two fronts are close, the large oscillation amplitude
leads to a strong effective potential generating a large coupling interval as
a function of the control parameter. Conversely, as the distance between the
fronts increases, this interval narrows. Ultimately, the LS can only form when
the fronts have null relative velocity, a condition that only occurs for
particular values of the parameters (at the Maxwell-Point). The above
ingredients induce a collapsed configuration in the bifurcation structure
known as collapsed snaking Ma et al. (2010). For the bistability of a PS and a
HS (anomalous GVD), the pinning of the two fronts, heteroclinic connections
between both states, does not depend on the number of oscillations of the LS
close to the PS but on the PS itself. Therefore, the generated branches
associated with different LSs span the same interval of the control parameter,
leading to what is known as standard homoclinic snaking. The invariance of the
system by a space-reversal symmetry, $x\rightarrow-x$, splits the LSs into two
disjointed sets with even and odd number of oscillations of the pattern,
respectively Avitabile et al. (2010); Burke and Knobloch (2007); Tlidi et al.
(1994).In general, the two snaking diagrams are found in different regimes due
to the distinct involved bistable states. These building mechanisms do not
coexist, and thus collapsed and standard snake branches are not connected.

Therefore, a natural way to observe the hybridisation of the system is by
inspecting the bifurcation diagram of these solutions. The study of the curve
connected to the homogeneous solutions is summarised in Fig. 2. Importantly,
in normal GVD regime, the stability of the upper homogeneous state is lost by
a Saddle Node (SN) bifurcation. This is opposed to the anomalous GVD where the
stability of the HS is lost at a TI point, arising due to the parametric four-
photon mixing. The introduction of the periodic potential takes the role of
Turing instability to modulate the system HSs. This induces the system to fold
the HS curve, and connect the upper and lower HS with periodic solutions that
become stable, see Fig. 2(a). The inset $i)$ zooms in the branch of the
unstable PS for $m=0$, which emerges from the middle HS, see inset $ii)$. With
the potential, this branch splits into two different solutions: one of them
being stable PS -now connected to the HSs-, and the other unstable - and now
part of a closed loop, that derives from the modulation of the middle branch-.
Note that the branches of the two (stable/unstable) modulated and the
unmodulated patterns are almost overlapping. Added to this and shown on the
insets $ii)$ and $iii)$, the middle HS becomes detached from SN1,2 and forms
an isola or closed loop.

Since the PSs can be understood as a train of solitons, accordingly the
branches follow the expected collapsed organisation until the width of a
single soliton reaches the width of the period of the modulation. This
organises the discrete periodic structures, $P_{n}$, where the label $n$
refers to the position in the collapsed snaking. A crucial effect of the
introduction of the potential is the fact that the snaking structure gains a
slant. This tilt is a finite-size effect, induced by the extra boundary
introduced by the potential at each period Kozyreff et al. (2009). Such
tilting grows with $m$, eventually exceeding the HS bistable region, see Fig.
2(a)-(d). The non-Hermiticity of the potential favours this tilt, see the
Supplemental material. This brings to light the hybridisation and dual
character of the system. Parts of the periodic solutions can no longer be
understood as a product of the HS bistability but achieve the pattern
character. The different families of solutions are displayed in the panels of
Fig. 2 and inset $iv)$, where $S_{n}$ refers to a single soliton from a
particular $P_{n}$ solution. Whenever the width of the soliton exceeds one
period of the modulation, we denote the soliton as $S^{Q}$, being $Q$ the
number of full oscillations of the upper HS. Importantly, the modulated upper
solution now takes the role of a pattern in the formation of LS. Panel $(6)$
in Fig. 2 shows solutions with an even number of oscillations, $S^{VI-0}$.

Figure 3: Bifurcation structure of LS in different hybridised scenarios. Norm
versus detuning diagram of solitonic solutions (orange). (a) Family of even
solitons for $m=0.1$ holding a tilted collapsed snaking followed by a standard
homoclinic bifurcation structure. The inset zooms in the narrower soliton
states $S_{n}$. (b) Family of odd solitons for $m=0.1$ with an analogous
structure but with a collapsed snaking structure of two solitons molecule
states of the type $\\{S_{n},S_{n}\\}$, which profiles are provided on the
panels. Lower-left inset: magnification of the SN2 of the lower homogeneous
solution. (c) Splitting of the even soliton solutions branch. Molecules of the
type $\\{S^{IV},S_{n},S_{m},S_{m},S_{n}\\}$ and $\\{S^{II},S_{n},S_{n}\\}$,
with field profiles provided on the right-hand panels. Inset: zoom in on the
collapsed snaking structure of narrow solitons. In all plots, the thick/thin
lines indicate stable/unstable states. Periodic solutions of Fig. 2 are
included for completeness.

In this study, multistability of two background solutions and patterns is
achieved by the introduction of a modulation. Note that while such background
solutions are not strictly homogeneous, the effect of the potential introduces
relatively weak modulations, see the dotted curves on panels $(1)-(6)$ in
Fig.2. With the potential, a new hybrid scheme where the basal solutions
acquire a dual character is generated and summarised in Fig. 3. The figure
captures the complexity and richness of the modulated system, where the
typical collapsed and standard homoclinic snaking link. We differentiate two
situations depending whether the PS existence lies only within or exceeds the
HS bistable area. In the first case, depicted in Figs. 3(a) and 3(b) for
$m=0.1$, the LS branches containing the different families, bifurcates from
the two main SNs. In Fig. 3(a) we show the branches of $S_{n}$. As for
stabilised $P_{n}$, their existence is organised in a tilted collapsed
snaking, until its width exceeds a modulation period. Inspecting this
structure, one could conclude that switching waves connecting both quasi-
homogeneous states are the building blocks of these solutions. Yet, once the
$P_{n}$ existence limit is exceeded, the homogeneous upper solution changes
its role to a pattern. Broader solitons with even number of oscillations of
the upper state, $S^{0,II,IV}$, are now organised in a standard homoclinic
snaking. This transition happens due to the duality of the upper HS induced by
the potential. In turn, Fig. 3(b) provides the analogous situation for even
solitons, $S^{I,III,V}$, also organised in an intertwined homoclinic snaking.
This existence (and detachment) of the branches of even an odd solitons
highlights that the physics from the case of anomalous GVD is also present in
the system. Interestingly, this hybridisation enlarges the stability area of
bright solitons (for instance $S^{III-0}$). For no modulation, and without
higher order terms, they are found just for the parameters of the Maxwell
Point Tlidi and Gelens (2010); Parra-Rivas et al. (2017). Notice, on the lower
inset, the branch approaching the SN2, and its difficult connection due to the
symmetry of the system. At the top, we show how the odd soliton branch
connects to a collapsed snaking structure of a molecule corresponding to two
solitons in two neighbouring periods $\\{S_{n},S_{n}\\}$, see the inset and
the profiles on the corresponding panels. Solitons can now bond due to the
oscillation of the upper state, induced by the potential. Although not shown,
stable molecules with different number of solitons, not necessarily laying in
neighbouring periods exist and are stable.

Figure 4: Temporal manipulation of frequency combs. Dynamical change of width,
$w$, (black curves) of the propagating field and change of the detuning (red
curves), for a) deterministic transitions between different families of
Solitons $S_{n}$; b) turning on the system and switching between different
patterns $P_{n}$. Shaded areas correspond to stable areas of the different
states. Panels $\left(i\right)$-$\left(v\right)$ depict stable switchable comb
spectra associated to solitons $S_{1-5}$.

An example of a case where the PS exceeds the bistable area is illustrated in
Fig. 3(c). As commented earlier, this introduces unequivocally the pattern
character to these periodic solutions. In this case, for solitons with even
number of oscillations of the upper state, it is interesting to note two
different organisations of the bifurcation structures. The previous branch
splits. The collapsed snaking structure of $S_{n}$ remains, although
connecting the upper HS to a PS, see inset. In turn, for broader solitons we
can observe a reminiscence of the standard homoclinic snaking surpassing the
HS bistable area, see $S_{0}$. Moreover, multiple branches of molecule-type
solution develop within the snaking. As an example, these solutions associated
with $S_{0}$ are of two types. First, solutions of the form
$\\{S^{II},S_{n},S_{n}\\}$, containing two narrow solitons $S_{n}$ within a
wide soliton $S^{II}$. Second, we find solutions
$\\{S^{IV},S_{n},S_{m},S_{m},S_{n}\\}$, containing four solitons. Due to the
symmetry of the system, we can have two pairs of different solitons embedded
in the upper state of $S_{IV}$. As expected, cases with $n\neq m$ exist in
areas of bistability of $P_{n}$ and $P_{m}$. In this figure, stable solutions
are only partially shown to illustrate the complexity and richness of the
hybrid system.

The search for real-time manipulation of frequency combs is a crucial property
that is being explored in different situations He et al. (2019); Nazemosadat
et al. (2021); Talla Mbé et al. (2017); Ivars et al. (2021, 2023b). In this
work, we uncover a mechanism for a deterministic and reversible switch between
states with different widths. Such control relies on the tilt of snaking
structures of different solutions. This phenomenon is exemplified for solitons
$S_{n}$ as provided in Fig. 4(a). The depicted numerical simulation is
performed assuming the stable soliton $S_{1}$ as the initial condition. The
detuning $\delta$ is adiabatically increased in time, up to the value of
stability of $S_{5}$. Tracking the width of the LS, it is possible to observe
the dynamic transition between them. Such transition is rapidly triggered when
the different solitons reach their existence limit. Indeed, it is the shift of
their existence intervals that allows for the reversible switching amongst
different solitonic families. Panels $i)-v)$ in Fig. 4 show the frequency
combs associated to solitons $S_{1-5}$.

Moreover, for normal GVD cavities, access to the homogeneous bistable area is
essential to start the generation of frequency combs. Recently, excitation of
dark solitons via self-injection-locking has been demonstrated Wang et al.
(2022); Jin et al. (2021). In the present work, we take advantage of the
hybridisation of the system which displaces the existence region of some
states outside the inaccessible uniform bistable area. The recession of the
first SN1 introduces an access point to the homogeneous bistable area which
can be excited by an adiabatic blue shift of the laser, as depicted in Fig.
4(b). The upper homogeneous solution is expected to evolve towards a $P_{n}$
rather than to the lower HS since the potential has an effect as the Turing
instability. Notably, patterns exhibit a deterministic switching mechanism
analogous to that observed in solitons.

To conclude, we here present a novel scenario for the formation of localised
structures. While it has been long-established that the building blocks for
solitons in non-linear systems are states that are either patterns or flat
states, we introduce a most unusual scenario in which the basal states
(patterns / flat) acquire a dual behaviour so that the system may use one
state as ’quasi-flat’ or as a pattern. Throughout the paper, such duality is
discussed in the context of a driven dissipative cavity with normal GVD,
described by the damped-driven NSE, under the introduction of a periodic non-
Hermitian modulation. The reported results represent a general new paradigm to
understand the formation of nonlinear localised structures, where physics from
normal and anomalous GVD regimes blend. In particular, this unconventional
hybridised scenario gives birth to an extraordinarily rich landscape where
many different species of frequency combs coexist, being stable and
accessible. Furthermore, the hybridisation of the LS formation scheme has
important practical consequences for real-time reshaping of combs. Beyond the
fundamental reported findings, these results pave the route for flexible
frequency comb generation.

###### Acknowledgements.

This work was supported by the Spanish government via PID2022-138321NB-C21 and
Generalitat de Catalunya via 2021 SGR 00606; C.M. acknowledges
PID2021-124618NB-C21, funded by MCIN/ AEI/10.13039/501100011033 and ‘ERDF: a
way of mak- ing Europe’ of the European Union, and Generalitat Valenciana via
PROMETEO/2021/082.

## References

  * Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, Reviews of modern physics 65, 851 (1993).
  * Arecchi et al. (1999) F. T. Arecchi, S. Boccaletti, and P. Ramazza, Physics Reports 318, 1 (1999).
  * Staliunas and Sanchez-Morcillo (2003) K. Staliunas and V. J. Sanchez-Morcillo, _Transverse patterns in nonlinear optical resonators_ , vol. 183 (Springer Science & Business Media, 2003).
  * Lugiato et al. (2015) L. Lugiato, F. Prati, and M. Brambilla, _Nonlinear optical systems_ (Cambridge University Press, 2015).
  * Turing (1990) A. M. Turing, Bulletin of mathematical biology 52, 153 (1990).
  * Knobloch (2015) E. Knobloch, conmatphys 6, 325 (2015).
  * Schäpers et al. (2000) B. Schäpers, M. Feldmann, T. Ackemann, and W. Lange, Physical Review Letters 85, 748 (2000).
  * Meron (1992) E. Meron, Physics reports 218, 1 (1992).
  * Pomeau (1986) Y. Pomeau, Physica D: Nonlinear Phenomena 23, 3 (1986).
  * Dumeige and Féron (2011) Y. Dumeige and P. Féron, Physical Review A 84, 043847 (2011).
  * Anderson et al. (2017) M. Anderson, Y. Wang, F. Leo, S. Coen, M. Erkintalo, and S. G. Murdoch, Physical Review X 7, 031031 (2017).
  * Milián et al. (2018a) C. Milián, Y. V. Kartashov, D. V. Skryabin, and L. Torner, Optics letters 43, 979 (2018a).
  * Tlidi et al. (1993) M. Tlidi, M. Georgiou, and P. Mandel, Physical Review A 48, 4605 (1993).
  * Spinelli et al. (1998) L. Spinelli, G. Tissoni, M. Brambilla, F. Prati, and L. Lugiato, Physical Review A 58, 2542 (1998).
  * Zelnik et al. (2018) Y. R. Zelnik, P. Gandhi, E. Knobloch, and E. Meron, Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (2018).
  * Al Saadi and Champneys (2021) F. Al Saadi and A. Champneys, Philosophical Transactions of the Royal Society A 379, 20200277 (2021).
  * Parra-Rivas et al. (2022) P. Parra-Rivas, A. Champneys, F. Al-Sahadi, D. Gomila, and E. Knobloch, arXiv preprint arXiv:2208.04009 (2022).
  * Kaup and Newell (1978) D. J. Kaup and A. C. Newell, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 361, 413 (1978).
  * Barashenkov and Smirnov (1996) I. Barashenkov and Y. S. Smirnov, Physical Review E 54, 5707 (1996).
  * Chembo and Menyuk (2013) Y. K. Chembo and C. R. Menyuk, Physical Review A 87, 053852 (2013).
  * Lugiato and Lefever (1987) L. A. Lugiato and R. Lefever, Physical review letters 58, 2209 (1987).
  * Haelterman et al. (1992) M. Haelterman, S. Trillo, and S. Wabnitz, Optics communications 91, 401 (1992).
  * Bao et al. (2020) H. Bao, L. Olivieri, M. Rowley, S. T. Chu, B. E. Little, R. Morandotti, D. J. Moss, J. S. T. Gongora, M. Peccianti, and A. Pasquazi, Physical Review Research 2, 023395 (2020).
  * Brasch et al. (2016) V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, Science 351, 357 (2016).
  * Kippenberg et al. (2011) T. J. Kippenberg, R. Holzwarth, and S. A. Diddams, science 332, 555 (2011).
  * Matsko et al. (2012) A. Matsko, A. Savchenkov, and L. Maleki, Optics letters 37, 43 (2012).
  * Copie et al. (2016) F. Copie, M. Conforti, A. Kudlinski, A. Mussot, and S. Trillo, Physical Review Letters 116, 143901 (2016).
  * Todd et al. (2023) C. Todd, Z. Li, S. Coen, S. G. Murdoch, G.-L. Oppo, and M. Erkintalo, Physical Review A 107, 013506 (2023).
  * Hendry et al. (2018) I. Hendry, W. Chen, Y. Wang, B. Garbin, J. Javaloyes, G.-L. Oppo, S. Coen, S. G. Murdoch, and M. Erkintalo, Physical Review A 97, 053834 (2018).
  * Tabbert et al. (2019) F. Tabbert, T. Frohoff-Hülsmann, K. Panajotov, M. Tlidi, and S. V. Gurevich, Physical Review A 100, 013818 (2019).
  * Tusnin et al. (2020) A. K. Tusnin, A. M. Tikan, and T. J. Kippenberg, Physical Review A 102, 023518 (2020).
  * Englebert et al. (2023) N. Englebert, N. Goldman, M. Erkintalo, N. Mostaan, S.-P. Gorza, F. Leo, and J. Fatome, Nature Physics pp. 1–8 (2023).
  * Dutt et al. (2020) A. Dutt, Q. Lin, L. Yuan, M. Minkov, M. Xiao, and S. Fan, Science 367, 59 (2020).
  * Sun et al. (2023) Y. Sun, P. Parra-Rivas, C. Milián, Y. V. Kartashov, M. Ferraro, F. Mangini, R. Jauberteau, F. R. Talenti, and S. Wabnitz, Physical Review Letters 131, 137201 (2023).
  * Ivars et al. (2022) S. B. Ivars, M. Botey, R. Herrero, and K. Staliunas, Chaos, Solitons & Fractals 165, 112774 (2022).
  * Akhter et al. (2023) M. N. Akhter, S. B. Ivars, M. Botey, R. Herrero, and K. Staliunas, Physical Review Letters 131, 043604 (2023).
  * Ivars et al. (2023a) S. B. Ivars, M. Botey, R. Herrero, and K. Staliunas, Chaos, Solitons & Fractals 168, 113089 (2023a).
  * Ahmed et al. (2018) W. W. Ahmed, S. Kumar, J. Medina, M. Botey, R. Herrero, and K. Staliunas, Optics letters 43, 2511 (2018).
  * Makris et al. (2015) K. G. Makris, Z. H. Musslimani, D. N. Christodoulides, and S. Rotter, Nature communications 6, 7257 (2015).
  * Barashenkov et al. (2016) I. Barashenkov, D. Zezyulin, and V. Konotop, New Journal of Physics 18, 075015 (2016).
  * Feng et al. (2014) L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, Science 346, 972 (2014).
  * Ashida et al. (2020) Y. Ashida, Z. Gong, and M. Ueda, Advances in Physics 69, 249 (2020).
  * Konotop et al. (2016) V. V. Konotop, J. Yang, and D. A. Zezyulin, Reviews of Modern Physics 88, 035002 (2016).
  * Suchkov et al. (2016) S. V. Suchkov, A. A. Sukhorukov, J. Huang, S. V. Dmitriev, C. Lee, and Y. S. Kivshar, Laser & Photonics Reviews 10, 177 (2016).
  * Özdemir et al. (2019) Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nature materials 18, 783 (2019).
  * Cheng et al. (2023) D. Cheng, E. Lustig, K. Wang, and S. Fan, Light: Science & Applications 12, 158 (2023).
  * Hu et al. (2024) Y. Hu, D. Zhu, S. Lu, X. Zhu, Y. Song, D. Renaud, D. Assumpcao, R. Cheng, C. Xin, M. Yeh, et al., arXiv preprint arXiv:2404.06398 (2024).
  * Lugiato and Prati (2010) L. Lugiato and F. Prati, Physical review letters 104, 233902 (2010).
  * Hachair et al. (2004) X. Hachair, S. Barland, L. Furfaro, M. Giudici, S. Balle, J. R. Tredicce, M. Brambilla, T. Maggipinto, I. M. Perrini, G. Tissoni, et al., Physical Review A 69, 043817 (2004).
  * Brambilla et al. (1996) M. Brambilla, L. Lugiato, and M. Stefani, Europhysics Letters 34, 109 (1996).
  * Peschel et al. (2003) U. Peschel, D. Michaelis, and C. O. Weiss, IEEE journal of quantum electronics 39, 51 (2003).
  * Brambilla et al. (1997) M. Brambilla, L. Lugiato, F. Prati, L. Spinelli, and W. Firth, Physical review letters 79, 2042 (1997).
  * Cole et al. (2018) D. C. Cole, A. Gatti, S. B. Papp, F. Prati, and L. Lugiato, Physical Review A 98, 013831 (2018).
  * Wildi et al. (2023) T. Wildi, M. A. Gaafar, T. Voumard, M. Ludwig, and T. Herr, Optica 10, 650 (2023).
  * Hult (2007) J. Hult, Journal of Lightwave Technology 25, 3770 (2007).
  * Milián et al. (2018b) C. Milián, Y. V. Kartashov, D. V. Skryabin, and L. Torner, Physical review letters 121, 103903 (2018b).
  * Ma et al. (2010) Y.-P. Ma, J. Burke, and E. Knobloch, Physica D: Nonlinear Phenomena 239, 1867 (2010).
  * Parra-Rivas et al. (2016) P. Parra-Rivas, E. Knobloch, D. Gomila, and L. Gelens, Physical Review A 93, 063839 (2016).
  * Avitabile et al. (2010) D. Avitabile, D. J. Lloyd, J. Burke, E. Knobloch, and B. Sandstede, SIAM Journal on Applied Dynamical Systems 9, 704 (2010).
  * Burke and Knobloch (2007) J. Burke and E. Knobloch, Chaos: An Interdisciplinary Journal of Nonlinear Science 17 (2007).
  * Tlidi et al. (1994) M. Tlidi, P. Mandel, and R. Lefever, Physical review letters 73, 640 (1994).
  * Kozyreff et al. (2009) G. Kozyreff, P. Assemat, and S. J. Chapman, Physical review letters 103, 164501 (2009).
  * Tlidi and Gelens (2010) M. Tlidi and L. Gelens, Optics letters 35, 306 (2010).
  * Parra-Rivas et al. (2017) P. Parra-Rivas, D. Gomila, and L. Gelens, Physical Review A 95, 053863 (2017).
  * He et al. (2019) Y. He, Q.-F. Yang, J. Ling, R. Luo, H. Liang, M. Li, B. Shen, H. Wang, K. Vahala, and Q. Lin, Optica 6, 1138 (2019).
  * Nazemosadat et al. (2021) E. Nazemosadat, A. Fülöp, Ó. B. Helgason, P.-H. Wang, Y. Xuan, D. E. Leaird, M. Qi, E. Silvestre, A. M. Weiner, et al., Physical Review A 103, 013513 (2021).
  * Talla Mbé et al. (2017) J. H. Talla Mbé, C. Milián, and Y. K. Chembo, The European Physical Journal D 71, 1 (2017).
  * Ivars et al. (2021) S. B. Ivars, Y. V. Kartashov, L. Torner, J. A. Conejero, and C. Milián, Physical review letters 126, 063903 (2021).
  * Ivars et al. (2023b) S. B. Ivars, Y. V. Kartashov, P. F. de Córdoba, J. A. Conejero, L. Torner, and C. Milián, Nature Photonics pp. 1–8 (2023b).
  * Wang et al. (2022) H. Wang, B. Shen, Y. Yu, Z. Yuan, C. Bao, W. Jin, L. Chang, M. A. Leal, A. Feshali, M. Paniccia, et al., Physical Review A 106, 053508 (2022).
  * Jin et al. (2021) W. Jin, Q.-F. Yang, L. Chang, B. Shen, H. Wang, M. A. Leal, L. Wu, M. Gao, A. Feshali, M. Paniccia, et al., Nature Photonics 15, 346 (2021).


# Closely spaced co-rotating helical vortices: General solutions

A. Castillo-Castellanos  Aix Marseille Université, CNRS, Centrale Marseille,
IRPHE, Marseille, France<EMAIL_ADDRESS>E.
Durán Venegas  Aix Marseille Université, CNRS, Centrale Marseille, IRPHE,
Marseille, France  Departamento de Bioingeniería e Ingeniería Aeroespacial,
Universidad Carlos III de Madrid, Spain  S. Le Dizès  Aix Marseille
Université, CNRS, Centrale Marseille, IRPHE, Marseille, France

###### Abstract

In this work, we present general solutions for closely spaced co-rotating
helical vortices using a filament approach. For these vortex structures,
helical symmetry is broken, but solutions maintain a form of spatial
periodicity. We show there exists a moving frame where vortex elements move
along the structure without distorting it. These solutions can be used to
characterize the evolution of the twin-tip vortices produced in the wake of a
tip-splitting rotor blade. The resulting wake has a dual nature. Locally, the
structure behaves much like an helical pair aligned with the locally tangent
flow. However, as we move away from the vortex structure the induced flow is
reminiscent of an ‘equivalent’ helical vortex with thicker core.

###### pacs:

##  I Introduction

Rotating blades, such as those of a helicopter rotor or a horizontal-axis wind
turbine, generate concentrated vortices at their tips, transported downstream,
creating a persistent helical pattern. These concentrated vortices can have a
core size as small as 1% of the rotor radius [1]. The structure and stability
of these helical vortices are associated with several practical issues
actively investigated today. One of these issues concerns the interaction
between a tip-vortex and a downstream surface (Vortex-Surface Interactions),
which causes significant noise and vibration problems. In particular, the
interaction between a tip-vortex and a rotor blade (Blade-Vortex Interaction,
BVI) is responsible for reduced power output and premature structural
component fatigue in wind turbines [2, 3], and undesirable noise in helicopter
rotors during low speed and descending flight [4]. BVI mitigation methods can
be classified as active and passive [5, 6, 7]. Active methods include
dynamical modifications of the wake, for instance, by steering the wake
towards a particular direction [8, 9], or by introducing a periodic
perturbation [10, 11, 12]. Passive methods use a modified blade geometry to
enhance particular flow features, for instance, by introducing serrations on
the leading edge [13, 14] or through the use of winglets [15, 16] and slotted
tips [17]. One alternative proposed by Brocklehurst and Pike [18] is to use a
modified airfoil to split the tip vortex into two closely spaced vortices. As
the tip vortices are advected downstream, they interact with each other and
eventually merge back into a large more diffuse vortex [19]. This approach
offers significant noise reduction and has no adverse effects on control load
or performance [20, 21, 22, 23], but little is known about the wake structure
and the mechanism leading to a larger vortex. One of our motivations is to
provide information on the vortical wake that could exist before the merging
process. We consider an ideal framework of vortex filaments to illustrate the
type of structures that can be created by the interaction of two closely
spaced but distinct helical vortices.

A wide array of models, including analytical, experimental, and numerical
approaches for rotor wake aerodynamics are available [24]. Direct numerical
simulations, which require solving the Navier-Stokes equations for all scales
of fluid motion up to the Kolmogorov scale remain prohibitively expensive for
realistic configurations. These numerical simulations require a turbulence
model, such as the Reynolds-averaged Navier-Stokes (RANS) or large eddy
simulation (LES), to take into account the effect of the unresolved small
scales [24, 25, 26, 27]. These simulations may take into account the rotor
directly, by solving the boundary layer along the rotor blade, or indirectly,
by representing the rotor blades as body forces as in the actuator line method
(ACL). For applications, the design and optimization of rotors often rely on
computationally inexpensive models that describe the aerodynamic loads and
wake dynamics. In the wind turbine community, blade loading is usually
obtained using the Blade Element Momentum (BEM) theory, which is based on the
actuator disc principle, corresponding to an infinitesimally thin rotor with
an infinite number of blades. This approach has been progressively extended to
include a finite number of blades, three-dimensional effects, unsteady flow
loading [28], and blade thickness [29]. But using BEM requires a wake model
and airfoil data [30, 31, 28]. Wakes with a prescribed helical symmetry are
often used [32, 33]. Such wakes do not account for the inboard motion and
subsequent radial contraction (resp. expansion) observed in helicopter flight
(resp. wind turbine) regimes. Additionally, a prescribed helical wake simply
cannot describe the tip-vortex interaction considered here.

A free vortex method is more adequate to describe the structure and dynamics
of the wake (see for instance [34]). Although the inboard motion associated
with vortex formation is not easily described, the free-vortex method captures
the contraction/expansion of the wake. This method has been used to extend the
analytical solution of a single infinite helix [35, 36] to multiple vortices
[37] and semi-infinite configurations [38, 39]. Using Joukowski’s wake model,
for which there is only one free tip-vortex per blade, Durán Venegas  _et al._
[39] were able to find steady solutions describing the wake in all the
vertical flight regimes of a helicopter (hover, climb, descent) and all wind
turbine regimes, including cases for which the vortex structure is strongly
deformed in the near wake and crosses the rotor plane.

For Joukowski’s model, the solution matches in the far-field a uniform helix
solution. When two (or more) free vortices are emitted by each blade, the
solution becomes more complex even in the far-field. As shown by [37] for two
counter-rotating vortices, steady solutions can still be obtained. These
solutions are no longer uniform but still exhibit a spatial periodicity
property. In the present work, we shall use the method developed in [39, 37]
to describe the wake of tip-splitting rotors.

As for counter-rotating vortices, we expect contraction/expansion of the wake
close to the rotor as well as a complex far-field. The closest canonical
example to these closely spaced co-rotating helical vortex pairs is a pair of
equistrength parallel vortices. In such a case, the parameters driving the
vortex dynamics are the circulation $\Gamma$, separation distance $d$, and
vortex core size $a$. In general, $a$ grows with time due to viscous diffusion
of vorticity with $a(t)=\sqrt{4\nu t+a_{0}^{2}}$ where $\nu$ is the kinematic
viscosity, while $d$ remains virtually constant. Simultaneously, the vortex
pair rotates around the total vortex centroid with a constant rotation rate
$\omega=\Gamma/(\pi d^{2})$ [40]. When the ratio $a/d$ exceeds a threshold
value around 0.23, $d$ suddenly decreases and the vortices merge into a single
entity (see, for instance [41, 42, 43, 44, 45]). In the present work, we
neglect viscous diffusion and consider a sufficiently small value of $a/d$,
typically 0.1 such that the vortices remain distinct and the filament approach
justified. We shall see that the rotation motion of the two vortices around
each other is also present in our solutions.

The paper is organized as follows. Section II briefly introduces the vortex
filament framework used. Section III is concerned with the far-field. We show
how steady spatially-periodic solutions are obtained from a few geometric
parameters. In section IV, these solutions are then used to compute the far-
field contribution for the near-field solutions obtained close to a tip-
splitting rotor. We present solutions representative of both wind turbine and
helicopter regimes and discuss the rotor influence on the wake structure. We
also provide a simple model explaining the observed spatial variations of the
parameters.

##  II Framework

Our problem is defined by three spatial scales. The radius $R$ of the rotor,
the separation distance $d$, and the vortex core size $a$. The filament
approach is perfectly justified as long as the core size remains small
compared to other spatial scales (separation distances, local curvature
radius). Closely spaced filament vortices correspond to situations for which
$a\ll d\ll R$. We use the same framework as in [37, 39]. All the vorticity is
considered as being concentrated along lines which move as material lines in
the fluid according to

$\frac{\mathrm{d}\mathbf{X}_{i}}{\mathrm{d}t}=\mathbf{U}(\mathbf{X}_{i})+\mathbf{U}^{\infty}$
(1)

where $\mathbf{X}_{i}$ is the position vector of the $i$-th vortex filament,
$\mathbf{U}$ the velocity induced by the vortices and $\mathbf{U}^{\infty}$ an
external velocity field. The induced velocity $\mathbf{U}$ is given by the
Biot-Savart law

$\mathbf{U}(\mathbf{X}_{i})=\sum_{j=1}^{n}\frac{\Gamma_{j}}{4\pi}\int\frac{(\mathbf{X}_{j}-\mathbf{X}_{i})\times\mathrm{d}\mathbf{T}_{j}}{|(\mathbf{X}_{j}-\mathbf{X}_{i})|^{2}}$
(2)

where the integrals cover each vortex filament defined by its circulation
$\Gamma_{j}$, position vector $\mathbf{X}_{j}$ and corresponding tangent unit
vector $\mathbf{T}_{j}$. Filaments are discretised in straight segments
$[\mathbf{X}_{i}^{n},\mathbf{X}_{i}^{n+1}]$ in order to compute the velocity
field and follow its evolution. The divergence in (2) is treated using the
cut-off method with a Gaussian vorticity profile [46]: to determine the local
contribution to the velocity field at $\mathbf{X}_{i}^{n}$ from the
neighboring segments $[\mathbf{X}_{i}^{n-1},\mathbf{X}_{i}^{n}]$ and
$[\mathbf{X}_{i}^{n},\mathbf{X}_{i}^{n+1}]$, we replace the two segments by an
arc of circle passing through
$[\mathbf{X}_{i}^{n-1},\mathbf{X}_{i}^{n},\mathbf{X}_{i}^{n+1}]$ and use the
cut-off formula. The discretised expressions that we shall use for the induced
velocity are given in [37]. When local contributions are not taken into
account, the number of segments required to obtain a good approximation of the
induced velocity increases by a factor $O(2\pi\rho/a)$, where $\rho$ is the
local curvature radius [37].

Our objective is to find a vortex structure that is steady in the frame of the
rotor. This condition of steadyness can be written as

$\left(\mathbf{U}(\mathbf{X}_{j})-\mathbf{U}^{F}\right)\times\mathbf{T}_{j}=0$
(3)

where $\mathbf{U}^{F}$ is the frame velocity, which simply indicates that
vortices are moving along the vortex structure. In practice, we shall write
the steadiness condition as a system of ordinary equations. In a cylindrical
frame, it can be written for each vortex as

$\frac{\mathrm{d}r_{j}}{\mathrm{d}\zeta}=U_{r}(\mathbf{X}_{j}),\quad\quad\frac{\mathrm{d}\theta_{j}}{\mathrm{d}\zeta}=\Omega(\mathbf{X}_{j})-\Omega^{F},\quad\quad\frac{\mathrm{d}z_{j}}{\mathrm{d}\zeta}=U_{z}(\mathbf{X}_{j})-U_{z}^{F},$
(4)

where $(r_{j}(\zeta),\theta_{j}(\zeta),z_{j}(\zeta))$ are the radial, angular,
and axial coordinates of the $j$-th vortex, $(U_{r},r_{j}\Omega,U_{z})$ the
corresponding velocity components, while $\Omega^{F}$ and $U_{z}^{F}$ are the
angular and axial velocity of the rotor frame. Each vortex curve is
parametrized by $\zeta$. These systems are to be solved with boundary
conditions on the rotor, at $\zeta=0$, where the position of each vortex is
prescribed, and far-field boundary conditions at $\zeta\rightarrow\infty$.

These last boundary conditions are not trivial as the far-field is a priori
unknown. Therefore, our first task is to characterize the far-field.

##  III Far-field solutions

We are interested in steady solutions created by the emission of two closely-
spaced co-rotating vortices from a rotating blade tip under an external axial
flow. Based on the observation of the tip vortices [47], one may naturally
expect close to the rotor, a global contraction of the structure in a climbing
helicopter regime, and a global expansion in a wind turbine regime. In the
far-field, a quasi-uniform regime is expected. However, this regime is not as
simple as for a single-tip vortex. We shall see that while not perfectly
uniform, the vortices deform but exhibit a certain periodicity induced by
their mutual interaction. The description of this far-field can be analyzed
using the method introduced in [37] for counter-rotating vortices. These
solutions are defined only by the geometrical parameters introduced in the
following section.

###  III.1 Simplified approach

####  III.1.1 Geometrical parameters

(a) (b)

Figure 1: Geometric parameters of closely spaced tip vortices: separation
distance $d$, radius $R$, vortex core size $a$, axial pitches $H$ and $h$, and
twist parameter $\beta=\pm 2\pi/\tilde{\theta}=H/h$.

(a) (b) (c)

Figure 2: Geometry described by (6) for $H=10$, $R=10$, $d=1$ and: (a)
$h_{\tau}=128$ (or $\beta\approx 0.5$); (b) $h_{\tau}=32$ (or $\beta\approx
2$); and (c) $h_{\tau}=2$ (or $\beta\approx 31.8$).

In the present situation, the far-field structure is expected to be close to a
double-helix inscribed on a “larger” underlying helix. Such idealized
structure is defined by the geometric parameters identified in figure 1. The
double-helix has a radius $d/2$ and pitch $h_{\tau}$ (corresponding to an
axial pitch $h$) while the large helix, noted $\mathscr{H}$, has a radius $R$
and axial pitch $H$. Another important parameter is the ratio $\beta\equiv
H/h=2\pi/\tilde{\theta}$ (see figure 1). Depending on the values of $H$, $R$,
$h_{\tau}$, and $d$, the double-helix structure may describe: (i) a leapfrog-
type pattern, where vortices trade places every $1/\beta$ turns (figure 2a);
(ii) a relatively sparse braid (figure 2b) ; or (iii) a dense ‘telephone
cord’-type pattern (figure 2c). We are typically in situation (i) when
$\beta<1$ and $h_{\tau}/d>10$, and in situation (iii) when $h_{\tau}/d<5$ and
$\beta>5$.

The curve $\mathbf{X}_{0}$ described by $\mathscr{H}$ is defined in a
Cartesian frame as a function of the angular coordinate $\theta_{0}$ by

$\displaystyle x_{0}$ $\displaystyle=$ $\displaystyle R\cos\theta_{0}$ (5a)
$\displaystyle y_{0}$ $\displaystyle=$ $\displaystyle R\sin\theta_{0}$ (5b)
$\displaystyle z_{0}$ $\displaystyle=$ $\displaystyle H\theta_{0}/(2\pi)$ (5c)

while the curves $(\mathbf{X}_{i}$, $i=1,2)$ described by the double-helix are
defined by

$\displaystyle x_{i}$ $\displaystyle=$ $\displaystyle
x_{0}+(d/2)\cos\phi_{i}\cos\theta_{0}-c_{\tau}(d/2)\sin\phi_{i}\sin\theta_{0}$
(6a) $\displaystyle y_{i}$ $\displaystyle=$ $\displaystyle
y_{0}+(d/2)\cos\phi_{i}\sin\theta_{0}+c_{\tau}(d/2)\sin\phi_{i}\cos\theta_{0}$
(6b) $\displaystyle z_{i}$ $\displaystyle=$ $\displaystyle
z_{0}-c_{\kappa}(d/2)\sin\phi_{i}$ (6c)

where $\phi_{1}\equiv\beta\theta_{0}$ and $\phi_{2}\equiv\beta\theta_{0}+\pi$
define the orientation of the helices relative to the chord plane. The torsion
coefficient $c_{\tau}$ and the curvature coefficient $c_{\kappa}$ of
$\mathscr{H}$ are given by,

$c_{\tau}\equiv\frac{|H|}{\sqrt{4\pi^{2}R^{2}+H^{2}}}~{},\quad\quad\quad\quad
c_{\kappa}\equiv\frac{2\pi R}{\sqrt{4\pi^{2}R^{2}+H^{2}}}~{}.$ (7)

We also note that $h_{\tau}$ and $h$ are related with each other by the
relation

$h_{\tau}=\frac{h}{c_{\tau}}.$ (8)

The above structure satisfies the following form of spatial periodicity: it is
invariant by the double operation $z\rightarrow z+h$ and
$\theta_{0}\rightarrow\theta_{0}+\tilde{\theta}$. Moreover, since we consider
vortices of equal core size $a$ and circulation $\Gamma$, both vortices are
deemed as interchangeable, such that our periodic domain is further reduced to
a domain of axial length $h/2$ and azimuthal angle $\tilde{\theta}/2$. In the
following, we keep this property to obtain the steady solutions. In
particular, we assume that there is a single location over an axial period
$h/2$ where both vortices are at the same azimuth. We chose this azimuth to
define the radius $R$ and separation distance $d$. The axial period $h$ then
defines a mean axial pitch while $H$ is obtained from $H=2\pi
h/\tilde{\theta}$. In addition to the core size $a$ that is assumed constant,
these length scales can be used to define four dimensionless parameters

$R^{*}\equiv\frac{R}{d},\quad\quad H^{*}\equiv\frac{H}{d},\quad\quad
h_{\tau}^{*}\equiv\frac{h_{\tau}}{d},\quad\quad\varepsilon\equiv\frac{a}{d},$
(9)

that will characterize our solutions in the far-field.

####  III.1.2 Characteristics of the moving frame

A vortex ring and a helical vortex are examples of vortex structures that move
in space at a constant speed without deformation [48]. The double-helix
structure is also bound to move by the induction of velocity it generates. For
this structure, one can naturally assume that the self-induced velocity is
predominantly composed of:

  * •

a translation of velocity $u_{t}$ and rotation of angular velocity $\omega$
around the double-helix axis, i.e., $\mathscr{H}$.

  * •

a translation of velocity $U_{z}$ and rotation of angular velocity $\Omega$
around the large helix axis, i.e., the vertical axis.

Under this assumption, we can show there exists a unique frame rotating and
translating along the vertical axis where the double-helix geometry shown in
figure 1 remains steady. The idea is to use the property that a helix of pitch
$H$ is unperturbed by a rotation of angular velocity $\Omega^{a}$ and
translation of axial velocity $U_{z}^{a}$ provided that

$U_{z}^{a}/\Omega^{a}=\pm H/(2\pi),$ (10)

where the sign is $+$ for right-handed helices, and $-$ for left-handed
helices. The self-rotation of the double-helix can then be cancelled by adding
a motion of velocity $u_{t}^{a}=\mp h_{\tau}\omega/(2\pi)$ along the double-
helix axis. This motion corresponds to an axial rotation of angular velocity
$c_{\kappa}u_{t}^{a}/R$ plus an axial translation of velocity
$c_{\tau}u_{t}^{a}$. These velocities add up to the self-induced velocities of
the large helix. Both sums can be cancelled if one chooses the frame
velocities such that

$\Omega^{F}=\Omega+\frac{c_{\kappa}u_{t}}{R}\mp\frac{h}{H}\omega~{},\quad\quad\quad\quad
U_{z}^{F}=U_{z}+c_{\tau}u_{t}\mp\frac{c_{\tau}h_{\tau}}{2\pi}\omega~{}.$ (11)

These expressions will be used as initial guess values for the frame velocity
in the numerical procedure described in the next section. If the vortices have
a circulation $\Gamma$, we shall use

  * •

for the self-induced motion ($\omega$ and $u_{t}$) of the double-helix, the
values obtained by the cut-off method with a Gaussian vortex core of radius
$a$ for a straight double-helix of circulation $\Gamma$ and the same
geometrical parameters (pitch $h_{\tau}$, radius $d/2$).

  * •

for the self-induced motion ($\Omega$ and $U_{z}$) of the large helix
$\mathscr{H}$, the values obtained by the cut-off method with a Gaussian
vortex core of “effective” radius $d/2$ for a helix of circulation $2\Gamma$,
pitch $H$ and radius $R$.

Typical values for these velocities are provided as function of the non-
dimensional pitch in figure 3. In the following, we shall refer to (11) with
these estimates as the rectilinear approximation for the frame velocity.

(a) (b)

Figure 3:  Dimensionless rotation rate and axial speed of the vortex elements
on a helix of radius $R$, circulation $2\Gamma$ and core size $d=0.1R$ as a
function of the pitch $H/R$ (in blue), and a double-helix of radius $d/2$,
circulation $\Gamma$ and core size $a=0.1d$ as a function of the pitch
$2h_{\tau}/d$ (in orange). The dotted line in (a) indicates the value
$d^{2}\omega^{h}/\Gamma=1/\pi$ for $h_{\tau}\gg 1$.

###  III.2 Numerical solutions

In the previous section, a prescribed shape for the solutions has been
assumed. This hypothesis will turn out to be too strong. We shall see that the
helices must deform to be compatible with the self-induced motion. For the
numerical procedure, we use the separation distance $d$ and the vortex
circulation $\Gamma$ to normalize time and length scales. The solutions will
then only depend on the four dimensionless parameters $R^{*}$, $H^{*}$,
$h_{\tau}^{*}$ and $\varepsilon$ introduced in (9). To reduce the size of the
parameter space, only $H^{*}$ and $h_{\tau}^{*}$ will be varied. In most
cases, we shall use $R^{*}=9.5$ and $\varepsilon=0.1$, for which the filament
approach is expected to provide reliable results.

####  III.2.1 Governing equations in helical coordinates

(a) (b)

Figure 4:  (a) Local helical basis following $\mathscr{H}$ and (b) 2-D
coordinates in the $(r^{\prime},b^{\prime})$ plane.

To describe the far-field, it is useful to introduce a local Frenet-Serret
basis $(\hat{\mathbf{e}}_{r0},\hat{\mathbf{e}}_{t},\hat{\mathbf{e}}_{b})$
based on $\mathscr{H}$ and centered at $\mathbf{X}_{0}$ as illustrated in
figure 4a. The position vector $\mathbf{X}_{j}$ will then be defined by its
local coordinates $r^{\prime}_{j}=(r_{j}-R^{*})$ and $b^{\prime}_{j}$ in the
local perpendicular basis $(\hat{\mathbf{e}}_{r0},\hat{\mathbf{e}}_{b})$, see
figure 4b. Developing $\mathbf{T}_{j}$ in terms of
$(\hat{\mathbf{e}}_{r0},\hat{\mathbf{e}}_{t},\hat{\mathbf{e}}_{b})$, the
steadiness condition (3) may be written as

$\displaystyle\frac{dr^{\prime}_{j}}{d\theta_{0}}=\left(\frac{R^{*}}{c_{\kappa}}+c_{\kappa}r^{\prime}_{j}\right)\frac{U_{r_{0}}-U_{r_{0}}^{F}}{U_{t}-U_{t}^{F}}+c_{\tau}b^{\prime}_{j}\quad\quad\frac{db^{\prime}_{j}}{d\theta_{0}}=\left(\frac{R^{*}}{c_{\kappa}}+c_{\kappa}r^{\prime}_{j}\right)\frac{U_{b}-U_{b}^{F}}{U_{t}-U_{t}^{F}}-c_{\tau}r^{\prime}_{j}.$
(12)

Because of the spatial periodicity and interchangeability conditions, the
calculation domain is comprised between $\theta_{0}=0$ and
$\theta_{0}=\tilde{\theta}/2$, with

$\tilde{\theta}=2\pi/\beta$ (13)

being the half-period of the vortex pair. Our description is completed by the
following boundary conditions

$\begin{cases}r^{\prime}_{1}(0)=r^{\prime}_{2}(\tilde{\theta}/2)=+\frac{1}{2},\\\
r^{\prime}_{2}(0)=r^{\prime}_{1}(\tilde{\theta}/2)=-\frac{1}{2},\end{cases}\quad\quad\begin{cases}b^{\prime}_{1}(0)=b^{\prime}_{2}(\tilde{\theta}/2)=0,\\\
b^{\prime}_{2}(0)=b^{\prime}_{1}(\tilde{\theta}/2)=0.\end{cases}$ (14)

Equation (12) is solved numerically as a non-linear minimization problem.
Starting from an initial guess for $(r^{\prime}_{j},b^{\prime}_{j})$ obtained
from (6)

$r^{\prime}_{1}=\frac{1}{2}\cos(\beta\theta_{0}),\quad\quad
b^{\prime}_{1}=\frac{1}{2}\sin(\beta\theta_{0}),\quad\quad
r^{\prime}_{2}=\frac{1}{2}\cos(\beta\theta_{0}+\pi),\quad\quad
b^{\prime}_{2}=\frac{1}{2}\sin(\beta\theta_{0}+\pi),$ (15)

and an initial guess for $\Omega^{F}$ and $U_{z}^{F}$ obtained from (11), one
applies an iterative procedure. At step $n$ of this procedure, one integrates
(12) between $\theta_{0}=0$ and $\theta_{0}=\tilde{\theta}/2$ to compute the
new values of $\Omega^{F}$ and $U_{z}^{F}$ which satisfy the boundary
conditions. Then, the next iteration of $(r^{\prime}_{j},b^{\prime}_{j})$ is
computed from the residual value of (12) using Newton’s method. Steps are
repeated until the current and subsequent iterations converge. To fix the
spatial resolution, we check for grid convergence using the primitive
variables, local curvature coefficient and frame velocities. In practice, for
most calculations we use $p_{n}=48$ segments per period.

####  III.2.2 Geometry of the solutions

(a) (b)

Figure 5:  Close-up of the deformed vortex structure for $H^{*}=15$ and
$h_{\tau}^{*}=5$ over a single period. Numerical solution (in colour)
superimposed to the undeformed initial guess (in dashed lines).

In this section, we characterize the geometry of the solutions, especially
their departure from the (undeformed) initial guess. In general, solutions
keep the overall form of the approximated solutions, but the helical
symmetries are now broken (figure 5). Both the double-helix structure and the
large underlying helix deform. Vortex segments placed in the first and second
half-periods are usually drawn in opposite directions. This combined motion
affects the local curvature, as it straightens vortex filaments around the
nodes, and bends them around the anti-nodes. Some useful indicators of this
deformation are the relative changes of the separation distance and of the
mean radius defined by

$\Delta
d(\theta_{0})=\sqrt{(r^{\prime}_{1}-r^{\prime}_{2})^{2}+(b^{\prime}_{1}-b^{\prime}_{2})^{2}}-1,\quad\quad\Delta
R(\theta_{0})=\frac{\sqrt{x_{1}^{2}+y_{1}^{2}}+\sqrt{x_{2}^{2}+y_{2}^{2}}}{2R^{*}}-1,$
(16)

respectively. A positive (resp. negative) $\Delta d$ indicates that the
distance between the two vortices is increasing (resp. decreasing). Similarly,
a positive (resp. negative) $\Delta R$ indicates an increase (resp. decrease)
of the mean radius. By construction, $\Delta d(0)=\Delta
d(\tilde{\theta}/2)=\Delta R(0)=\Delta R(\tilde{\theta}/2)=0$. As seen in
figure 6a-b, both quantities are shown to vary over a single period and their
evolution varies with the type of solution. For leapfrogging wakes (i),
$\Delta d$ remains relatively constant while $\Delta R$ is positive with local
maxima near $\tilde{\theta}/4$. For sparsely braided wakes (ii), $\Delta d$
reaches a plateau while $\Delta R$ becomes close to zero. For densely coiled
wakes (iii), $\Delta R$ may take positive and negative values, while the
maximum in $\tilde{\theta}/4$ becomes a local extremum. The maximum $|\Delta
d|$ and $|\Delta R|$ initially contracts, before increasing once again as
$H^{*}$ decreases (figure 7). Below this point, the mutually-induced velocity,
which drives the rotation of the vortex pair, is no longer the dominant
contribution to the total induced velocity and numerical convergence becomes
increasingly hard to obtain for this family of solutions.

(a) (b)

(c) (d)

Figure 6:  (a,c) Change in the separation distance $\Delta d$ and (b,d) mean
radius $\Delta R$ as function of $\theta_{0}$ for: (a,b) $H^{*}=15$ and
different $h_{\tau}^{*}$; and (c,d) $h_{\tau}^{*}=10$ and different $H^{*}$.

(a) (b)

Figure 7:  Maximum change in (a) the separation distance $|\Delta d|$ and (b)
mean radius $|\Delta R|$ as function of $H^{*}$ for $h_{\tau}^{*}=2$ and $N$
vortex pairs with and without a straight hub vortex of circulation
$-2N\Gamma$.

For applications, it is interesting to compare these solutions to a similar
configuration composed of $N$ pairs of vortices of circulation $\Gamma$. In
this context, this wake geometry would represent the far-field in a Joukowski
rotor wake model [49], but produced by a $N$-bladed tip-splitting rotor.
Figure 8a (resp. 8b) displays the wake geometry obtained for two (resp. three)
vortex pairs for $h_{\tau}^{*}=2$. In general, for equal
$(H^{*},h_{\tau}^{*})$, deviations from (6) are found to increase with $N$
(figure 7). However, if we base our comparison on the separating distance
between neighbouring spires, the deformation actually decreases with $N$. A
similar effect has been reported for counter-rotating helical vortices [39].
Finally, we note that including a straight vortex hub of circulation
$-2N\Gamma$ was shown to have little effect on the wake geometry. Unless
otherwise stated, in the sequel, we shall focus on the case for $N=1$ without
a hub vortex.

(a) (b) (c)

Figure 8:  Deformed structures for (a) $N=1$, (b) $N=2$ and (b) $N=3$ vortex
pairs for $H^{*}/N=15$ and $h_{\tau}^{*}=2$.

####  III.2.3 Frame and tangential velocities

(a) (b)

Figure 9:  Frame velocities as function of $h_{\tau}^{*}$. Numerical solutions
(in solid lines) are compared to the rectilinear approximation (in dash-dotted
lines) and to a first order approximation based on the undeformed initial
guess (in dashed lines).

(a) (b)

Figure 10:  Frame velocities as function of $H^{*}$. Numerical solutions (in
solid lines) are compared to the rectilinear approximation (in dash-dotted
lines) and to a first order approximation based on the undeformed initial
guess (in dashed lines).

(a) (b)

Figure 11:  Evolution of $\overline{U}_{tan}$ (in solid lines) and
$\overline{U}_{t}$ (in dotted lines) as function of (a) $H^{*}$ and (b)
$h_{\tau}^{*}$. Numerical solutions are also compared to the velocity
associated with the moving frame $\overline{U}^{F}_{t}$ (in dashed lines).

The numerical procedure also provides the frame velocities $\Omega^{F}$ and
$U_{z}^{F}$ where the structure is steady. These quantities are plotted as a
function of $h_{\tau}^{*}$ (resp. $H^{*}$) in figure 9 (resp. figure 10). They
are compared to the rectilinear approximation obtained in §III.1.2 (in dash-
dotted lines) and to a first order approximation obtained by numerically
computing the frame velocities with an undeformed solution (in dashed lines).
Both approximations work well for large $H^{*}$ and large $h_{\tau}^{*}$,
i.e., when self-induction is expected to be small. The first order
approximation, which takes into account the mutual- and self-induced
velocities but neglects the deformation, is naturally better and provides
accurate estimates whenever the deformations are weak.

Now let us consider the velocities relative to the moving frame. By
construction, the vortex elements are advected along the stationary vortex
structure, where positive tangential velocity $U_{tan}$ indicates an advection
in the positive axial direction and vice versa. In general, values of
$U_{tan}$ are positive (resp. negative) for $h^{*}_{\tau}<0$ (resp.
$h^{*}_{\tau}>0$). The tangential velocity varies with respect to $\theta_{0}$
identically for the internal and external vortices but is shifted by a half-
period. Variations are generally small with respect to mean velocity
$\overline{U}_{tan}$. As seen on figures 11a and 11b, the mean velocity
displays a linear dependency on $h_{\tau}$ with only a weak dependency on
$H^{*}$ for $H^{*}<10$. Figure 11a also shows that, whatever $h_{\tau}^{*}$
and $H^{*}$, the main contribution to the tangential velocity is the
translation speed $U_{t}$ along $\mathscr{H}$. We suspect that this comes from
the large value of $R^{*}$ that we have considered ($R^{*}=9.5$).

It is also interesting to compare $U_{t}$ to the velocity associated with the
moving frame

$U^{F}_{t}=c_{\tau}U_{z}^{F}+c_{\kappa}R^{*}\Omega^{F}$ (17)

since the difference between the two is issued from vortex induction. For
leapfrogging and weakly braided configurations ($h_{\tau}^{*}>10$), most of
the tangential velocity is associated with the moving frame and the role of
vortex induction is generally small. Note in particular that when $H^{*}$ is
also large, we recover the estimate $\overline{U}_{tan}\approx
h_{\tau}\Gamma/(2\pi^{2}d^{2})$ obtained from the rectilinear approximation in
this limit. By contrast, for densely coiled wakes, such as $h_{\tau}^{*}=2$,
$\overline{U}_{tan}$ and $\overline{U}_{t}$ can be significantly larger than
the frame velocity and even point opposite directions (figure 11b).

####  III.2.4 Properties of the induced flow

(a) (b) (c)

Figure 12:  Axial velocity $\frac{d}{\Gamma}U_{z}$ (top) and angular velocity
$\frac{d}{\Gamma}U_{\theta}$ at the plane $(x,y,z=0)$ for $H^{*}=15$. (a)
Leapfrogging wake for $\beta=1$, (b) sparsely braided wake for $\beta=4$, and
(c) densely coiled wake for $\beta=16$. Solid lines indicate the intersection
with the vortex cores, while dashed lines correspond to a helical vortex tube
enclosing the twin-vortex.

(a) (b)

Figure 13:  Azimuthally averaged velocity profiles for $H^{*}=15$ and
$h_{\tau}^{*}=2$. Fig. (a) (resp. (b)) Rotation rate and axial velocity
component in the global cylindrical (resp. local helical) frame. Dashed lines
correspond to expressions (18) (resp. (19)).

For applications, it is useful to know the velocity field induced by the
vortex structure. In this section, the induced flow is provided in cross
planes and compared to the predictions obtained from the Kawada-Hardin
solutions for perfect helices [35, 36, 50]. Figures 12a-c show the axial and
angular velocity components taken at the plane $z=0$ (where both vortices have
the same azimuth) for a leapfrogging, sparsely braided and densely coiled
structures, respectively. In these plots, we subtract the frame velocities
such that the velocity field vanishes far from the center. Vortex cores
intersect the plane on numerous occasions resulting in inhomogeneous velocity
and vorticity distributions. For densely coiled structures, this results in
regions of very intense angular velocity and opposite signed vorticity.
Nevertheless, away from the cores (outside the region enclosed in dashed
lines) the induced velocity is reminiscent to that of a helical vortex.

To highlight this similarity, consider the azimuthally averaged velocity
profiles presented in figure 13a. Inside the region
$R^{*}-\frac{1}{2}<r<R^{*}+\frac{1}{2}$, the angular velocity displays a
maximum around $r=R^{*}$, while the axial velocity can be roughly approximated
by a linear function. Outside the same region the azimuthally averaged
profiles approach the ideal profiles of a uniform helix obtained from the
Kawada-Hardin solutions [35, 36, 50]

$\left\langle U_{z}^{H}\right\rangle_{\theta}=\begin{cases}2\Gamma/H&\mbox{ if
}r<R^{*}-1/2\\\ 0&\mbox{ if
}r>R^{*}+1/2\end{cases}~{};\quad\quad\left\langle\Omega^{H}\right\rangle_{\theta}=\begin{cases}0&\mbox{
if }r<R^{*}-1/2\\\ \Gamma/(\pi r^{2})&\mbox{ if }r>R^{*}+1/2\end{cases}.$ (18)

Close to the vortex cores, the velocity field can be better understood by
projecting the velocity components into the local helical frame
$(\rho^{\prime},\phi^{\prime})$. This is shown in figure 13b, where the
$\phi^{\prime}$-averaged profiles also approach the mean profiles of a perfect
double-helix locally aligned with $\mathscr{H}$

$\left\langle
U_{t}^{h}\right\rangle_{\phi^{\prime}}=\begin{cases}2\Gamma/h_{\tau}&\mbox{ if
}\rho^{\prime}<(1-\varepsilon)/2\\\ 0&\mbox{ if
}\rho^{\prime}>(1+\varepsilon)/2\end{cases}~{};\quad\quad\left\langle\omega^{h}\right\rangle_{\phi^{\prime}}=\begin{cases}0&\mbox{
if }\rho^{\prime}<(1-\varepsilon)/2\\\ \Gamma/(\pi{\rho^{\prime}}^{2})&\mbox{
if }\rho^{\prime}>(1+\varepsilon)/2\end{cases}.$ (19)

The above results highlight the dual nature exhibited by this vortex
structure. Solutions for closely spaced co-rotating helical vortices issue
from a balance between contributions at different scales with long-distance
effects primarily governed by the large-scale helix and more local features
defined by the double-helix geometry.

##  IV Near-field solutions

In the previous section, we have described the solutions in the far-field. In
this section, we consider the vortex structure close to the rotor plane. The
objective is not to perform a complete parametric study but to illustrate how
these solutions can describe the wake of a tip-splitting rotor under different
flight regimes.

###  IV.1 Numerical method

There is one notable difference with respect to the analysis performed in the
previous section. Instead of fixing the geometric parameters to find the
associated frame velocities, now we solve the inverse problem. That is, we fix
the operating conditions, i.e., the rotation rate $\Omega^{F}=\Omega_{R}$ and
external velocity $U_{z}^{F}=-U_{\infty}$, and compute the corresponding wake
geometry. Now the prescribed geometrical parameters are the radial coordinates
$r=R_{fin}=R_{0}-d_{0}/2$ and $r=R_{tip}=R_{0}+d_{0}/2$ where the vortices are
emitted on the rotor blade (located at $\theta_{0}=0$ and $z=0$), and the
vortex core size $a$. As previously, we define four non-dimensional parameters

$\lambda\equiv\frac{R_{0}\Omega_{R}}{U_{\infty}},\quad\quad\eta\equiv\frac{\Gamma}{R_{0}^{2}\Omega_{R}},\quad\quad
R_{0}^{*}\equiv\frac{R_{0}}{d_{0}},\quad\quad\varepsilon_{0}^{*}\equiv\frac{a}{d_{0}},$
(20)

where $\lambda$ is known as the tip-speed ratio and $\eta$ represents the
relative vortex strength. As in [39], we choose a convention such that $\eta$
remains positive, while $\lambda$ may change sign. Wind turbine regimes will
always correspond to positive values of $\lambda$. By contrast, helicopters
may correspond to either negative values of $\lambda$ in climbing flight or
positive values of $\lambda$ in descending flight. In the present study, only
the parameters $\eta$ and $\lambda$ are varied. In most cases,
$\varepsilon_{0}^{*}$ and $R_{0}^{*}$ are fixed to 0.1 and 9.5, respectively.
Note that these parameters are different from the parameters $R^{*}$ and
$\varepsilon^{*}$ defined above for the far-field, as the mean radius and the
vortex separation distance in the far-field are now the results of the
calculation.

Figure 14:  Schematic representation of the semi-infinite case: two co-
rotating vortices are emitted from a rotor blade at radial positions $R_{fin}$
and $R_{tip}$, respectively. The calculation domain and the prescribed far-
wake are shown in different colors.

The near-field solution satisfies (4) with the boundary condition at
$\zeta=0$: $r_{1}=R_{0}^{*}-\frac{1}{2}$ and $r_{2}=R_{0}^{*}+\frac{1}{2}$ and
$\theta_{j}=z_{j}=0$. As we go away from the rotor, it should match a far-
field solution. To implement this condition, we follow the numerical procedure
used by Durán Venegas  _et al._ [39] for Joukowski’s rotor wake model. We
decompose the induced velocity into contributions from the near-field and far-
field. Contributions from the far-field are modeled by imposing that after a
certain distance from the rotor plane, the wake adopts the geometry of a far-
field solution (see figure 14). The geometrical parameters of this far-wake
solution are estimated from the near-field solution at the end of calculation
domain. The computational domain must be large enough for the wake to develop
and match the far-field. In practice, the size of the domain and of the
corresponding far-wake range from $30\pi$ to $60\pi$.

###  IV.2  Description of the steady solutions

(a) (b)

Figure 15:  Schematic representation of the vortex wake and evolution of the
radial coordinates $r_{j}$ as function of $z$ for $\eta=0.04$,
$R_{0}^{*}=9.5$, $\varepsilon_{0}^{*}=0.1$. (a) $\lambda=5.5$ in a wind
turbine regime; (b) $\lambda=40$ in a helicopter descending flight regime.
Markers indicate the positions where both vortices have the same azimuth.

Figure 15 displays the wake structure with the corresponding radial
coordinates as function of the axial coordinate $z$. Figure 15a (resp. 15b) is
representative of the wind turbine (resp. helicopter) regime, which displays a
radial expansion (resp. contraction), while the two vortices trade places
continuously. Simultaneously, the separation distance may also contract or
expand due to the combined effects of the vortex pair and the rotor’s
presence. Additionally, the structure may exhibit local variations of both
axial pitches. For clarity, we characterize the wake in terms of the geometric
parameters introduced in §III.

####  IV.2.1 Evolution of the wake in the near-field

(a) (b)

Figure 16:  Evolution of (a) $R_{n}/R_{0}$, and (b) $H_{n}/R_{0}$ as function
of $|z_{n}|/R_{0}$ for different values of $\lambda$ for $\eta=0.04$,
$R_{0}^{*}=9.5$, $\varepsilon_{0}^{*}=0.1$. For each case, we display the
corresponding evolution obtained for an ‘equivalent’ helical vortex in dashed
lines.

As in §III, we introduce the curve $\mathbf{X}_{mid}(\zeta)$ equidistant to
both vortices and the separation distance $d(\zeta)$ defined by

$\mathbf{X}_{mid}(\zeta)=(\mathbf{X}_{tip}(\zeta)+\mathbf{X}_{fin}(\zeta))/2,$
(21) $d(\zeta)=\lVert\mathbf{X}_{tip}(\zeta)-\mathbf{X}_{fin}(\zeta)\rVert.$
(22)

We also perform a series of discrete measurements at the points along the
curve $(\zeta_{n},n=1,2,\cdots)$ where the internal and external vortices have
the same azimuth (indicated by solid marks in figure 15). In this way, the
large helix $\mathscr{H}$ is characterized by a set of radii and pitches

$R_{n}\equiv r_{mid}(\zeta_{n}),\quad\quad\quad\quad\quad\quad\quad
H_{n}\equiv\frac{z_{mid}(\zeta_{n+1})-z_{mid}(\zeta_{n})}{\theta_{mid}(\zeta_{n+1})-\theta_{mid}(\zeta_{n})},\quad\mbox{
for }n=0,1,2,\cdots$ (23)

which converge towards the values at the far-field $R_{\infty}$ and
$H_{\infty}$, respectively.

Figure 16a displays the evolution of $R_{n}/R_{0}$ as a function of
$z(\zeta_{n})/R_{0}$ for cases representative of wind turbine and helicopter
regimes. Far-field solutions suggest that $\mathscr{H}$ can be modeled by a
single tip vortex of circulation $2\Gamma$ and core size $d$ emitted from the
radial coordinate $R_{0}$ under the same operating conditions. A good
agreement is observed between the evolution of $R_{n}/R_{0}$ and the
equivalent helix approximation. A similar agreement is observed between
$H_{n}/R_{0}$ and the ‘equivalent’ pitch in figure 16b.

Conversely, the rotation of the vortex pair is characterized by a set of mean
separation distances and axial pitches

$d_{n}\equiv d(\zeta_{n}),\quad
h_{n}\equiv\frac{z_{mid}(\zeta_{n+1})-z_{mid}(\zeta_{n})}{\pi}\quad,\mbox{ for
}n=0,1,2,\cdots$ (24)

which converge towards the values at the far-field $d_{\infty}$ and
$h_{\infty}$, respectively. Combining both sets of measurements, we obtain a
set for the mean pitch and twist parameter

$h_{\tau,n}\equiv
h_{n}\frac{\sqrt{H_{n}^{2}+4\pi^{2}R_{n}^{2}}}{H_{n}},\quad\beta_{n}\equiv\frac{H_{n}}{h_{n}}\quad,\mbox{
for }n=1,2,\cdots$ (25)

which converge towards the values at the far-field $h_{\tau,\infty}$ and
$\beta_{\infty}$.

(a) (b) (c)

Figure 17:  Evolution of (a) $d_{n}/d_{0}$, (b) $h_{\tau,n}/d_{0}$ and
$\beta_{n}$ as a function of $z_{n}/R_{0}$ for the cases displayed in figure
16. Dashed lines in figures (a) and (b), correspond to (28) and (29),
respectively. For reference, the evolution of $d(z)/d_{0}$ is shown in solid
lines in fig. (a).

The evolution of these parameters is shown in figures 17a-c. The first
remarkable feature is a dynamical similarity between local parameters’
evolution and those of the large helix. Variations are limited to the very
near wake, and are rapid whenever those of the large helix are rapid. However,
variations seem to be opposite for the small and the large helix: the
separation distance $d_{n}$ increases (resp. decreases) when the large helix
radius $R_{n}$ decreases (resp. increases). The local pitch $h_{\tau,n}$ is
also found to have an opposite sign to $H_{n}$, meaning that the small and
large helices are oriented differently. Additionally, both pitches remain
constant in the wake. It is also important to note the large absolute value of
$h_{\tau,n}/d_{0}$, indicating that configurations are weakly twisted where
self-induced effects are generally small. This has some consequences that can
be used to derive a model explaining the variations of the local properties.

The evolution of the separation distance $d_{n}$ can be understood by the
following arguments. As the radius $R_{n}$ of the large helix varies, we
expect variations of the advection speed $U_{t}$ of the vortices along
$\mathscr{H}$. These variations are responsible of a local strain that tends
to contract the vortex structure (that is, decrease $d$) when $U_{t}$
increases, and expand it when $U_{t}$ decreases. The variations of $d$ are
then obtained by the conservation of the mass flux in the vortex stream tube
which prescribes that

$U_{t}(\zeta_{n})d^{2}(\zeta_{n})=cst.$ (26)

An estimate for $U_{t}(\zeta_{n})$ can be obtained when the self-induction is
weak. We have seen in §III that in that case, $U_{t}(\zeta_{n})$ is mainly
provided by the frame velocity
$U_{t}^{F}(\zeta_{n})=\pm(U_{\infty}^{2}+R_{n}^{2}\Omega_{R}^{2})^{1/2}$ where
the $\pm$ sign is defined by the sign of $U_{\infty}$. In terms of the non-
dimensionalized parameters $\lambda$ and $\eta$, this gives

$U_{t}(\zeta_{n})\frac{R_{0}}{\Gamma}\approx\frac{\sqrt{1+\lambda^{2}R_{n}^{2}/R_{0}^{2}}}{\eta\lambda}$
(27)

Combining this expression with (26), we obtain the following relation between
$d_{n}$ and $R_{n}$

$\frac{d_{n}}{d_{0}}=\left[\frac{1+\lambda^{2}\quad\quad\quad}{1+\lambda^{2}R_{n}^{2}/R_{0}^{2}}\right]^{1/4}_{.}$
(28)

This expression is tested in figure 17a, where a contraction (resp. expansion)
of the separation distance is observed for the wind turbine (resp. helicopter)
regimes. While equation (28) gives correct trends and orders of magnitude, it
tends to overestimate the expansion and underestimate the contraction. This
difference could be explained in different ways. First, the contributions from
the double-helix to $U_{t}$ estimated at the far-field from (19) as
$2R_{0}^{*}{d_{0}}/{h_{\tau,n}}$ have been neglected in (27) since vortex
induction is rather weak. Second, a radial contraction of the double-helix
close to the blade that should also be present but that has not been taken
into account. Such contraction could be estimated by analogy to a vortex pair
under a locally axial flow $U_{t}$. Considering these two effects, would give
an additional contraction as observed in the numerics, but still would be too
small to account for the discrepancy. Another possibility is the differential
expansion/contraction experienced by the two vortices close to the rotor plane
due to their radial position. The outermost vortex would experience a larger
expansion/contraction than the innermost vortex, thus changing the separation
distance. However, since the two vortices continuously trade places with one
another, this effect is hard to quantify.

In a similar vein, the pitch of the double-helix can be estimated using
$h_{\tau}=2\pi{U_{t}}/{\omega}$ where $\omega$ is the rotation rate of the
vortex pair. As shown in figure 3, for large values of $h_{\tau}/d$ we can
approximate $\omega$ as $\omega=\Gamma/\pi d^{2}$. Combining both expressions
shows that $h_{\tau,n}$ should be constant and given by

$\frac{h_{\tau,n}}{d_{0}}\approx\frac{2\pi^{2}}{\eta
R_{0}^{*}}\left[\frac{\sqrt{1+\lambda^{2}}}{\lambda}\right]$ (29)

which is verified in figure 17b. Equation (29) also suggests that
$h_{\tau,n}/d_{0}$ is primarily governed by $R_{0}^{*}$ and $\eta$, since the
part shown in brackets remains close to 1 (in absolute value). The condition
$h_{\tau,n}/d_{0}>10$ for being weakly twisted then reduces to $\eta
R_{0}^{*}\lessapprox 2$. For $\eta=0.04$, this condition is satisfied as soon
as $R_{0}^{*}\lessapprox 50$.

(a) (b) (c)

Figure 18:  (a) Sparsely braided wakes obtained for $\eta=0.04$,
$R_{0}^{*}=19.5$, $\varepsilon_{0}^{*}=0.1$ in a wind turbine regime
($\lambda=5$) and a helicopter climbing flight regime ($\lambda=-40$). (b,c)
Evolution of $R_{n}/R_{0}$ and $\beta_{n}$ compared to the equivalent helix
approximation in dashed lines.

From (29), we can also compute the twist parameter
$\beta_{n}=H_{n}/(h_{\tau,n}c_{\tau}(\zeta_{n}))$ as

$\beta_{n}\approx-\frac{\eta(R_{0}^{*})^{2}}{\pi}\left[\frac{|\lambda|}{\sqrt{1+\lambda^{2}}}\frac{\sqrt{H_{n}^{2}+4\pi^{2}R_{n}^{2}}}{2\pi
R_{0}}\right].$ (30)

which has been compared to the data in figure 17c. As the part shown in
brackets remains in general close to 1, the typical order of $|\beta|$ is then
given by the product $\eta(R_{0}^{*})^{2}$. In particular, we expect to be in
a leapfrogging situation, i.e. $|\beta|<1$, as soon as
$\eta(R_{0}^{*})^{2}\lessapprox\pi$. Equation (30) suggests a simple way to
control the twist parameter by increasing $R_{0}^{*}$ (i.e. reducing the
separation distance $d_{0}$). In order to illustrate this point, we consider
the same values of $(\lambda,\eta)$ for $R_{0}^{*}=19.5$ in different flight
regimes and measure $\beta_{n}$. As expected, the evolution of the large-scale
pattern is comparable to the equivalent single tip rotor, while $\beta_{n}$ is
roughly four times larger with respect to the case with $R_{0}^{*}=9.5$
(compare figure 17c and figure 18c).

(a) (b)

(c) (d)

Figure 19:  Evolution of the far-wake parameters as a function of $1/\lambda$
for $\eta=0.04$, $\varepsilon_{0}^{*}=0.1$. Dashed lines correspond to the
equivalent vortex approximations in (a,b) and to equations (29) and (30) in
(c,d). For reference, fig. (c) also displays the minimum and maximum values of
$d(\zeta)/d_{0}$ for $\zeta\to\infty$ as a vertical bar. Helicopter regimes
are on the left, wind turbine regimes on the right.

(a) (b)

(c) (d)

Figure 20:  Evolution of the far-wake parameters as a function of $\eta$ for
$\lambda=2$ (wind turbine regime), $R_{0}^{*}=9.5$ and
$\varepsilon_{0}^{*}=0.1$. Dashed lines correspond to the equivalent vortex
approximations in (a,b) and to equations (29) and (30) in (c,d). For
reference, fig. (c) also displays the minimum and maximum values of
$d(\zeta)/d_{0}$ for $\zeta\to\infty$ as a vertical bar. The vertical dotted
line in (c) indicates the value $\eta=\pi/(R_{0}^{*})^{2}$ where
$\beta_{\infty}\approx 1$.

The above approximations can also be used to predict the characteristics of
the far-wake. In figures 19 and 20, we have compared the model with the
numerics as both $\lambda$ and $\eta$ vary. We can observe in figures 19a,b
and 20a,b that for the radius $R_{\infty}$ and the pitch $H_{\infty}$, the
agreement with the equivalent vortex approximation is excellent for all
regimes. For the separation distance $d_{\infty}$ (figures 19c and 20c), the
trend is qualitatively good but there is a constant systematic over-estimation
in the model, that has already been noticed above. The fact that $d_{\infty}$
changes little with $R_{0}^{*}$ is also consistent with (29). Interestingly,
the twist parameter $\beta_{\infty}$ (figures 19d and 20d) is always found to
be well-predicted, probably because it does not directly depend on the
separation distance.

##  V Conclusions

In this article, we have provided a vortex wake model for a tip-splitting
rotor in any incident normal wind, which applies to all the vertical
helicopter flight and wind turbine regimes. This model is based on a
simplified description of the wake with a pair of co-rotating vortex filaments
attached to each blade at different radial positions. A free-vortex method
together with Biot & Savart law has been used to obtain steady solutions in a
frame moving with the rotor blades.

We have first analyzed the structure of the solutions in the far-field. As for
counter-rotating vortices [37], we have shown that steady periodic solutions
can be obtained in specific reference frames and that the properties of these
solutions are governed by 4 purely geometrical dimensionless parameters
($R/d$, $H/d$, $h_{\tau}/d$ and $a/d$). We have documented the fluctuations of
the separation distance and of the mean radius, as well as the variations of
the frame velocities with respect to two of these parameters ($H/d$ and
$h_{\tau}/d$). We have shown that these solutions can be understood in terms
of two interlaced helical vortices inscribed on top of a larger helical
structure, which allowed us to obtain good estimates for the wake geometry and
the frame velocities.

The far-field solutions have then been used to get near-wake solutions. We
have solved the semi-infinite free-vortex problem with a prescribed condition
on the blade (the vortices are attached at given positions) and a condition of
matching with a far-field solution away from the rotor. We have been able to
obtain steady solutions in the frame rotating with the rotor, that capture the
expected global expansion/contraction of the wake as well as the deformations
induced by the mutual interaction of the vortices. Interestingly, we have also
been able to show that the solution is well described by a small-scale twisted
structure inscribed onto a large-scale pattern obtained for a single large
tip-vortex with the total circulation. We have further shown how the
variations of the parameters of the small-scale structure are driven by the
large-scale motion.

In applications, the core size of tip vortices is sufficiently small to make
the filament approach relevant. For instance, in the MEXICO project, a wind
turbine of diameter 4.5 m creates a tip vortex of radius 2 cm giving a ratio
$a/R\sim 0.01$ [51]. Similar values can be found in helicopter rotor wakes
[52]. In such cases, there is clearly enough room to generate a second tip
vortex at a small distance $d$ from the tip such that $a/d$ remains small. We
therefore claim that the typical parameters $a/d=0.1$ and $R/d=10$ that we
have considered in the present study can be obtained in real applications.
However, it does not mean that our analysis will automatically apply. The
roll-up process giving rise to the vortices could be significantly more
complex. In a recent experimental study, Schröder  _et al._ [19] succeeded in
creating a second vortex close to the tip using a fin. However, the fin vortex
was also found to contain negative vorticity contributions coming from the
vorticity sheet shed by the modified blade, that made it unstable with respect
to the centrifugal instability. In that case, the second vortex is thus
rapidly disrupted and our solution cannot be formed. We suspect that if the
blade was more tapered close to the tip, the fin vortex would not have been
unstable and a different evolution would have been observed.

Our solutions are nevertheless expected to be unstable. Helical solutions are
indeed known to be unstable with respect to a long-wavelength instability, the
so-called Widnall instability [53, 54]. As shown in Quaranta  _et al._ [11]
and Huang  _et al._ [12], this instability can be excited to accelerate the
destruction of the vortices. In a follow-up paper [55], the linear stability
of our solutions is analyzed using the filament framework. Owing to the higher
complexity of our solutions, we demonstrate that there exist different types
of instability modes. As for the Widnall instability [10, 11], each
instability mode is nevertheless shown to be associated with a specific
pairing event.

Short-wavelength instabilities such as the curvature instability [56] and the
elliptic instability [57, 58] are also expected to develop in the core of the
vortices. These instabilities cannot be described by the filament approach
which neglects everything occurring in the vortex cores. It requires
monitoring the core deformations induced by curvature, torsion and strain
fields which are responsible for these instabilities. This can be done by
matched asymptotic techniques [59, 60] or by direct numerical simulations when
the base flow exhibits particular symmetries as shown by [40] for 2D vortex
pairs, [61] for rings and [62] for helices. The internal core structure that
is obtained can then be used for the short-wavelength stability study (see for
instance [63, 61]).

Note that the determination of a quasi-steady solution may not always be
necessary for the stability study if the time scale of the instability is long
compared to the relaxation time scale needed to get the correct core
structure. Brynjell-Rahkola and Henningson [64] have indeed shown that they
could obtain the correct stability properties of a helix with a Gaussian
vorticity core by just analysing the temporal evolution of the numerical
solution obtained from the initial condition formed by the filament structure
with the prescribed axisymmetric vortex profile in the core. This is an
interesting result as it provides a simple way to analyse the short-wavelength
stability properties of our solutions.

## Acknowledgment

This work is part of the French-German project TWIN-HELIX, supported by the
Agence Nationale de la Recherche (grant no. ANR-17-CE06-0018) and the Deutsche
Forschungsgemeinschaft (grant no. 391677260).

## References

  * Young [2003] L. A. Young, Vortex core size in the rotor near-wake, NASA Technical Report TM-2003-212275 (2003).
  * Hansen  _et al._ [2012] K. S. Hansen, R. J. Barthelmie, L. E. Jensen, and A. Sommer, The impact of turbulence intensity and atmospheric stability on power deficits due to wind turbine wakes at horns rev wind farm, Wind Energy 15, 183 (2012).
  * Barthelmie and Jensen [2010] R. J. Barthelmie and L. Jensen, Evaluation of wind farm efficiency and wind turbine wakes at the nysted offshore wind farm, Wind Energy 13, 573 (2010).
  * Gandhi and Tauszig [2000] F. Gandhi and L. Tauszig, Influence of individual interactions on helicopter blade-vortex interaction noise, in  _Proceedings of the 26th European Rotorcraft Forum_ (2000).
  * Hardin and Lamkin [1987] J. C. Hardin and S. L. Lamkin, Concepts for reduction of blade/vortex interaction noise, J. Aircr. 24, 120 (1987).
  * Yu [2000] Y. H. Yu, Rotor blade-vortex interaction noise, Prog. Aerosp. Sci. 36, 97 (2000).
  * Brocklehurst and Barakos [2013] A. Brocklehurst and G. Barakos, A review of helicopter rotor blade tip shapes, Prog. Aerosp. Sci. 56, 35 (2013).
  * Fleming  _et al._ [2017] P. Fleming, J. Annoni, J. J. Shah, L. Wang, S. Ananthan, Z. Zhang, K. Hutchings, P. Wang, W. Chen, and L. Chen, Field test of wake steering at an offshore wind farm, Wind Energy 2, 229 (2017).
  * Campagnolo  _et al._ [2020] F. Campagnolo, R. Weber, J. Schreiber, and C. L. Bottasso, Wind tunnel testing of wake steering with dynamic wind direction changes, Wind Energy 5, 1273 (2020).
  * Quaranta  _et al._ [2015] H. U. Quaranta, H. Bolnot, and T. Leweke, Long-wave instability of a helical vortex, J. Fluid Mech. 780, 687 (2015).
  * Quaranta  _et al._ [2019] H. U. Quaranta, M. Brynjell-Rahkola, T. Leweke, and D. S. Henningson, Local and global pairing instabilities of two interlaced helical vortices, J. Fluid Mech. 863, 927 (2019).
  * Huang  _et al._ [2019] X. Huang, S. A. Moghadam, P. Meysonnat, M. Meinke, and W. Schröder, Numerical analysis of the effect of flaps on the tip vortex of a wind turbine blade, Int J Heat Fluid Flow 77, 336 (2019).
  * Ito [2009] S. Ito, Aerodynamic influence of leading-edge serrations on an airfoil in a low Reynolds number, J. Biomech. Sci. Eng. 4, 117 (2009).
  * Pang  _et al._ [2018] E. Pang, A. Cambray, D. Rezgui, M. Azarpeyvand, and S. A. Showkat Ali, Investigation towards a better understanding of noise generation from UAV propellers, in  _2018 AIAA/CEAS Aeroacoustics Conference_ (2018) p. 3450.
  * Chattot [2009] J.-J. Chattot, Effects of blade tip modifications on wind turbine performance using vortex model, Comput Fluids 38, 1405 (2009).
  * Ebrahimi and Mardani [2018] A. Ebrahimi and R. Mardani, Tip-vortex noise reduction of a wind turbine using a winglet, J. Energy Eng. 144, 04017076 (2018).
  * Han and Leishman [2004] Y. O. Han and J. G. Leishman, Investigation of helicopter rotor-blade-tip-vortex alleviation using a slotted tip, AIAA journal 42, 524 (2004).
  * Brocklehurst and Pike [1994] A. Brocklehurst and A. Pike, Reduction of BVI noise using a vane tip, in  _AHS Aeromechanics Specialists Conference_ (American Helicopter Society, 1994).
  * Schröder  _et al._ [2021] D. Schröder, T. Leweke, R. Hörnschemeyer, and E. Stumpf, Experiments on helical vortex pairs in the wake of a rotor, in  _AIAA Scitech 2021 Forum_ (2021) p. 1088.
  * Copland  _et al._ [1998] C. Copland, F. Coton, and R. M. Galbraith, An experimental study of the idealised vortex system of a novel rotor blade tip, Aeronaut. J 102, 385 (1998).
  * Coton  _et al._ [2005] F. N. Coton, R. B. Green, and R. M. Galbraith, Analysis of model rotor blade pressures during parallel interaction with twin vortices, J. Aircr. 42, 1552 (2005).
  * Cho  _et al._ [2006] C.-H. Cho, D.-J. Lee, C. Hwang, and Kari, Numerical investigations of parallel BVI noise with split tip vortices, in  _AHS Internaitonal 62nd Annual Forum_ , Vol. I (American Helicopter Society, 2006) pp. 13 – 19.
  * Jung  _et al._ [2008] J. Jung, W. Y. Choi, and S. Lee, Effect of the blade inner-tip position on the generation of twin vortices, J. Aircr. 45, 728 (2008).
  * Sanderse  _et al._ [2011] B. Sanderse, S. Van der Pijl, and B. Koren, Review of computational fluid dynamics for wind turbine wake aerodynamics, Wind energy 14, 799 (2011).
  * Stevens and Meneveau [2017] R. J. Stevens and C. Meneveau, Flow structure and turbulence in wind farms, Annual review of fluid mechanics 49, 311 (2017).
  * Caprace  _et al._ [2020] D.-G. Caprace, P. Chatelain, and G. Winckelmans, Wakes of rotorcraft in advancing flight: A large-eddy simulation study, Phys. Fluids 32, 087107 (2020).
  * Branlard [2020] E. Branlard,  _Wind turbine aerodynamics and vorticity-based methods_, Vol. 7 (Springer, 2020).
  * Sørensen [2016] J. N. Sørensen,  _General momentum theory for horizontal axis wind turbines_, Vol. 4 (Springer, 2016).
  * Van Kuik  _et al._ [2015] G. Van Kuik, J. N. Sørensen, and V. Okulov, Rotor theories by Professor Joukowsky: momentum theories, Prog. Aerosp. Sci. 73, 1 (2015).
  * Shen  _et al._ [2009] W. Z. Shen, M. O. Hansen, and J. N. Sørensen, Determination of the angle of attack on rotor blades, Wind Energy 12, 91 (2009).
  * Sørensen [2011] J. N. Sørensen, Aerodynamic aspects of wind energy conversion, Annu. Rev. Fluid Mech. 43, 427 (2011).
  * Goldstein [1929] S. Goldstein, On the vortex theory of screw propellers, Proc. R. Soc. Lond. A 123, 440 (1929).
  * Okulov and Sørensen [2007] V. L. Okulov and J. N. Sørensen, Optimum operating regimes for the ideal wind turbine, in  _J. Phys. Conf. Ser._, Vol. 75 (IOP Publishing, 2007) p. 012009.
  * Leishman  _et al._ [2002] J. G. Leishman, M. J. Bhagwat, and A. Bagai, Free-vortex filament methods for the analysis of helicopter rotor wakes, J. Aircr. 39, 759 (2002).
  * Kawada [1936] S. Kawada, Induced velocity by helical vortices, J. Aeronaut. Sci. 3, 86 (1936).
  * Hardin [1982] J. C. Hardin, The velocity field induced by a helical vortex filament, Phys. Fluids 25, 1949 (1982).
  * Durán Venegas and Le Dizès [2019] E. Durán Venegas and S. Le Dizès, Generalized helical vortex pairs, J. Fluid Mech. 865, 523 (2019).
  * Miller and Bliss [1993] W. O. Miller and D. B. Bliss, Direct periodic solutions of rotor free wake calculations, J Am Helicopter Soc 38, 53 (1993).
  * Durán Venegas  _et al._ [2021] E. Durán Venegas, P. Rieu, and S. Le Dizès, Structure and stability of Joukowski’s rotor wake model, J. Fluid Mech. 911, A6 (2021).
  * Le Dizes and Verga [2002] S. Le Dizes and A. Verga, Viscous interactions of two co-rotating vortices before merging, J. Fluid Mech. 467, 389 (2002).
  * Melander  _et al._ [1988] M. Melander, N. Zabusky, and J. McWilliams, Symmetric vortex merger in two dimensions: causes and conditions, J. Fluid Mech. 195, 303 (1988).
  * Meunier  _et al._ [2002] P. Meunier, U. Ehrenstein, T. Leweke, and M. Rossi, A merging criterion for two-dimensional co-rotating vortices, Phys. Fluids 14, 2757 (2002).
  * Cerretelli and Williamson [2003] C. Cerretelli and C. Williamson, The physical mechanism for vortex merging, J. Fluid Mech. 475, 41 (2003).
  * Josserand and Rossi [2007] C. Josserand and M. Rossi, The merging of two co-rotating vortices: a numerical study, Eur. J. Mech. B Fluids 26, 779 (2007).
  * Leweke  _et al._ [2016] T. Leweke, S. Le Dizes, and C. H. Williamson, Dynamics and instabilities of vortex pairs, Annu. Rev. Fluid Mech. 48, 507 (2016).
  * Saffman [1992] P. G. Saffman,  _Vortex dynamics_ (Cambridge university press, 1992).
  * Leishman [2016] J. Leishman,  _Principles of Helicopter Aerodynamics_ , Cambridge Aerospace Series (Cambridge University Press, 2016).
  * Lamb [1945] H. Lamb,  _Hydrodynamics_ , Dover Books on Physics (Dover publications, 1945).
  * Zhukovskiĭ and Vetchinkin [1929] N. E. Zhukovskiĭ and V. P. Vetchinkin,  _Théorie tourbillonnaire de l’hélice propulsive_ (Gauthier-Villars, 1929).
  * Fukumoto  _et al._ [2015] Y. Fukumoto, V. L. Okulov, and D. H. Wood, The contribution of kawada to the analytical solution for the velocity induced by a helical vortex filament, Appl. Mech. Rev. 67 (2015).
  * Snel  _et al._ [2007] H. Snel, J. Schepers, and B. Montgomerie, The MEXICO project (model experiments in controlled conditions): The database and first results of data processing and interpretation, in  _J Phys Conf Ser_, Vol. 75 (IOP Publishing, 2007) p. 012014.
  * Bauknecht  _et al._ [2017] A. Bauknecht, B. Ewers, O. Schneider, and M. Raffel, Blade tip vortex measurements on actively twisted rotor blades, Exp Fluids 58, 1 (2017).
  * Widnall [1972] S. E. Widnall, The stability of a helical vortex filament, J. Fluid Mech. 54, 641 (1972).
  * Gupta and Loewy [1974] B. Gupta and R. Loewy, Theoretical analysis of the aerodynamic stability of multiple, interdigitated helical vortices, AIAA journal 12, 1381 (1974).
  * Castillo-Castellanos  _et al._ [2021] A. Castillo-Castellanos, E. Durán Venegas, and S. Le Dizès, Closely spaced co-rotating helical vortices: Long-wave instability, Manuscript in preparation (2021).
  * Blanco-Rodríguez and Le Dizès [2017] F. J. Blanco-Rodríguez and S. Le Dizès, Curvature instability of a curved Batchelor vortex, J. Fluid Mech. 814, 397 (2017).
  * Kerswell [2002] R. R. Kerswell, Elliptical instability, Annu. Rev. Fluid Mech. 34, 83 (2002).
  * Blanco-Rodríguez and Le Dizès [2016] F. J. Blanco-Rodríguez and S. Le Dizès, Elliptic instability of a curved Batchelor vortex, J. Fluid Mech. 804, 224 (2016).
  * Callegari and Ting [1978] A. Callegari and L. Ting, Motion of a curved vortex filament with decaying vortical core and axial velocity, SIAM J. Appl. Math 35, 148 (1978).
  * Blanco-Rodríguez  _et al._ [2015] F. J. Blanco-Rodríguez, S. Le Dizès, C. Selçuk, I. Delbende, and M. Rossi, Internal structure of vortex rings and helical vortices, J. Fluid Mech. 785, 219 (2015).
  * Hattori  _et al._ [2019] Y. Hattori, F. Blanco-Rodríguez, and S. Le Dizès, Numerical stability analysis of a vortex ring with swirl, J. Fluid Mech. 878, 5 (2019).
  * Selçuk  _et al._ [2017] C. Selçuk, I. Delbende, and M. Rossi, Helical vortices: linear stability analysis and nonlinear dynamics, Fluid Dyn. Res. 50, 011411 (2017).
  * Roy  _et al._ [2008] C. Roy, N. Schaeffer, S. Le Dizès, and M. Thompson, Stability of a pair of co-rotating vortices with axial flow, Phys. Fluids 20, 094101 (2008).
  * Brynjell-Rahkola and Henningson [2020] M. Brynjell-Rahkola and D. S. Henningson, Numerical realization of helical vortices: application to vortex instability, Theor. Comp. Fluid Dyn 34, 1 (2020).