# Mode selection in concentric jets. The steady-steady 1:2 resonant mode
interaction with O(2) symmetry.

A. Corrochano  School of Aerospace Engineering, Universidad Politécnica de
Madrid, Madrid 28040, Spain  J. Sierra-Ausín  Institut de Mécanique des
Fluides de Toulouse (IMFT), Toulouse 31400, France  DIIN, Universit´a degli
Studi di Salerno, Via Giovanni Paolo II, 84084 Fisciano (SA), Italy  J. A.
Martin  School of Aerospace Engineering, Universidad Politécnica de Madrid,
Madrid 28040, Spain  D. Fabre  Institut de Mécanique des Fluides de Toulouse
(IMFT), Toulouse 31400, France  S. Le Clainche Email address for
correspondence<EMAIL_ADDRESS>School of Aerospace Engineering,
Universidad Politécnica de Madrid, Madrid 28040, Spain

###### Abstract

In this article, a thorough characterization of the configuration composed by
two concentric jets at a low Reynolds number is presented. The analysis
comprises a layout with a wide range for the velocity ratio between the inner
and outer jets, defined within the interval $[0,2]$, and also details the
influence of the distance between jets, where the wall thicknesses separating
the two jets is $[0.5,4]$. Global linear stability analysis identifies the
most significant modes driving the changes in the flow dynamics. The neutral
lines revealing the critical Reynolds number connected to the presence of the
main (steady and unsteady) flow bifurcations, which are presented by global
azimuthal modes, show the high complexity of the problem under study, where
hysteresis and other types of complex cycles are pointed out. Finally, the
mode interaction is analysed, highlighting the presence of travelling waves
emerging from the interaction of steady states, and the existence of robust
heteroclinic cycles that are asymptotically stable. The high level of detail
in the results presented, makes this work as a reference for future research
development in the field of concentric jets.

##  1 Introduction

Double concentric jets is a configuration enhancing the turbulent mixing of
two jets, which is used in several industrial applications where the breakup
of the jet into droplets due to flow instabilities is presented as the key
technology. Combustion (i.e.: combustion chamber of rocket engines, gas
turbine combustion, internal combustion engines, etc.) and noise reduction
(e.g.: in turbofan engines) are the two main applications of this geometry,
although the annular jets can also be found in some other relevant
applications such as ink-jet printers or spray coating.

The qualitative picture emerging from this type of flow divides the inner
field of concentric jets in three different regions: (i) initial merging zone,
(ii) transitional zone and (iii) merged zone, as presented in fig. 1, that
follows the initial sketch presented by [12].

Figure 1:  Sketch representing the three flow regimes in the near field of
double concentric jets. Figure based on the sketch presented in [12, 31].

In the initial merging zone (i), just at the exit of the two jets, two
axisymmetric shear layers (inner and outer boundary layer) develop and start
to merge. In this region, we distinguish the inner and outer shear layers,
related with the inner and outer jet stream. Then, most of the mixing occurs
in the transitional zone (ii), that extends until the external shear layer
reaches the centreline. Finally, in the merged zone (iii), the two jets are
totally merged, modelling a single jet flow.

Several parameters define the characteristic of this flow: the inner and outer
jet velocities, the jet diameters, the shape and thickness of the wall
separating both jets, the Reynolds number, the boundary layer state and
thickness at the jet exit and the free stream turbulence. Based on these
parameters, it is possible to identify several types of flow behaviour, which
can be related with the presence of flow instabilities.

Ko & Kwan (1976) [12] postulated that the double concentric jet configuration
could be considered as a combination of single jets. Nevertheless, Dahm et al.
(1992) [7] revealed by means of flow visualizations, several topology patterns
as function of the outer/inner jet velocity ratio, reflecting that the
dynamics of the inner and outer jet shear layers were different from that in a
single jet. Moreover, this study exhibited a complex interaction between
vortices identified in both shear layers, affecting the instability mechanism
of the flow. Buresti et al. (1992) [4] found that the outer shear layer
dominated the flow dynamics for cases in which the outer velocity was much
larger than the inner velocity. These authors also detected the presence of an
alternate vortex shedding when the wall thickness between the two jets was
sufficiently large. The same mechanism was recognised by other authors [7,
20]. Rehab et al. (1997) [25] studied in detail the flow differences as
function of the outer/inner velocity ratio, finding two different flow regimes
when the external jet diameter is much larger than the internal jet one. When
the outer/inner velocity ratio was larger than a critical value, the authors
spotted a low frequency recirculation bubble at the jet outlet. On the
contrary, for outer/inner velocity ratio smaller than such critical value, the
outer (still fast) jet excites the inner jet, which ends oscillating at the
same frequency as the external jet. This is known as the lock-in phenomenon.
Moreover, the oscillation frequency detected was similar to the one defined by
a Kelvin-Helmholtz flow instability, which is generally encountered in single
jets. This lock-in phenomenon was also identified by other authors [7, 6].

Following previous works [4, 7, 20] and paying especial attention to the
separating wall thickness and the vortex shedding located behind the wall,
Wallace & Redekopp (1992) [33] showed that the wall thickness and sharpness
change the characteristic of the jet. Segalini & Talamelli (2011) [26]
performed experiments to inspect in detail the effects of the outer/inner
velocity ratio and the wall thickness in double concentric jets. These authors
found that for small outer/inner velocity ratios, the inner jet presents its
own flow instability in the shear layer, while a different flow instability
was identified in the outer jet. On the contrary, for large outer/inner
velocity ratios, the outer shear layer drives the flow dynamics, forcing the
inner shear layer to oscillate with the same frequency, occurring in the lock-
in phenomenon previously mentioned. Finally, for similar outer/inner velocity
ratios, a Von Kármán vortex street was detected near the separating wall, as
also depicted by other authors [4, 7, 20]. A wake instability affected the
inner and outer shear layers, reversing the lock-in phenomenon.

Different configurations can also be found, changing the velocity ratio
between jets. Williams et al. (1969) [34] worked on the influence of the
exterior/interior velocity ratio on noise attenuation, which was analysed
experimentally. It was observed that for some given configurations, more noise
attenuation was present than for the others, with a maximum between $12$ and
$15dB$.

Talamelli & Gavarini (2006) [31] performed a local linear stability analysis,
finding that for specific wall thickness, the vortex shedding identified
behind the wall, can be related with an absolute instability that exists for
some specific outer/inner velocity ratios. The authors explained that this
absolute instability may trigger the destabilization of the flow field. This
theoretical work was verified experimentally by [21]. These authors showed
once more that the wake behind the wall separating the two jets creates a
vortex shedding driving the frequency of the external shear layer also
controlling the evolution of the inner shear layer, which can be the mechanism
that triggers a global absolute instability. This passive mechanism can be
considered as a potential tool for flow control, delaying the transition to
turbulence by means of controlling the near field of the jet. Recently, Canton
et al. (2017) [5] performed a global linear stability analysis to study more
in detail this vortex shedding mechanism behind the wall. They examined a
concentric jet configuration with a very small wall thickness ($0.1D_{i}$,
with $D_{i}$ the inner jet diameter), but the authors selected an outer/inner
velocity ratios where it was known that the alternate vortex shedding behind
the wall was driving the flow. A global unstable mode (absolute instability)
with azimuthal wavenumber $m=0$ was found, confirming that the primary
instability was axisymmetric (the modes with $m=1,2$ were stable at the flow
conditions at which the study was carried out). The highest intensity of the
global mode was located in the wake of the jet, composed by an array of
counter-rotating vortex rings. The shape of the mode changes when moving along
its neutral curve, revealing through the numerical simulations a Kelvin-
Helmholtz instability over the shear-layer between the two jets and in the
outer jet at high Reynolds numbers. Nevertheless, the authors showed that the
wavemaker was located in the bubble formed upstream the separating wall, in
good agreement with the results presented by [32], who performed a similar
stability analysis in a two-dimensional configuration (wakes with co-flow).

The stability of annular jets, a limit case where the inner jets have zero
velocity, has also been investigated. In different analysis of annular jets
[3, 17], it has been illustrated that this type of axisymmetric configuration
does not behave as it appears. The $m=0$ modes studied have been shown to be
stable, and the dominant mode found by both studies is helical ($m=1$). In
addition, to characterise the annular jet, these investigations analyse the
behaviour of the case by adding an azimuthal component to the inflow velocity,
making the discharge of the annular jet eddy-like, comparing the evolution of
the frequency and growth rate of this $m=1$ mode.

This paper expands on the work done by [5], where they use a specific geometry
and vary the outer/inner velocity ratio. This paper presents a complete
characterisation of the main global modes identified in two concentric jets.
The wall thicknesses separating the two jets are defined in the interval
$L\in[0.5,4]$, and the flow is simulated for different inside/outside velocity
ratios in the interval $U_{i}/U_{o}\in[0,2]$, where the case with
$U_{i}/U_{o}=0$ represents an annular jet. Global modes with azimuthal
wavenumber $m=0$ (axisymmetric modes), $m=1$ and $m=2$ will be searched for.
As identified in the literature [17, 3], no axisymmetric modes ($m=0$) could
be identified for any of the distances, as this is a helical case. This paper
expands the conclusions found in these two previous works, extending the
results to different wall thicknesses between jets. This part of the paper
studies in detail the configuration of two concentric jets at low Reynolds
numbers. Using a linear approximation of the equations that model the flow,
the base flows will be obtained on which to apply the linear stability
analysis, by means of which it is possible to identify the most relevant modes
that influence the flow dynamics.

This work also performs a study of mode selection, as some configurations
presents interactions between different modes. Different analysis have been
done to know the different coherent structures when there is an interaction
between modes. [30] conducted the study on the flow past a rotating sphere,
finding different coherent structures on a triple-Hopf bifurcation. Some of
these configurations are steady states, travelling waves or rotating waves.

To the authors’ knowledge, this is the first time that the characterisation of
two concentric jets is presented with such level of detail, presenting neutral
curves for a wide range of different configurations, as well as providing a
deep understanding of the flow physics through the interaction between the
different modes.

The article is organized as follows. Section 2 defines the problem and the
governing equations for the double concentric jets, as well as the linear
stability equations and the methodology for mode selection. The axisymmetric
steaty-state is characterised in Section 3. In Section 4, we perform a
parametric exploration in terms of the velocity ratio between the jets and the
jet distance in order to determine the neutral curves of global stability. The
results about the mode selection are discussed in Section 5. Finally, Section
6 summarises the main conclusions.

##  2 Problem formulation

###  2.1 Computational domain and general equations

The computational domain, represented in fig. 2, models a coaxial flow
configuration, which is composed of two inlet regions, an inner and outer
pipe, both of diameter $D$ and length $5D$, i.e. $z_{min}=-5D$. The
computational domain has an extension of $z_{max}=50D$ and $r_{max}=25D$. The
distance between the pipes is equal to $L$, measured from the inner face of
the outer tube to the face of the inner jet.

Figure 2: Computational domain of the configuration of two concentric jets,
used in StabFem.

The governing equations of the flow within the domain are the incompressible
Navier–Stokes equations. These are written in cylindrical coordinates
$(r,\theta,z)$, which are made dimensionless by considering $D$ as the
reference length scale and $W_{o,max}$ as the reference velocity scale, which
is the maximum velocity in the outer pipe at $z=z_{min}$.

$\displaystyle\frac{\partial\mathbf{U}}{\partial
t}+\mathbf{U}\cdot\nabla\mathbf{U}=-\nabla
P+\nabla\cdot\tau(\mathbf{U}),\qquad\nabla\cdot\mathbf{U}=0,$ (1a)
$\displaystyle\text{with
}\tau(\mathbf{U})=\frac{1}{\mbox{{Re}}}(\nabla\mathbf{U}+\nabla\mathbf{U}^{T}),\qquad\mbox{{Re}}=\frac{W_{o,max}D}{\nu}.$
(1b)

The dimensionless velocity vector $\mathbf{U}=(U,V,W)$ is composed of the
radial, azimuthal and axial components, $P$ is the dimensionless-reduced
pressure, the dynamic viscosity $\nu$ and the viscous stress tensor
$\tau(\mathbf{U})$.

The incompressible Navier–Stokes equations eq. 1 are complemented with the
following boundary conditions

$\mathbf{U}=(0,0,W_{i})\text{ on }\Gamma_{in,i}\text{ and
}\mathbf{U}=(0,0,W_{o})\text{ on }\Gamma_{in,o},$ (2)

where

$W_{i}=\delta_{u}\tanh{\big{(}b_{i}(1-2r)\big{)}}\text{ and
}W_{o}=\tanh\left[b_{o}\left(1+\left\lvert\frac{r-(R_{outer,1}+R_{outer,2})}{D}\right\rvert\right)\right].$

The parameter $\delta_{u}$ corresponds to the velocity ratio between the two
jets, defined as $\delta_{u}=W_{i,max}/W_{o,max}$. The parameters $b_{o}$ and
$b_{i}$ represent the boundary layer thickness within the nozzle, which are
fixed equal to $5$ (as in [5]). There is a weak influence of the boundary
layer thickness on the stability properties of the jet, and it is related to
the vortex shedding regime developed upstream the separation wall (more
details may be found in [31]). Finally, no-slip boundary condition is set on
$\Gamma_{wall}$ and stress-free
($\big{(}\frac{1}{Re}\tau(\mathbf{U})-P\big{)}\cdot\mathbf{n}={\bf 0}$)
boundary condition is set on $\Gamma_{top}$ and $\Gamma_{out}$, as shown in
fig. 2.

In the sequel, Navier–Stokes equations eq. 1 and the associated boundary
conditions will be written symbolically under the form

$\mathbf{B}\frac{\partial\mathbf{Q}}{\partial
t}=\mathbf{F}(\mathbf{Q},\mathbf{\bm{missing}})\equiv\mathbf{L}\mathbf{Q}+\mathbf{N}(\mathbf{Q},\mathbf{Q})+\mathbf{G}(\mathbf{Q},\bm{\eta}),$
(3)

with the flow state vector $\mathbf{Q}=\left[\mathbf{U},{P}\right]^{T}$,
$\bm{\eta}=\left[\mbox{{Re}},\delta_{u}\right]^{T}$. Such a form of the
governing equations takes into account a linear dependency on the state
variable $\mathbf{Q}$ through $\mathbf{L}$. And a quadratic dependency on the
parameters and the state variable through operators $\mathbf{G}(\cdot,\cdot)$
and $\mathbf{N}(\cdot,\cdot)$.

###  2.2 Asymptotic stability

####  2.2.1 Linear stability analysis

In this study, the authors attempt to characterize the stable asymptotic state
from the spectral properties of the Navier–Stokes equations eq. 1. First, let
us consider the stability of an axisymmetric steady-state solution named
$\mathbf{Q}_{0}$, which will be also referred to as trivial steady-state. For
that purpose, let evaluate a solution of eq. 1 in the neighbourhood of the
trivial steady state, i.e., a perturbed state as follows,

$\mathbf{Q}(\mathbf{x},t)=\mathbf{Q}_{0}(\mathbf{x},t)+\varepsilon\hat{\mathbf{q}}(r,z)\text{e}^{-\mathrm{i}(\omega
t-m\theta)}.$ (4)

The next step consists in the characterization of the dynamics of small-
amplitude perturbations around this base flow by expanding them over the basis
of linear eigenmodes (4). If there is a pair
$\left[\mathrm{i}\omega_{\ell},\hat{\mathbf{q}}_{\ell}\right]$ with
$\text{Im}(\omega_{\ell})>0$ (resp. the spectrum is contained in the half of
the complex plane with negative real part) there exists a basin of attraction
in the phase space where the trivial steady-state $\mathbf{Q}_{0}$ is unstable
(resp. stable) [11]. The eigenpair
$\left[\mathrm{i}\omega_{\ell},\hat{\mathbf{q}}_{\ell}\right]$ is determined
as a solution of the following eigenvalue problem,

$\displaystyle\mathbf{J}_{(\omega_{\ell},m_{\ell})}\hat{\mathbf{q}}_{(z_{\ell})}=\Big{(}i\omega_{\ell}\mathbf{B}-\frac{\partial\mathbf{F}}{\partial\mathbf{q}}|_{\mathbf{q}=\mathbf{Q}_{0},\bm{\eta}=\bm{0}}\Big{)}\hat{\mathbf{q}}_{(z_{\ell})},$
(5)

where
$\Big{(}\frac{\partial\mathbf{F}}{\partial\mathbf{q}}|_{\mathbf{q}=\mathbf{Q}_{0},\bm{\eta}=\bm{0}}\Big{)}\hat{\mathbf{q}}_{(z_{\ell})}=\mathbf{L}_{m_{\ell}}\hat{\mathbf{q}}_{(z_{\ell})}+\mathbf{N}_{m_{\ell}}(\mathbf{Q}_{0},\hat{\mathbf{q}}_{(z_{\ell})})+\mathbf{N}_{m_{\ell}}(\hat{\mathbf{q}}_{(z_{\ell})},\mathbf{Q}_{0})$.
The subscript $m_{\ell}$ indicates the azimuthal wavenumber used for the
evaluation of the operator. In the following, we account for eigenmodes
$\hat{\mathbf{q}}_{(z_{\ell})}(r,z)$ that have been normalised in such a way
$\langle\mathbf{\hat{u}}_{(z_{\ell})},\mathbf{\hat{u}}_{(z_{\ell})}\rangle_{L^{2}}=1$.

####  2.2.2 Methodology for the study of mode selection

In the following, we briefly outline the main aspects of the methodology
employed in the study of mode interaction, a comprehensive explanation is left
to appendix A. The determination of the attractor or coherent structure is
explored within the framework of equivariant bifurcation theory. The trivial
steady-state is axisymmetric, i.e. the symmetry group is the orthogonal group
$O(2)$. Near the onset of the instability, dynamics can be reduced to those of
the centre manifold. Particularly, due to the non-uniqueness of the manifold
one can always look for its simplest polynomial expression, which is known as
the normal form of the bifurcation. The reduction to the normal form is
carried out via a multiple scales expansion of the solution $\mathbf{Q}$ of
eq. 3. The expansion considers a two scale development of the original time
$t\mapsto t+\varepsilon^{2}\tau$, here $\varepsilon$ is the order of magnitude
of the flow disturbances, assumed to be small $\varepsilon\ll 1$. In this
study we carry out a normal form reduction via a weakly non-linear expansion,
where the small parameters are

$\varepsilon_{\delta_{u}}^{2}=\delta_{u,c}-\delta_{u}\
\sim\varepsilon^{2}\text{ and
}\varepsilon_{\nu}^{2}=\big{(}\nu_{c}-\nu\big{)}=\big{(}\mbox{{Re}}_{c}^{-1}-\mbox{{Re}}^{-1}\big{)}\sim\varepsilon^{2}.$

A fast timescale $t$ of the self-sustained instability and a slow timescale of
the evolution of the amplitudes $z_{i}(\tau)$ are also considered in eq. 10,
for $i=1,2,3$. The ansatz of the expansion is as follows

$\displaystyle\mathbf{Q}(t,\tau)=\mathbf{Q}_{0}+\varepsilon\mathbf{q}_{(\varepsilon)}(t,\tau)+\varepsilon^{2}\mathbf{q}_{(\varepsilon^{2})}(t,\tau)+O(\varepsilon^{3}).$
(6)

Herein, we evaluate the mode interaction between two steady symmetry breaking
states with azimuthal wave number $m_{1}=1$ and $m_{2}=2$, that is,

$\begin{array}[]{rl}\mathbf{q}_{(\varepsilon)}(t,\tau)&=\big{(}z_{1}(\tau)\hat{\mathbf{q}}_{(z_{1})}(r,z)e^{-\mathrm{i}m_{1}\theta}+\text{c.c.}\big{)}\\\
&+\big{(}z_{2}(\tau)\hat{\mathbf{q}}_{(z_{2})}(r,z)e^{-\mathrm{i}m_{2}\theta}+\text{c.c.}\big{)}.\end{array}$
(7)

Note that the expansion of the LHS of eq. 3 up to third order is as follows

$\displaystyle\varepsilon\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon)}}{\partial
t}+\varepsilon^{2}\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon^{2})}}{\partial
t}+\varepsilon^{3}\big{[}\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon^{3})}}{\partial
t}\big{]}+O(\varepsilon^{4}),$ (8)

and the RHS respectively,

$\mathbf{F}(\mathbf{q},\bm{\eta})=\mathbf{F}_{(0)}+\varepsilon\mathbf{F}_{(\varepsilon)}+\varepsilon^{2}\mathbf{F}_{(\varepsilon^{2})}+\varepsilon^{3}\mathbf{F}_{(\varepsilon^{3})}+O(\varepsilon^{4}).$
(9)

Then, the problem up to third order in $z_{1}$ and $z_{2}$ can be reduced to
[1]

$\begin{array}[]{rl}\dot{z}_{1}&=\lambda_{1}z_{1}+e_{3}\overline{z}_{1}z_{2}+z_{1}\big{(}c_{(1,1)}|z_{1}|^{2}+c_{(1,2)}|z_{2}|^{2}\big{)},\\\
\dot{z}_{2}&=\lambda_{2}z_{2}+e_{4}z^{2}_{1}+z_{2}\big{(}c_{(2,1)}|z_{1}|^{2}+c_{(2,2)}|z_{2}|^{2}\big{)}.\\\
\end{array}$ (10)

An exhaustive analysis of the nonlinear implications of this normal form on
dynamics is left to section 5. The procedure followed for the determination of
the coefficients $c_{(i,j)}$ for $i,j=1,2$ and $e_{3}$ and $e_{4}$ is left to
Appendix A.

####  2.2.3 Numerical methodology for stability tools

Results presented herein follow the same numerical approach adopted by [9, 28,
27, 30], where a comparison with DNS can be found. The calculation of the
steady-state, the eigenvalue problem and the normal form expansion are
implemented in the open-source software FreeFem++. Parametric studies and
generation of figures are collected by StabFem drivers, an open-source project
available in https://gitlab.com/stabfem/StabFem. For steady-state, stability
and normal form computations, we set the stress-free boundary condition at the
outlet, which is the natural boundary condition in the variational
formulation.

The resolution of the steady nonlinear Navier-Stokes equations is tackled by
means of the Newton method. While, the generalised eigenvalue problem (eq. 24)
is solved following the Arnoldi method with spectral transformations. The
normal form reduction procedure of section 2.2.2 only requires to solve a set
of linear systems, which is also carried out within StabFem. On a standard
laptop, every computation considered below can be attained within a few hours.

##  3 Characterisation of the axisymmetric steady-state

###  3.1 Velocity ratio effects

We begin by characterizing the development of the axisymmetric steady-state
with varying $\delta_{u}$ at a constant Reynolds number fixed to
$\text{Re}=100$. Figure 3 synthesises the main topological changes experiences
by the steady-state. At $\delta_{u}=0$, the solution (point (a) in fig. 3)
represents an annular jet, which diffuses as it travels downstream and enters
the ambient fluid. This figure illustrates that the solution curve can be
divided into three segments. The first segment comprised between
$0\leq\delta_{u}<\delta_{u}^{1}$ is characterised by an inner jet nearly
trapped by a large recirculation region with a characteristic length $L_{r}$,
which remains almost constant with the velocity ratio.

Figure 3:  Evolution of the recirculation length ($L_{r}$) of the
recirculating bubble with respect to the velocity ratio $\delta_{u}$ between
the inner and outer jet. The diagram of the second row on the left displays
the minimum value within the domain of the axial velocity. It is spatially
localised within the recirculating region for $\delta_{u}<0.5$ and near the
middle wall for larger values of the velocity ratio. Meridional projections of
the axisymmetric streamfunction isolines and the axial velocity contour in a
range of $(z,r)\in[-1,5]\times[0,5]$.

In the second region, which ranges between
$\delta_{u}^{1}<\delta_{u}<\delta_{u}^{2}$ and it is represented as a shaded
area in the figure, the recirculating region rapidly reduces its size. In this
region, the axial velocity of the inner jet is comparable with the axial
velocity observed in the recirculating region, which promotes mixing between
both regions. As the velocity ratio is increased, the inner jet is
sufficiently energetic to break the recirculating region, which occurs between
point (c) and (d) in fig. 3. The final segment, that ranges between
$\delta_{u}>\delta_{u}^{2}$, is characterised by two quasi-planar jets that
rapidly mix to form a larger one at around $z\approx 5$.

##  4 Linear stability analysis

We explore the parameter space (Re, $\delta_{u}$, $L$). Herein, we examine the
velocity ratio between the jets ($0<\delta_{u}<2$) and the distance between
the jets ($0.5<L<4$). Within this range of parameters, we have analysed the
linear stability properties of the flow configuration. For this purpose, we
first investigate the influence of the jet distance on the stability for the
case of the annular jet ($\delta_{u}=0$).

(a)  (b)  (c)  (d)

(e)  (f)

(g)  (h)

(i)

(j)

Figure 4:  Cross-section view at $z=1$ of the four unstable modes at
criticality for the annular jet case ($\delta_{u}=0$). The streawise component
of the vorticity vector $\varpi_{z}$ is visualised by colours. (a) Mode
$S_{1}$ for $L=0.5$, (b) Mode $S_{2}$ for $L=0.5$, (c) Mode $F_{1}$ for $L=3$
and (d) Mode $F_{2}$ for $L=3$. (e) Linear stability boundaries for the
annular jet ($\delta_{u}=0$). (f) Frequency evolution of the unsteady modes.
Legend: $S_{1}$ mode is displayed with a solid black line, $S_{2}$ with a
solid red line and $F_{1}$ and $F_{2}$ modes are depicted with dashed black
and red lines, respectively. Streamwise velocity of the neutral modes for
$L=3$ and $\delta_{u}=0$ (i) $F_{1}$, (h) $F_{2}$ .

These findings are summarized in fig. 4 which displays the evolution of the
critical Reynolds number with respect to the distance ($L$) for the four most
unstable modes: two steady modes with azimuthal wavenumber $m=1$ and $m=2$,
hereinafter referred to as modes $S_{1}$ and $S_{2}$, respectively. A cross-
section view at $z=1$ is displayed in fig. 4 (a-b). The other two unsteady
modes, named $F_{1}$ and $F_{2}$ have respectively azimuthal wavenumbers $m=1$
and $m=2$. A cross-section view of these two modes is displayed in fig. 4
(c-d). Please note that for the chosen set of parameters the axisymmetric
unsteady mode $F_{0}$, is always found at larger Reynolds numbers than the
aforementioned modes. This is one of the major differences with the case
studied by [5], for small values of the jet distance $L$, the dominant
instability is an unsteady axisymmetric one, which would be named $F_{0}$ with
our nomenclature. Thus, in the following, we only include the results for the
$S_{1}$, $S_{2}$, $F_{1}$ and $F_{2}$ modes. The primary instability of the
annular jet is then a steady symmetry-breaking bifurcation that leads to a jet
flow with a single symmetry plane, displayed in fig. 4 (a). On the contrary,
bifurcations that lead to the mode $S_{2}$ possess two orthogonal symmetry
planes, see fig. 4 (b). As indicated in fig. 4 (g-h), these two stationary
modes $S_{1}$ and $S_{2}$ are localised within the recirculation bubble. For
jet distances $L<2$, the second mode that bifurcates is $F_{1}$ mode, depicted
in fig. 4 (i). This situation corresponds to a bifurcation scenario similar to
other axisymmetric flow configurations, such as the flow past a sphere or a
disk [2, 14]. For larger distances between jets, the scenario changes. The
second bifurcation from the axisymmetric steady-state is the $F_{2}$,
displayed in fig. 4 (j). Other configurations where the primary or secondary
instability involves modes with azimuthal component $m=2$ are swirling jets
[15] and the wake flow past a rotating sphere [29]. The unsteady modes $F_{1}$
and $F_{2}$ possess a much larger spatial support than $S_{1}$ and $S_{2}$.
They are formed by an array of counter-rotating vortex spirals developing in
the wake of the separating duct wall. For the mode $F_{2}$ the amplitude of
these structures grows downstream of the nozzle, in the axial direction, with
a maximum around $z\approx 10$, after which they slowly decay. The mode
$F_{1}$ grows further downstream, with a maximum around $z\approx 50$. The
spatial structure of these eigenmodes resembles the axisymmetric mode of
Figure 9 in [5]. Thus, the steady modes and unsteady modes differ in their
spatial support, that is, even though both steady and unsteady modes are
localised in space, the support of the steady ones is confined within the
recirculation bubble. Instead, the unsteady modes are convected far downstream
until they reach a maximum. This latter characteristic is classical of modes
with a large transient growth, as it was noticed by [5]. On the other hand,
the nature of the steady modes is similar to the symmetry-breaking
instabilities behind the disk or the sphere. These modes are far less
sensitive to transient growth and are observable with direct numerical
simulations and experiments.

###  4.1 Fixed distance between jets and variable velocity ratio $\delta_{u}$

(a)  (b)

(c)  (d)

Figure 5:  Linear stability boundaries for the concentric jets (a) $L=0.5$ and
(b) $L=1$. Same legend as fig. 4. Visualizations of real part of the
streamwise axial velocity of the critical modes (c) $S_{1}$ and (d) $S_{2}$.

In the following, we focus on the influence of the velocity ratio $\delta_{u}$
between jets for fixed jet distances $L$. Figure 5 displays the neutral curve
of stability for jet distances (a) $L=0.5$ and (b) $L=1$. One may observe that
the primary bifurcation is not always associated to the mode $S_{1}$ as it is
the case for $\delta_{u}=0$. For sufficiently large velocity ratios, the
primary instability leads to a non-axisymmetric steady-state with a double
helix. Another interesting feature, which could motivate a control strategy,
is the occurrence of vertical asymptotes. This sudden change in the critical
Reynolds number is due to the retraction and disappearance of the
recirculating region. For $L=0.5$, this sudden change occurs for
$\delta_{u}\approx 0.25$, and for higher values of $\delta_{u}$ the critical
Reynolds number is around twice larger than the one of the annular jet
($\delta_{u}=0$). The case of jet distance $L=1$ was discussed in section 3.
The sudden change in the stability of the branch $S_{1}$ occurs between
$\delta_{u}\in[0.25,0.5]$. Within this narrow interval, the primary branch of
instability is the $F_{1}$. At around $\delta_{u}=0.4$, the primary
bifurcation is again the branch $F_{1}$, which becomes secondary at around
$\delta_{u}\approx 0.8$ in favour of the branch $S_{2}$. In fig. 5 we have
highlighted the codimension two point interaction between the $S_{1}-S_{2}$
modes, whose modes are depicted in fig. 5 (c-d), which will be analysed in
detail in section 5. Around this point, we can observe the largest
stabilisation ratio between the annular jet ($\delta_{u}=0$) and a
configuration of concentric jets ($\delta_{u}\neq 0$).

###  4.2 Fixed velocity ratio $\delta_{u}$ and variable distance between jets

a) b)

c) d)

e) f)

Figure 6: Neutral lines of the four modes found studying the configuration of
two concentric jets fixing the velocity ratio. (a-b) $\delta_{u}=0.5$, (c-d)
$\delta_{u}=1$, (e-f) $\delta_{u}=2$. Black lines: modes with $m=1$, red
lines: modes with $m=2$. Straight lines: steady modes, dashed lines: unsteady
modes.

Figure 6 compares the results obtained for a constant velocity ratio when
varying the distance between jets. As observed before, the solution becomes
more unstable by increasing the distance between jets. The largest critical
Reynolds number is found at $\delta_{u}=0$. The critical Reynolds decreases
with the distance between jets $L$. The points of codimension two, i.e., the
points where mode switching occurs, are highlighted in fig. 6. We can
appreciate that the interaction between the branch $S_{1}$ and $S_{2}$ happens
for every velocity ratio $\delta_{u}$ explored, and it is the mode interaction
associated to the smallest distance between jets. Additionally, for a velocity
ratio $\delta_{u}=0.5$ there exists two points where the branches of the
linear modes $S_{1}$ and $F_{1}$ intersect. Another feature of the neutral
curves is the existence of turning points, which are associated to the
existence of saddle node bifurcations of the axisymmetric steady-state. The
saddle node bifurcations of the steady-state induces the existence of regions
in the neutral curves with a tongue shape. These saddle node bifurcations are
also responsible for the formation of the vertical asymptotes observed in fig.
5. Finally, it is of interest the transition of the modes $S_{1}$ and $S_{2}$,
which induce the symmetry breaking of the axisymmetric steady state to slow
low frequency spiralling structures. These modes have been identified for
$\delta_{u}=0.5$ for $m=1$, $\delta_{u}=1$ for $m=2$, and $\delta_{u}=2$ for
both $m=1$ and $m=2$. As it will be clarified in section 5, these oscillations
are issued from the non-linear interaction of modes, emerging simultaneously
for a specific Reynolds number, and changing their position as the most
unstable global mode of the flow.

##  5 Mode interaction between two steady states. Resonance $1:2$

###  5.1 Normal form, basic solutions and their properties

The linear diagrams of section 4 have shown the existence of the mode
interaction between the modes $S_{1}$ and $S_{2}$. It corresponds roughly to
the mode interaction that occurs at the largest critical Reynolds number for
any value of $L$ herein explored. In this section, we analyse the dynamics
near the $S_{1}:S_{2}$ organizing centre. We perform a normal form reduction,
which allows us to predict non-axisymmetric steady, periodic, quasiperiodic
and heteroclinic cycles between non-axisymmetric states.
The mode interaction that is herein analysed corresponds to a steady-steady
bifurcation with $O(2)$ symmetry and with strong resonance $1:2$. Such a
bifurcation scenario has been extensively studied in the past by [8, 10, 23,
1] and the reflection symmetry breaking case (SO(2)) by [24]. In order to
unravel the existence and the stability of the nonlinear states near the
codimension two point, let write the flow field as

$\begin{array}[]{rl}\displaystyle\mathbf{q}&=\mathbf{Q}_{0}+\text{Re}\big{[}r_{1}(\tau)e^{i\phi_{1}(\tau)}e^{-i\theta}\hat{\mathbf{q}}_{s,1}\big{]}+\text{Re}\big{[}r_{2}(\tau)e^{i\phi_{2}(\tau)}e^{-2i\theta}\hat{\mathbf{q}}_{s,2}\big{]}\\\
\end{array}$ (11)

in polar coordinates for the complex amplitudes
$z_{1}=r_{1}\text{e}^{\mathrm{i}\phi_{1}}$ and
$z_{2}=r_{2}\text{e}^{\mathrm{i}\phi_{2}}$ where $r_{j}$ and $\phi_{j}$ for
$j=1,2$ are the amplitude and phase of the symmetry-breaking modes $m=1$ and
$m=2$, respectively. The complex-amplitude normal form eq. 10 is expressed in
this reduced polar notation as follows,

$\displaystyle\begin{split}\dot{r}_{1}=e_{3}r_{1}r_{2}\cos(\chi)+r_{1}\Big{(}\lambda_{(s,1)}+c_{(1,1)}r_{1}^{2}+c_{(1,2)}r_{2}^{2}\Big{)},\end{split}$
(12a)
$\displaystyle\begin{split}\dot{r}_{2}=e_{4}r_{1}^{2}\cos(\chi)+r_{2}\Big{(}\lambda_{(s,2)}+c_{(2,1)}r_{1}^{2}+c_{(2,2)}r_{2}^{2}\Big{)},\end{split}$
(12b)
$\displaystyle\begin{split}\dot{\chi}=-\Big{(}2e_{3}r_{2}+e_{4}\frac{r_{1}^{2}}{r_{2}}\Big{)}\sin(\chi),\end{split}$
(12c)

where the phase $\chi=\phi_{2}-2\phi_{1}$ is coupled with the amplitudes
$r_{1}$ and $r_{2}$ because of the existence of the $1:2$ resonance. The
individual phases evolve as

$\begin{array}[]{rl}\dot{\phi}_{1}&=e_{3}r_{2}\sin(\chi),\\\
\dot{\phi}_{2}&=-e_{4}\frac{r_{1}^{2}}{r_{2}}\sin(\chi).\end{array}$ (13)

Before proceeding to the analysis of the basic solutions of eq. 12, we can
simplify these equations by the rescaling

$\big{(}\frac{r_{1}}{|e_{3}e_{4}|^{1/2}},\frac{r_{2}}{e_{3}}\big{)}\to(r_{1},r_{2}),$

which yields the following equivalent system

$\displaystyle\begin{split}\dot{r}_{1}=r_{1}r_{2}\cos(\chi)+r_{1}\Big{(}\lambda_{(s,1)}+c_{11}r_{1}^{2}+c_{12}r_{2}^{2}\Big{)},\end{split}$
(14a)
$\displaystyle\begin{split}\dot{r}_{2}=sr_{1}^{2}\cos(\chi)+r_{2}\Big{(}\lambda_{(s,2)}+c_{21}r_{1}^{2}+c_{22}r_{2}^{2}\Big{)},\end{split}$
(14b)
$\displaystyle\begin{split}\dot{\chi}=-\frac{1}{r_{2}}\Big{(}2r_{2}^{2}+sr_{1}^{2}\Big{)}\sin(\chi),\end{split}$
(14c)

where the coefficients

$s=\text{sign}(e_{3}e_{4}),\quad c_{11}=\frac{c_{(1,1)}}{|e_{3}e_{4}|},\quad
c_{12}=\frac{c_{(1,2)}}{e_{3}^{2}},\quad
c_{21}=\frac{c_{(2,1)}}{|e_{3}e_{4}|},\quad
c_{22}=\frac{c_{(2,2)}}{e_{3}^{2}}.$

Finally, we consider a third normal form equivalent to the previous ones but
which removes the singularity of eqs. 12 and 14 when $r_{2}=0$. Standing waves
($\sin{\chi}=0$) naturally encounter this type of artificial singularity,
which manifests as in eq. 14 as an instantaneous jump from one standing
subspace to the other by a $\pi$-translation. This is the case of the
heteroclinic cycles, previously studied by [1, 23]. The third normal form,
which we shall refer to as reduced Cartesian normal form, takes advantage of
the simple transformation $x=r_{2}\cos(\chi)$, $y=r_{2}\sin(\chi)$ [24]:

$\displaystyle\begin{split}\dot{r}_{1}=r_{1}\Big{(}\lambda_{(s,1)}+c_{11}r_{1}^{2}+c_{12}(x^{2}+y^{2})+x\Big{)},\end{split}$
(15a)
$\displaystyle\begin{split}\dot{x}=sr_{1}^{2}+2y^{2}+x\Big{(}\lambda_{(s,2)}+c_{21}r_{1}^{2}+c_{22}(x^{2}+y^{2})\Big{)},\end{split}$
(15b)
$\displaystyle\begin{split}\dot{y}=-2xy+y\Big{(}\lambda_{(s,2)}+c_{21}r_{1}^{2}+c_{22}(x^{2}+y^{2})\Big{)},\end{split}$
(15c)

In this final representation standing wave solutions are contained within the
invariant plane $y=0$, and due to the invariance of eq. 15 under the
reflection $y\mapsto-y$, one can restrict attention, without loss of
generality, to solutions with $y\geq 0$, cf [23].

The system eq. 14 possess four types of fixed points, which are listed in
table 1.

Name | Definition | Bifurcations | Comments
---|---|---|---
O | $r_{1,O}=r_{2,O}=0$ | $-$ | Steady axisymmetric state
P |  $r_{2,P}^{2}=\frac{-\lambda_{(s,2)}}{c_{22}}$, $r_{1,P}=0$ | $\lambda_{(s,2)}=0$ | Bifurcation from O
| $r_{1,MM}=-\frac{\lambda_{(s,1)}\pm r_{2,MM}+c_{12}r_{2,MM}^{2}}{c_{11}}$ | $\lambda_{(s,1)}=0$ | Bifurcation from O
MM | $\text{P}_{\text{MM}}(r_{2,MM}\cos(\chi_{MM}))=0$ | $\sigma_{\pm}=0$ | Bifurcation from P
| $\cos(\chi_{MM})=\pm 1$ |  |
| $\cos(\chi_{TW})=\frac{(2c_{11}+c_{12})\lambda_{(s,2)}-(2c_{21}+c_{22})\lambda_{(s,1)}}{\Sigma_{TW}(2\lambda_{(s,1)}+\lambda_{(s,2)})}$ |  |
TW | $r_{2,TW}^{2}=\frac{-(2\lambda_{(s,1)}+\lambda_{(s,2)})}{\Sigma_{TW}}$ | $\cos(\chi_{TW})=\pm 1$ | Bifurcation from MM
| $r_{1,TW}^{2}=2r^{2}_{2,TW}$ |  |
Table 1: Definition of the fixed points of the reduced polar normal form eq.
14. $\sigma_{\pm}$ is defined in eq. 16, the polynomial $P_{MM}$ is defined in
eq. 17 and $\Sigma_{TW}\equiv 4c_{11}+2(c_{12}+c_{21})+c_{22}$.

First, the axisymmetric steady state (O) is represented by
$(r_{1},r_{2})=(0,0)$, so it is the trivial steady-state of the normal form.
The second steady-state is what it is denoted as pure mode (P). In the
original coordinates, it corresponds to the symmetry breaking structure
associated to the mode $S_{2}$. This state bifurcates from the axisymmetric
steady state (O) when $\lambda_{(s,2)}=0$. The third fixed point is the mixed
mode state (MM), which is listed in table 1. It corresponds to the reflection
symmetry preserving state associated to the mode $S_{1}$. It may bifurcate
directly from the trivial steady state O, when $\lambda_{(s,1)}=0$ or from P
whenever $\sigma_{+}=0$ or $\sigma_{-}=0$, where $\sigma_{\pm}$ is defined as

$\sigma_{\pm}\equiv\lambda_{(s,1)}-\frac{-\lambda_{(s,2)}c_{12}}{c_{22}}\pm\sqrt{\frac{-\lambda_{(s,2)}}{c_{22}}}.$
(16)

The representation in the reduced polar form is

$r_{1,MM}=-\frac{\lambda_{(s,1)}\pm
r_{2,MM}+c_{12}r_{2,MM}^{2}}{c_{11}},\qquad\cos(\chi_{MM})=\pm 1,$

and the condition $\text{P}_{\text{MM}}(r_{2,MM}\cos(\chi_{MM}))=0$, where
$\text{P}_{\text{MM}}$ is defined as

$\text{P}_{\text{MM}}(x)\equiv
s\mu_{1}+(s+c_{21}\lambda_{(s,1)}-c_{(1,1)}\lambda_{(s,2)})x+(c_{21}+sc_{12})x^{2}+(c_{12}c_{21}-c_{11}c_{22})x^{3}.$
(17)

Finally, the fourth fixed point of the system are travelling waves (TW). It is
surprising that the interaction between two steady-states causes a time-
periodic solution. The travelling wave emerges from MM in parity-breaking
pitchfork bifurcation that breaks the reflection symmetry when
$\cos(\chi_{TW})=\pm 1$. The TW drifts at a steady rotation rate $\omega_{TW}$
along the group orbit, i.e., the phases
$\dot{\phi}_{1}=r_{2,TW}\sin(\chi_{TW})$ and
$\dot{\phi}_{2}=-s\frac{r_{1,TW}^{2}}{r_{2,TW}}\sin(\chi_{TW})$ are non-null.

Name | Bifurcation condition | Comments
---|---|---
SW | $sr_{1}^{2}-2c_{11}r_{1}^{2}r_{2,MM}\cos(\chi_{MM})-2c_{22}r_{2,MM}^{3}\cos(\chi_{MM})^{3}=0$ | Bif. from MM
MTW |  $D_{TW}-T_{TW}I_{TW}=0$, $I_{TW}>0$ | Bif. from TW
Table 2: Definition of the limit cycles of the reduced polar normal form eq.
14.

Mixed modes and travelling waves may further bifurcate into standing waves
(SW) and modulated travelling waves (MTW), respectively. These are generic
features of the $1:2$ resonance for small values of $\lambda_{(s,1)}$ and
$\lambda_{(s,2)}$, when $s=-1$. In the original coordinates, SW are periodic
solutions, whereas MTW are quasiperiodic. Standing waves emerge via a Hopf
bifurcation from MM when the conditions
$\text{P}_{\text{SW}}\Big{(}r_{2,MM}\cos(\chi_{MM})\Big{)}>0$ for

$\text{P}_{\text{SW}}(x)\equiv(2c_{22}x^{3}-sr_{1}^{2})c_{11}-(2c_{12}x+1)(c_{21}x+s)x,$

and the one listed in table 2 are satisfied. MTW are created when a torus
bifurcation happens on the travelling wave branch when the conditions listed
in table 2 are satisfied.

Name | Condition | Comments
---|---|---
Ht AGH |  $\lambda_{(s,1)}>0$, $\lambda_{(s,2)}>0$, $c_{22}<0$ | Existence
| $\sigma_{+}>0,\sigma_{-}<0$ | Asymptotic stability
Table 3: Definition of the conditions for the existence of the Ht AGH (robust
heteroclinic cycles connecting pure modes) of the reduced polar normal form
eq. 14.  Figure 7:  Schematic representation of the basic solutions of eq. 12
and their bifurcation path.

Another remarkable feature of eq. 12 is the existence of robust heteroclinic
cycles that are asymptotically stable. When $s=-1$, there are open sets of
parameters (see table 3) where the reduced polar normal form exhibits
structurally stable connections between $\pi-$translations on the circle of
pure modes, cf [1]. These structures are robust and have been observed in a
large variety of systems, [19, 18, 16, 22, 13]. In addition to these robust
heteroclinic cycles connecting pure modes, there exists more complex limit
cycles connecting O, P, MM and SW, cf [23]. These cycles are located for
larger values of $\lambda_{(s,1)}$ and $\lambda_{(s,2)}$, with possibly
chaotic dynamics (Shilnikov type). In this study, we have not identified any
of these. Finally, a summary of the basic solutions and the bifurcation path
is sketched in fig. 7.

###  5.2 Results of the steady-steady $1:2$ mode interaction

Section 4.2 reported the location of mode interaction points for discrete
values of the velocity ratio $\delta_{u}$. The location of the mode
interaction between $S_{1}$ and $S_{2}$ is depicted in fig. 8. It shows that
the mode switching between the modes $S_{1}$ and $S_{2}$ is indeed stationary
only for $\delta_{u}<1.5$ and $L<1.3$. For larger values of the velocity ratio
and the jet distance, the interaction is not purely stationary; at least one
of the linear modes oscillates with a slow frequency. It implies that the mode
selection for large velocity ratios near the codimension two points is similar
to the one reported by [15] for swirling jets. However, even when the two
primary bifurcations are non-oscillating ($S_{1}$ and $S_{2}$), the $1:2$
resonance of the azimuthal wavenumbers induces a slow frequency, what we
denote as travelling wave solutions (TW).

a) b)

Figure 8:  Evolution of the codimension two interaction $S_{1}-S_{2}$ in the
space of parameters $(\text{Re},L,\delta_{u})$. Grey points denote the points
that were computed and the red point denotes the transition from steady to
unsteady with low frequency as reported in section 4.2.  Figure 9:  Phase
portrait at the codimension two point $S_{1}:S_{2}$ for parameter values
$(L,\delta_{u})=(1.15,1.0)$. Visualisations of blue and red surfaces in the
isometric views represent the respective positive and negative isocontour
values of the perturbative axial velocity indicated in the figure.

Figure 10:  Bifurcation diagram with respect to the Reynolds number for
$L=1.15$ and $\delta_{u}=0.8$. The left diagram reports the evolution of
$r_{2}$ for the fixed point solutions of the normal form. The right diagram
displays the bifurcation diagram in the Cartesian coordinates. Solid lines and
dashed lines denote stable attractors and unstable attractors, respectively.

We consider the bifurcation sequence for $\delta_{u}=1.0$ and $L=1.15$, which
is qualitatively similar to transitions in the range $0.5<\delta_{u}<1.5$,
near the codimension two points, which are depicted in fig. 8. At the
codimension two points for $\delta_{u}<0.5$, at least one of the two
bifurcations is sub-critical and a normal form reduction up to fifth order is
necessary. Subcritical transition was also noticed for a distance between jets
$L=0.1$ by [5], who reported high levels of the linear gain associated to
transient growth mechanisms. This last case is out of the scope of the present
manuscript. Figure 9 displays the phase portrait of the stable attractors near
the $S_{1}:S_{2}$ interaction. For values of $\delta_{u}>1.0$, the
axisymmetric steady-state loses its axisymmetry leading to a new steady-state
with symmetry $m=2$, herein denoted as pure mode (P). A reconstruction of the
fluctuating flow field of such a state is performed at the bottom right of
fig. 9, which shows that the state P possesses two orthogonal planes of
symmetry. Near the codimension two point, for values of the velocity ratio
$\delta_{u}<1.1$, the state P is only observable, that is non-linearly stable,
within a small interval with respect to the Reynolds number. For larger values
of the velocity ratio, the state P remains stable within the analysed interval
of Reynolds numbers. For values of the velocity ratio $\delta_{u}<1.0$, the
bifurcation diagram is more complex. Figure 10 displays the bifurcation
diagram of the fixed-point solutions of eq. 15 on the left diagram and the
full set of solutions of the normal form in the right diagram. The
axisymmetric steady-state first bifurcates towards a Mixed-Mode solution,
which is the solution in the $y=0$ plane for the right diagram of fig. 10. A
solution with a non-symmetric wake has been reconstructed in fig. 9. The
Mixed-Mode solution is only stable within a small interval of the Reynolds
number. A secondary bifurcation, denoted BifMM-TW, gives raise to a slowly
rotating wave of the wake. The TW and the MM solutions are identical at the
bifurcation point. The phase speed is zero at the bifurcation, thus this is
not a Hopf bifurcation. It corresponds to a drift instability that breaks the
azimuthal symmetry, i.e. it starts to slowly drift. This unusual feature, that
travelling waves bifurcate from a steady solution at a steady bifurcation, is
a generic feature of the $1:2$ resonance. A reconstruction of the travelling
wave solution is depicted on the top of fig. 9. It corresponds to the line
with non-zero $y$ component in the right diagram of fig. 10. The TW solution
loses its stability in a tertiary bifurcation, denoted as BifTW-MTW. It
conforms to a Hopf bifurcation of the TW solution, which gives birth to a
quasi-periodic solution name Modulated Travelling Wave (MTW). A representation
of this kind of solution in the Cartesian coordinates $(r_{1},x,y)$ is
depicted on the right image of fig. 10.

Figure 11:  Heteroclinic cycle solution for parameter values $\text{Re}=200$,
$\delta_{u}=0.8$. The top and bottom image sequences along the heteroclinic
cycle show (from left to right) an axial slice plane at $z=1$ of the
instantaneous fluctuations of the axial velocity of the flow field as viewed
from downstream, along with a three-dimensional isometric view (d on the top
and g on the bottom).The middle diagram displays the heteroclinic cycle in the
coordinates ($r_{1},x,y$).

Eventually, the Modulated Travelling Wave experiences a global bifurcation.
That occurs when the periodic MTW solution, in the $(r_{1},x,y)$ coordinates,
nearly intersects the invariant $r_{1}=0$ and $y=0$ planes. The transition
sequence is represented in the right image of fig. 10 in the Cartesian
coordinates ($r_{1},x,y$). The amplitude of the MTW limit cycle increases
until the MTW arising at the tertiary bifurcation BifTW-MTW are destroyed by
meeting a heteroclinic cycle at BifMTW-Ht. The locus of BifMTW-Ht is reported
in fig. 9 and in good agreement with [1]. The conditions for the existence of
the heteroclinic cycles are listed in table 3. When $\sigma_{-}$ becomes
negative, the cycle is attracting and robust heteroclinic cycles are observed.
It is destroyed when $\sigma_{+}$ becomes negative, in that case the pure
modes are no longer saddles which breaks the heteroclinic connection. Figure
11 displays the instantaneous fluctuation field from a heteroclinic orbit
connecting P and its conjugate solution P’, which is obtained by a rotation of
$\pi/2$, for parameter values $\text{Re}=200$ and $\delta_{u}=0.8$. The
dynamics of the cycle takes place in two phases. Figure 11 depicts the motion
of the coherent structure associated to the heteroclitic cycle. Starting from
the conjugated pure mode P’, the cycle leaves the point (a), located in the
vicinity of P’, along the unstable eigenvector $y$, which is the stable
direction of P. The first phase consists in a rapid rotation by $\pi/2$ of the
wake, it corresponds to the sequence a-b-c-d-e displayed in fig. 11. Then it
is followed by a slow approach following the direction $y$ and departure from
the pure mode state P along the direction $r_{1}$. The second phase consists
in a rapid horizontal motion of the wake, which is an evolution from P to P’
that takes place by the breaking of the reflectional symmetry with respect to
the vertical axis; it constitutes the sequence e-f-g-h-i-a. Please note that
equivalent motions are also possible. The first phase of rapid counter-
clockwise rotation by $\pi/2$ can be performed in the opposite sense. It
corresponds to a motion in the Cartesian coordinates along the plane $r_{1}$
along negative values of $y$. The sequence e-f-g-h-i-a can be replaced by a
horizontal movement in the opposite sense, which adjusts to connect the plane
$y=0$ corresponding to negative values of $r_{1}$,

##  6 Discussion & Conclusions

This article achieves a complete description of the configuration consisting
of two coaxial jets, broadly found in industrial processes, covering a wide
range of applications such as noise reduction, mixing enhancement, or
combustion control. The numerical procedure herein employed has been validated
with the existing literature in the case of the stability analysis (see B for
a detailed overview), and compared to DNS results, as done in [30]. The
analysis comprises a layout with a wide range of the velocity ratio
($\delta_{u}=U_{i}/U_{o}$) between the jets, from $\delta_{u}=0$ to
$\delta_{u}=2$, as well as the distance between jets ($L$) enclosing values
from $L=0.5$ to $L=4$, substantially expanding the work of [5].

A linear stability analysis reveals the most significant modes, consisting of
two steady modes ($S_{1}$ and $S_{2}$, located within the recirculation
bubble) and two unsteady ones ($F_{1}$ and $F_{2}$, evolving as a transient
growth in the downstream direction). The critical Reynolds number is
determined for a wide range of velocity and distance ratios, starting with the
influence of the velocity ratio. As the relation between inner and outer
velocities grows, the flow is stabilised, increasing the critical Reynolds
number. The primary instability swaps from mode $S_{1}$, characterised with
one symmetry plane, to mode $S_{2}$ that possesses two symmetry planes. An
abrupt divergence in the critical Reynolds number is captured, associated with
the vanishing of the recirculation region, that could suggest a stability
control strategy. Subsequently, the effect of the distance $L$ between jets is
analysed. The primary effect of increasing this distance is a decrease in the
critical Reynolds number for all values of $\delta_{u}$ investigated.
Additionally, the existence of saddle node bifurcations, that swap the most
unstable mode of the flow, generates turning points in the neutral curve.

The investigated bifurcation scenario starts from the codimension point, with
an axisymmetric steady state located at a velocity ratio $\delta_{u}=1.0$ and
distance between jets of $L=1.15$. It is qualitatively equivalent to
transitions found in the range $0.5<\delta_{u}<1.5$. It reveals a break of the
axisymmetry for values higher than $\delta_{u}=1.0$, presenting a steady state
as a pure mode P with two orthogonal planes of symmetry. For values lower than
$\delta_{u}=1.0$, the bifurcation diagram exhibits a slightly complicated
path. Firstly, it drives into a Mixed-Mode (MM) solution presenting a non-
symmetric wake, that is only stable for a small range of the Reynolds number.
Subsequently, a slowly rotating wake is triggered in the form of a Travelling
Wave (TW). This unusual feature, an unsteady state emerging from a steady
state, corresponds to a drift instability commonly found at $1:2$ resonance.
Then, the TW solution encounters a Hopf bifurcation, developing a quasi-
periodic solution in the form of a Modulated Travelling Wave (MTW). Finally,
the MTW solution undergoes a global bifurcation meeting a heteroclinic cycle
(Ht). This heteroclinic orbit links the solution P with its conjugate solution
P’, spinning the wake from P’ to P, and moving it horizontally from P to P’.

## Acknowledgments

A.C., J.A.M. and S.L.C. acknowledge the grant PID2020-114173RB-I00 funded by
MCIN/AEI/ 10.13039/501100011033. S.L.C., J.A.M. and A.C. acknowledge the
support of Comunidad de Madrid through the call Research Grants for Young
Investigators from Universidad Politécnica de Madrid. AC also acknowledges the
support of Universidad Politécnica de Madrid, under the programme ‘Programa
Propio’.

##  Appendix A Normal form reduction

In this section we provide a detail explanation of the normal form reduction
to obtain the coefficients of eq. 10, we define the terms of the compact
notation of the governing equations eq. 3, which is reminded here, for the
sake of conciseness,

$\mathbf{B}\frac{\partial\mathbf{Q}}{\partial
t}=\mathbf{F}(\mathbf{Q},\bm{\eta})\equiv\mathbf{L}\mathbf{Q}+\mathbf{N}(\mathbf{Q},\mathbf{Q})+\mathbf{G}(\mathbf{Q},\bm{\eta}).$
(18)

The nonlinear convective operator
$\mathbf{N}(\mathbf{Q}_{1},\mathbf{Q}_{2})=\mathbf{U}_{1}\cdot\nabla\mathbf{U}_{2}$
accounts for the quadratic interaction on the state variable. The linear
operator on the state variable is $\mathbf{L}\mathbf{Q}=[\nabla
P,\nabla\cdot\mathbf{U}]^{T}$. The remaining term accounts for the linear
variations in the state variable and the parameter vector. It is defined as
$\mathbf{G}(\mathbf{Q},\bm{\eta})=\mathbf{G}(\mathbf{Q},[\eta_{1},0]^{T})+\mathbf{G}(\mathbf{Q},[0,\eta_{2}]^{T})$
where
$\mathbf{G}(\mathbf{Q},[\eta_{1},0]^{T})=\eta_{1}\nabla\cdot(\nabla\mathbf{U}+\nabla\mathbf{U}^{T})$
and $\mathbf{G}(\mathbf{Q},[0,\eta_{2}]^{T})$. The former operator shows the
dependency on the parameter $\eta_{1}$, which accounts for the viscous
effects. The latter operator depends on the parameter $\eta_{2}$, which
accounts for the velocity ratio between jets and it is used to impose the
boundary condition
$\mathbf{U}=(0,\eta_{2}\tanh{\big{(}b_{i}(1-2r)\big{)}},0)\text{ on
}\Gamma_{in,i}$. In addition, we consider the following splitting of the
parameters $\bm{\eta}=\bm{\eta}_{c}+\Delta\bm{\eta}$. Here $\bm{\eta}_{c}$
denotes the critical parameters
$\bm{\eta}_{c}\equiv[Re_{c}^{-1},\delta_{u,c}]^{T}$ attained when the spectra
of the Jacobian operator posses at least an eigenvalue whose real part is
zero. The distance in the parameter space to the threshold is represented by
$\Delta\bm{\eta}=[Re_{c}^{-1}-Re^{-1},\delta_{u,c}-\delta_{u}]^{T}$.

###  A.1 Multiple scales ansatz

The multiple scales expansion of the solution $\mathbf{q}$ of eq. 3 is

$\displaystyle\mathbf{q}(t,\tau)=\mathbf{Q}_{0}+\varepsilon\mathbf{q}_{(\varepsilon)}(t,\tau)+\varepsilon^{2}\mathbf{q}_{(\varepsilon^{2})}(t,\tau)+O(\varepsilon^{3}),$
(19)

where $\varepsilon\ll 1$ is a small parameter. The distance in the parameter
space to the critical point
$\Delta\bm{\eta}=[Re_{c}^{-1}-Re^{-1},\delta_{u,c}-\delta_{u}]^{T}$ is assumed
to be of second order, i.e. $\Delta{\eta}_{i}=O(\varepsilon^{2})$ for $i=1,2$.
The expansion eq. 19 considers a two scale expansion of the original time
$t\mapsto t+\varepsilon^{2}\tau$. A fast timescale $t$ and a slow timescale of
the evolution of the amplitudes $z_{i}(\tau)$ in eq. 19, for $i=1,2$. Note
that the expansion of the LHS eq. 3 up to third order is as follows

$\displaystyle\varepsilon\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon)}}{\partial
t}+\varepsilon^{2}\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon^{2})}}{\partial
t}+\varepsilon^{3}\big{[}\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon^{3})}}{\partial
t}+\mathbf{B}\frac{\partial\mathbf{q}_{(\varepsilon)}}{\partial\tau}\big{]},$
(20)

and the RHS respectively,

$\mathbf{F}(\mathbf{q},\bm{\eta})=\mathbf{F}_{(0)}+\varepsilon\mathbf{F}_{(\varepsilon)}+\varepsilon^{2}\mathbf{F}_{(\varepsilon^{2})}+\varepsilon^{3}\mathbf{F}_{(\varepsilon^{3})}.$
(21)

The expansion eq. 21 will be detailed at each order.

####  A.1.1 Order $\varepsilon^{0}$

The zeroth order $\mathbf{Q}_{0}$ of the multiple scales expansion eq. 19 is
the steady state of the governing equations evaluated at the threshold of
instability, i.e. $\bm{\eta}=\bm{\eta}_{c}$,

$\displaystyle\mathbf{0}=\mathbf{F}(\mathbf{Q}_{0},\bm{\eta}_{c}).$ (22)

####  A.1.2 Order $\varepsilon^{1}$

The first order $\mathbf{q}_{(\varepsilon)}(t,\tau)$ of the multiple scales
expansion of eq. 19 is composed of the eigenmodes of the linearized system

$\mathbf{q}_{(\varepsilon)}(t,\tau)\equiv\big{(}z_{1}(\tau)e^{-\mathrm{i}m_{1}\theta}\hat{\mathbf{q}}_{1}+z_{2}(\tau)e^{\mathrm{i}-m_{2}\theta}\hat{\mathbf{q}}_{2}+\text{c.
c.}\big{)}.$ (23)

in our case, $m_{1}=1$ and $m_{2}=2$. Each term $\hat{\mathbf{q}}_{\ell}$ of
the first order expansion eq. 23 is a solution of the following linear
equation

$\mathbf{J}_{(\omega_{\ell},m_{\ell})}\hat{\mathbf{q}}_{\ell}=\Big{(}i\omega_{\ell}\mathbf{B}-\frac{\partial\mathbf{F}}{\partial\mathbf{q}}|_{\mathbf{q}=\mathbf{Q}_{0},\bm{\eta}=\bm{\eta}_{c}}\Big{)}\hat{\mathbf{q}}_{\ell},$
(24)

where
$\frac{\partial\mathbf{F}}{\partial\mathbf{q}}|_{\mathbf{q}=\mathbf{Q}_{0},\bm{\eta}=\bm{\eta}_{c}}\hat{\mathbf{q}}_{\ell}=\mathbf{L}_{m_{\ell}}\hat{\mathbf{q}}_{\ell}+\mathbf{N}_{m_{\ell}}(\mathbf{Q}_{0},\hat{\mathbf{q}}_{\ell})+\mathbf{N}_{m_{\ell}}(\hat{\mathbf{q}}_{\ell},\mathbf{Q}_{0})$.
The subscript $m_{\ell}$ indicates the azimuthal wavenumber used for the
evaluation of the operator.

####  A.1.3 Order $\varepsilon^{2}$

The second order expansion term $\mathbf{q}_{(\varepsilon^{2})}(t,\tau)$ is
determined from the resolution of a set of forced linear systems, where the
forcing terms are evaluated from first and zeroth order terms. The expansion
in terms of amplitudes $z_{i}(\tau)$ ($i=1,2$) of
$\mathbf{q}_{(\varepsilon^{2})}(t,\tau)$ is assessed from term-by-term
identification of the forcing terms at the second order. Non-linear second
order terms in $\varepsilon$ are

$\begin{array}[]{rl}\displaystyle\mathbf{F}_{(\varepsilon^{2})}&\equiv\displaystyle\sum_{j,k=1}^{2}\Big{(}z_{j}z_{k}\mathbf{N}(\hat{\mathbf{q}}_{j},\hat{\mathbf{q}}_{k})e^{-\mathrm{i}(m_{j}+m_{k})\theta}+\text{
c.c.}\Big{)}\\\
&+\displaystyle\sum_{j,k=1}^{2}\Big{(}z_{j}\overline{z}_{k}\mathbf{N}(\hat{\mathbf{q}}_{j},\overline{\hat{\mathbf{q}}}_{k})e^{-\mathrm{i}(m_{j}-m_{k})\theta}+\text{
c.c.}\Big{)}\\\
&+\displaystyle\sum_{\ell=0}^{2}\eta_{\ell}\mathbf{G}(\mathbf{Q}_{0},\mathbf{e}_{\ell}),\end{array}$
(25)

where the terms proportional to $z_{j}z_{k}$ are named
$\hat{\mathbf{F}}_{(\epsilon^{2})}^{(z_{j}z_{k})}$ and $\mathbf{e}_{\ell}$ is
an element of the orthonormal basis of $\mathbb{R}^{2}$.

Then, we look for a second order term expanded as follows

$\mathbf{q}_{(\varepsilon^{2})}\equiv\displaystyle\sum_{\begin{subarray}{c}j,k=1\\\
k\leq
j\end{subarray}}^{2}\big{(}z_{j}z_{k}\hat{\mathbf{q}}_{z_{j}z_{k}}+z_{j}\overline{z}_{k}\hat{\mathbf{q}}_{z_{j}\overline{z}_{k}}+\text{
c.c
}\big{)}+\displaystyle\sum_{\ell=1}^{2}\eta_{\ell}\mathbf{Q}_{0}^{(\eta_{\ell})}.$
(26)

Terms $\hat{\mathbf{q}}_{z^{2}_{j}}$ are azimuthal harmonics of the flow. The
terms $\hat{\mathbf{q}}_{z_{j}z_{k}}$ with $j\neq k$ are coupling terms, and
$\hat{\mathbf{q}}_{|z_{j}|^{2}}$ are harmonic base flow modification terms.
Finally, $\mathbf{Q}_{0}^{(\eta_{\ell})}$ are base flow corrections due to a
variation of the parameter $\eta_{\ell}$ from the critical point.

At this order, there exists two resonant terms, the terms proportional to
$\overline{z}_{1}z_{2}$ and $z_{1}^{2}$, which are associated with the
singular Jacobian $\mathbf{J}_{(0,m_{k})}$ for $k=1,2$. To ensure the
solvability of these terms, we must enforce compatibility conditions, i.e. the
Fredholm alternative. The resonant terms are then determined from the
resolution of the following set of bordered systems

$\begin{pmatrix}\mathbf{J}_{(0,m_{k})}&\hat{\mathbf{q}}_{k}\\\
\hat{\mathbf{q}}_{k}^{\dagger}&0\end{pmatrix}\begin{pmatrix}\hat{\mathbf{q}}_{(\mathbf{z}^{(R)})}\\\
e\end{pmatrix}=\left(\begin{array}[]{c}\hat{\mathbf{F}}^{(\mathbf{z}^{(R)})}_{(\varepsilon^{2})}\\\
0\end{array}\right),\text{
}\mathbf{z}^{(R)}\in[\overline{z}_{1}z_{2},z_{1}^{2}]^{T},$ (27)

where $e=e_{3}$ for $\mathbf{z}^{(R)}=\overline{z}_{1}z_{2}$ and $e=e_{4}$ for
$\mathbf{z}^{(R)}=z_{1}^{2}$. The non-resonant terms are computed by solving
the following non-degenerated forced linear systems

$\displaystyle\mathbf{J}_{(0,m_{j}+m_{k})}\hat{\mathbf{q}}_{z_{j}z_{k}}=\hat{\mathbf{F}}_{(\epsilon^{2})}^{(z_{j}z_{k})},$
(28)

and

$\displaystyle\mathbf{J}_{(0,0)}\mathbf{Q}_{0}^{(\eta_{\ell})}=\mathbf{G}(\mathbf{Q}_{0},\mathbf{e}_{\ell}).$
(29)

####  A.1.4 Order $\varepsilon^{3}$

At third order, there exists six degenerate terms. In our case, we are not
interested in solving for terms of third-order, instead, we will determine the
linear and cubic coefficients of the third order normal form eq. 10 from a set
of compatibility conditions.

The linear terms $\lambda_{(s,1)}$ and $\lambda_{s,2}$ and cubic terms
$c_{(i,j)}$ for $i=1,2$ are determined as follows

$\lambda_{(s,1)}=\frac{\langle\hat{\mathbf{q}}_{1}^{\dagger},\hat{\mathbf{F}}^{(z_{1})}_{(\varepsilon^{3})}\rangle}{\langle\hat{\mathbf{q}}_{z}^{\dagger},\mathbf{B}\hat{\mathbf{q}}_{z}\rangle},\text{
}\lambda_{(s,2)}=\frac{\langle\hat{\mathbf{q}}_{2}^{\dagger},\hat{\mathbf{F}}^{(z_{2})}_{(\varepsilon^{3})}\rangle}{\langle\hat{\mathbf{q}}_{2}^{\dagger},\mathbf{B}\hat{\mathbf{q}}_{2}\rangle},\text{
}c_{(i,j)}=\frac{\langle\hat{\mathbf{q}}_{2}^{\dagger},\hat{\mathbf{F}}^{(z_{i}|z_{j}|^{2})}_{(\varepsilon^{3})}\rangle}{\langle\hat{\mathbf{q}}_{i}^{\dagger},\mathbf{B}\hat{\mathbf{q}}_{i}\rangle}.$
(30)

The forcing terms for the linear coefficient are

$\hat{\mathbf{F}}^{(z_{j})}_{(\varepsilon^{3})}\equiv\sum_{\ell=1}^{2}\eta_{\ell}\Big{(}\big{[}\mathbf{N}(\hat{\mathbf{q}}_{j},\mathbf{Q}_{0}^{(\eta_{\ell}))}+\mathbf{N}(\mathbf{Q}_{0}^{(\eta_{\ell})},\hat{\mathbf{q}}_{j})\big{]}+\mathbf{G}(\hat{\mathbf{q}}_{j},\mathbf{e}_{\ell})\Big{)}.$
(31)

which allows the decomposition of
$\lambda_{(s,\ell)}=\lambda_{(s,\ell),\text{Re}}(\text{Re}^{-1}_{c}\text{Re}^{-1})+\lambda_{(s,\ell),{\delta_{u}}}(\delta_{u,c}-\delta_{u})$
for $\ell=1,2$.

The forcing terms for the cubic coefficients are

$\begin{array}[]{rl}\hat{\mathbf{F}}^{({z}_{j}|z_{k}|^{2})}_{(\varepsilon^{3})}&\equiv\big{[}\mathbf{N}(\hat{\mathbf{q}}_{j},\hat{\mathbf{q}}_{|z_{k}|^{2}})+\mathbf{N}(\hat{\mathbf{q}}_{|z_{k}|^{2}},\hat{\mathbf{q}}_{j})\big{]}\\\
&+\big{[}\mathbf{N}(\hat{\mathbf{q}}_{-k},\hat{\mathbf{q}}_{z_{j}z_{k}})+\mathbf{N}(\hat{\mathbf{q}}_{j,k},\hat{\mathbf{q}}_{-k})\big{]}\\\
&+\big{[}\mathbf{N}(\hat{\mathbf{q}}_{k},\hat{\mathbf{q}}_{z_{j}\overline{z}_{k}})+\mathbf{N}(\hat{\mathbf{q}}_{z_{j}\overline{z}_{k}},\hat{\mathbf{q}}_{k})\big{]}.\end{array}$
(32)

if $j\neq k$ and

$\begin{array}[]{rl}\hat{\mathbf{F}}^{({z}_{j}|z_{j}|^{2})}_{(\epsilon^{3})}&\equiv\big{[}\mathbf{N}(\hat{\mathbf{q}}_{j},\hat{\mathbf{q}}_{|z|_{j}^{2}})+\mathbf{N}(\hat{\mathbf{q}}_{|z|_{j}^{2}},\hat{\mathbf{q}}_{j})\big{]}\\\
&+\big{[}\mathbf{N}(\hat{\mathbf{q}}_{-j},\hat{\mathbf{q}}_{z_{j}^{2}})+\mathbf{N}(\hat{\mathbf{q}}_{z^{2}_{j}},\hat{\mathbf{q}}_{-j})\big{]},\end{array}$
(33)

for the diagonal forcing terms.

##  Appendix B Validation of the code - Comparison with the literature

The calculations made in StabFem in the sections at the main manuscript are
validated comparing the leading global mode in the geometry proposed by [5].
Moreover, the critical Reynolds number and the frequency associated are also
analysed. In the cited work, the authors use an analogous geometry with the
following parameters:

  * •

Radious of the inner jet $R_{inner}=0.5$

  * •

Diameter of the outer jet $D=0.4$

  * •

Distance between jets $L=0.1$

  * •

Ratio between velocities $\delta_{u}=1$

The linear stability analysis has been carried out imposing $m=0$, as done by
[5], so the leading global mode will be axisymmetric. The critical Reynolds
number $Re_{c}$ and the frequency $\omega$ of the leading global mode are
compared in Tab. 4. As seen, few differences can be found on the critical
Reynolds number and the frequency. The relative error in the $Re_{c}$
calculation is $1.06\%$ and the one of the frequency is $0.17\%$.

| Canton et al. (2017) [5] | Present work
---|---|---
$Re_{c}$ | $1420$ | $1405$
$\omega$ | $5.73$ | $5.72$
Table 4: Comparison of $Re_{c}$ and $\omega$ between previous work and the
present one.

The global mode is now calculated using StabFem and compared with the one
calculated by [5]. This mode can be found in figures 9, 10 and 11 on the cited
paper. As it can be seen, there are not substantial differences between the
direct modes, being both of them a vortex street with their biggest amplitude
situated 10 units downstream the exit of the jets. The adjoint mode is
concentrated within the nozzle, with its biggest amplitude situated on the
sharp corners. There is not any difference between the adjoint mode calculated
with StabFem and the one in [5]. Finally, the structural sensitivity is
similar to the one computed by [5]. It is composed by two lobes in the space
between the exit of the two jets.

Figure 12: Direct mode, adjoint mode and sensitivity of the leading global
mode studied by [5] calculated using StabFem.

## References

  * [1] Dieter Armbruster, John Guckenheimer, and Philip Holmes.  Heteroclinic cycles and modulated travelling waves in systems with o (2) symmetry.  Physica D: Nonlinear Phenomena, 29(3):257–282, 1988.
  * [2] F. Auguste, D. Fabre, and J. Magnaudet.  Bifurcations in the wake of a thick circular disk.  Theoretical and Computational Fluid Dynamics, 24:305–313, 2010\.
  * [3] A. Bogulawski and K. Wawrzak.  Absolute instability of an annular jet: local stability analysis.  Meccanica, 55:2179–2198, 2020.
  * [4] G. Buresti, G. A. Talamelli, and P. Petagna.  Experimental characterization of the velocity field of a coaxial jet configuration.  Exp. Therm. Fluid Sci., 9:135, 1994.
  * [5] J. Canton, F. Auteri, and M. Carini.  Linear global stability of two incompressible coaxial jets.  J. Fluid Mech., 824:886–911, 2017.
  * [6] C. B. da Silva, G. Balarac, and O. Métais.  Transition in high velocity ratio coaxial jets analysed from direct numerical simulations.  J. Turbul., 24:1, 2003.
  * [7] W. J. A. Dahm, C. E. Frieler, and G. Tryggvason.  Vortex structure and dynamics in the near field of a coaxial jet.  J. Fluid Mech., 241:371, 1992.
  * [8] Gerhard Dangelmayr.  Steady-state mode interactions in the presence of 0 (2)-symmetry.  Dynamics and Stability of Systems, 1(2):159–185, 1986.
  * [9] D. Fabre, V. Citro, D. Ferreira Sabino, P. Bonnefis, J. Sierra, F. Giannetti, and M. Pigou.  A Practical Review on Linear and Nonlinear Global Approaches to Flow Instabilities.  Applied Mechanics Reviews, 70(6), 02 2019.  060802\.
  * [10] CA Jones and MRE Proctor.  Strong spatial resonance and travelling waves in bénard convection.  Physics Letters A, 121(5):224–228, 1987.
  * [11] Todd Kapitula and Keith Promislow.  Spectral and dynamical stability of nonlinear waves, volume 457\.  Springer, 2013.
  * [12] N.W.M. Ko and A.S.H. Kwan.  The initial region of subsonic coaxial jets.  J. Fluid Mech., 73:305, 1976.
  * [13] Paolo Maria Mariano and Furio Lorenzo Stazi.  Computational aspects of the mechanics of complex materials.  Archives of Computational Methods in Engineering, 12(4):391–478, 2005.
  * [14] Philippe Meliga, Jean-Marc Chomaz, and Denis Sipp.  Global mode interaction and pattern selection in the wake of a disk: a weakly nonlinear expansion.  Journal of Fluid Mechanics, 633:159–189, 2009.
  * [15] Philippe Meliga, François Gallaire, and Jean-Marc Chomaz.  A weakly nonlinear mechanism for mode selection in swirling jets.  Journal of Fluid Mechanics, 699:216–262, 2012.
  * [16] Isabel Mercader, Joana Prat, and Edgar Knobloch.  Robust heteroclinic cycles in two-dimensional rayleigh–bénard convection without boussinesq symmetry.  International Journal of Bifurcation and Chaos, 12(11):2501–2522, 2002.
  * [17] A. Michalke.  Absolute inviscid instability of a ring jet with back-flow and swirl.  Eur. J. Mech. B/Fluids, 18:3–12, 1999.
  * [18] C Nore, F Moisy, and L Quartier.  Experimental observation of near-heteroclinic cycles in the von kármán swirling flow.  Physics of Fluids, 17(6):064103, 2005.
  * [19] Caroline Nore, Laurette S Tuckerman, Olivier Daube, and Shihe Xin.  The 1 [ratio] 2 mode interaction in exactly counter-rotating von kármán swirling flow.  Journal of Fluid Mechanics, 477:51–88, 2003.
  * [20] W. Olsen and A. Karchmer.  Lip noise generated by flow separation from nozzle surfaces.  AIAA J., 76:3, 1976.
  * [21] R. Örlü, A. Segalini, P. H. Alfredsson, and A. Talamelli.  On the passive control of the near-field of coaxial jets by means of vortex shedding.  Proceedings of the International Conference on Jets, Wakes and Separated Flows, ICJWSF-2008, Technical University of Berlin, Berlin, Germany, Sept. 16–19, 2008.
  * [22] Antonio Palacios, Gemunu H Gunaratne, Michael Gorman, and Kay A Robbins.  Cellular pattern formation in circular domains.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 7(3):463–475, 1997.
  * [23] J Porter and E Knobloch.  New type of complex dynamics in the 1: 2 spatial resonance.  Physica D: Nonlinear Phenomena, 159(3-4):125–154, 2001.
  * [24] J Porter and E Knobloch.  Dynamics in the 1: 2 spatial resonance with broken reflection symmetry.  Physica D: Nonlinear Phenomena, 201(3-4):318–344, 2005.
  * [25] H. Rehab, E. Villermaux, and E. J. Hopfinger.  Flow regimes of largevelocity-ratio coaxial jets.  J. Fluid Mech., 345:357, 1997.
  * [26] A. Segalini and A. Talamelli.  Experimental analysis of dominant instabilities in coaxial jets.  Phys. Fluids, 23:024103, 2011.
  * [27] J. Sierra, D. Fabre, V. Citro, and F. Giannetti.  Bifurcation scenario in the two-dimensional laminar flow past a rotating cylinder.  Journal of Fluid Mechanics, 905:A2, 2020.
  * [28] Javier Sierra, David Fabre, and Vincenzo Citro.  Efficient stability analysis of fluid flows using complex mapping techniques.  Computer Physics Communications, 251:107100, 2020.
  * [29] J. Sierra-Ausin, Lorite-Diez M., V. Citro, Jimenez J.I., and D. Fabre.  Rotating sphere.  Journal of Fluid Mechanics, -:–, 2022.
  * [30] J. Sierra-Ausín, M. Lorite-Díez, J.I. Jiménez-González, V. Citro, and D. Fabre.  Unveiling the competitive role of global modes in the pattern formation of rotating sphere flows.  Journal of Fluid Mechanics, 942:A54, 2022.
  * [31] A. Talamelli and I. Gavarini.  Linear instability characteristics of incompressible coaxial jets.  Flow Turbulence Combust., 76:221–240, 2006.
  * [32] O. Tammisola.  Oscillatory sensitivity patterns for global modes in wakes.  J. Fluid Mech., 701:251–277, 2012.
  * [33] D. Wallace and L.G. Redekopp.  Linear instability characteristics of wake-shear layers.  Phys. Fluids, 4:189––191, 1992.
  * [34] T.J. Williams, M.R.M.H. Ali, and J.S. Anderson.  Noise and flow characteristics of coaxial jets.  J. Mech. Eng. Sci., 1:2, 1969.