Léopold Shaabani-Ardali
email: leopold.shaabani-ardali@ladhyx.polytechnique.fr

Denis Sipp

Lutz Lesshafft

Vortex pairing in jets as a global Floquet instability: modal and transient dynamics

The spontaneous pairing of rolled-up vortices in a laminar jet is investigated as a global secondary instability of a time-periodic spatially developing vortex street. The growth of subharmonic perturbations, associated with vortex pairing, is analysed both in terms of modal Floquet instability and in terms of transient growth dynamics. The article has the double objective to outline a toolset for global analysis of time-periodic flows, and to leverage such an analysis for a fresh view on the vortex pairing phenomenon.

Axisymmetric direct numerical simulations (DNS) of jets with single-frequency inflow forcing are performed, in order to identify combinations of the Reynolds and Strouhal numbers for which vortex pairing is naturally observed. The same DNS calculations are then repeated with an added time-delay control term, which artificially suppresses pairing, so as to obtain time-periodic unpaired base flows for linear stability analysis. It is demonstrated that the natural occurrence of vortex pairing in nonlinear DNS coincides with a linear subharmonic Floquet instability of the underlying unpaired vortex street. However, DNS results suggest that the onset of pairing involves much stronger temporal growth of subharmonic perturbations than what is predicted by modal Floquet analysis, as well as a spatial distribution of these fast-growing perturbation structures that is inconsistent with the unstable Floquet mode. Singular value decomposition of the phaseshift operator (the operator that maps a given perturbation field to its state one flow period later) is performed for an analysis of optimal transient growth in the vortex street. Non-modal mechanisms near the jet inlet are thus found to provide a fast route towards the limit-cycle regime of established vortex pairing, in good agreement with DNS observations.

It is concluded that modal Floquet analysis accurately predicts the parameter regime where sustained vortex pairing occurs, but that the bifurcation scenario under typical conditions is dominated by transient growth phenomena.

Introduction

When the shear layer of a jet is subjected to low-level forcing at the nozzle, perturbations within a band of unstable frequencies are amplified as they travel downstream [START_REF] Michalke | Instabilität eines kompressiblen runden Freistrahls unter Berücksichtigung des Einflusses der Strahlgrenzschichtdicke[END_REF]. When the perturbation amplitude reaches nonlinear levels, the shear layer rolls up into vortices. In the case of laminar jets, forced axisymmetrically at a single frequency, a regular street of ring vortices is formed, where the passage frequency of vortices is controlled by the nozzle forcing. It has long been observed that such vortices may, in certain parameter regimes, spontaneously undergo regular pairing: two neighbouring vortices then merge into one, such that a new vortex street is formed, with a periodicity that corresponds to the subharmonic of the applied forcing frequency. This Figure 1: Experimental visualisation of pairing in a jet, snapshot of a movie by [START_REF] Schram | Aeroacoustics of subsonic jets: prediction of the sound produced by vortex pairing based on particle image velocimetry[END_REF].

scenario is illustrated in figure 1, which shows a flow visualisation from the experiments by [START_REF] Schram | Aeroacoustics of subsonic jets: prediction of the sound produced by vortex pairing based on particle image velocimetry[END_REF]. The first systematic studies of the pairing scenario have been conducted by Zaman & Hussain (1980) and Hussain & Zaman (1980). While the merging of two vortices is clearly a nonlinear process, it remains to be clarified to what extent the occurrence of regular vortex pairing is governed by linear mechanisms, and how such linear dynamics may be properly formalised.

Numerous past studies have sought to explain the pairing phenomenon by investigating the consistency of the paired flow state itself. [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF] developed a theoretical framework based on weakly nonlinear interaction between the forced fundamental instability wave and its subharmonic. His formalism leads to a phase-relation criterion for subharmonic resonance, which has subsequently been validated experimentally by [START_REF] Husain | Subharmonic resonance in a shear layer[END_REF], [START_REF] Raman | Axisymmetric jet forced by fundamental and subharmonic tones[END_REF], [START_REF] Paschereit | Experimental investigation of subharmonic resonance in an axisymmetric jet[END_REF] and others.

A different approach consists in considering the bifurcation process from a harmonic unpaired to a subharmonic paired state. This viewpoint leads to the question whether an array of vortices convecting at the imposed forcing frequency is unstable with respect to subharmonic perturbations. In forced jets, vortex pairing occurs as a self-sustained process, and the pairing location is stationary in a spatially non-periodic flow. Therefore, the appropriate framework will have to be based on either locally absolute/convective, or fully global analysis. [START_REF] Brancher | Absolute and convective secondary instabilities in spatially periodic shear flows[END_REF] investigate the absolute nature of pairing instability in a periodic array of Stuart vortices, as a model for a rolled-up two-dimensional plane shear layer. In order to apply the notion of absolute/convective instability [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF], commonly used for steady configurations, they conduct their analysis in a co-moving frame of reference, where the vortices are stationary. By varying the concentration factor of the vortices, they show that the required backflow rate for absolute instability decreases as the vortex concentration increases. This formalism pertains to a spatially periodic array of vortices and does not account for viscous effects.

In the present study, we aim to describe the onset of vortex pairing as a secondary instability in a global analysis framework. As the underlying basic state, the unpaired vortex street resulting from the primary shear instability, is time-periodic, a classical Floquet formalism is employed. Modal instability as well as transient growth scenarios will be explored within this framework. Global instability analysis has become a standard approach in the context of steady base flows [START_REF] Theofilis | Global linear instability[END_REF]. Linearisation of the governing flow equations around a steady base flow yields an autonomous operator; the spectrum of this operator indicates the possibility of perturbation growth in the long-time limit. However, the base flow in our case is not steady but periodic in time, the periodicity being controlled by the harmonic forcing of the jet. Therefore, the global analysis method must be adapted to account for this periodicity, by use of the Floquet theorem [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF]. A global Floquet formalism has been successfully applied by [START_REF] Barkley | Three-dimensional Floquet stability analysis of the wake of a circular cylinder[END_REF] to the secondary instability of a cylinder wake around Re = 200 and more recently by [START_REF] Jallas | Linear and nonlinear perturbation analysis of the symmetry breaking in time-periodic propulsive wakes[END_REF] to a study of secondary instabilities in wake flows.

Another possible scenario for the onset of instability would be through transient growth; this phenomenon has first been described in parallel flows by [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF]: if the linearised flow operator is non-normal, then, even though the system is stable, in the sense that at an infinite time horizon any perturbation decays towards zero, intial perturbations can be greatly amplified over a short time. For stationary jets, Nichols & Lele (2011) and [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF] studied this phenomenon and shown its significant role in the bifurcation. However, we aim at studying transient growth of a timeperiodic jet. Several authors [START_REF] Barkley | Direct optimal growth analysis for timesteppers[END_REF][START_REF] Blackburn | Convective instability and transient growth in steady and pulsatile stenotic flows[END_REF][START_REF] Arratia | Transient perturbation growth in time-dependent mixing layers[END_REF][START_REF] Johnson | Crow instability: nonlinear response to the linear optimal perturbation[END_REF]) applied direct-adjoint looping in order to identify finite-time optimal perturbations in unsteady base flows. In the present paper, we use a technique that allows us to retrieve optimal perturbations based on direct time-stepping alone, and which takes full advantage of the time-periodicity of the underlying base flows.

As a prerequisite for our analysis, periodic base flow states without vortex pairing must be computed even in situations where such pairing arises naturally. Standard methods for the computation of periodic flow states may be based on Newton-Picard shooting methods, as described in [START_REF] Roose | A Newton-Picard shooting method for computing periodic solutions of large-scale dynamical systems[END_REF]; [START_REF] Lust | An adaptive Newton-Picard algorithm with subspace iteration for computing periodic solutions[END_REF] and applied to several flows in [START_REF] Sánchez | Newton-Krylov continuation of periodic orbits for Navier-Stokes flows[END_REF]; [START_REF] Sánchez | On the multiple shooting continuation of periodic orbits by Newton-Krylov methods[END_REF]. One alternative is the harmonic balance technique, in which several temporal Fourier components of a given periodicity are calculated simultaneously [START_REF] Hall | Computation of unsteady nonlinear flows in cascades using a harmonic balance technique[END_REF][START_REF] Thomas | Nonlinear inviscid aerodynamic effects on transonic divergence, flutter, and limit-cycle oscillations[END_REF]. In a recent publication [START_REF] Shaabani-Ardali | Time-delayed feedback technique for suppressing instabilities in time-periodic flow[END_REF], we describe how subharmonic fluctuations in time-stepping simulations can be efficiently suppressed by way of timedelayed feedback control ; this technique is employed here in order to construct the periodic base states.

The paper is organised as follows: a comprehensive literature review is provided in § 2, in order to delineate the context of our investigation. Section 3 presents a systematic study of the parameter regimes where vortex pairing is observed in direct numerical simulations (DNS) of laminar jets at moderate Reynolds numbers. Modal Floquet analysis is performed in § 4, and it is shown that the Floquet-unstable parameter regime coincides with the observation of pairing in the DNS. Non-modal transient growth dynamics are investigated in § 5, and their relevance for the bifurcation scenario is demonstrated.

Review of the literature on vortex pairing

Discovery

The pairing of vortex rings in jets was described for the first time by [START_REF] Becker | Vortex evolution in a round jet[END_REF]: in jet experiments at moderate Reynolds number, acoustic single-frequency forcing was observed to give rise to regular vortex formation and subsequent pairing. [START_REF] Winant | Vortex pairing: the mechanism of turbulent mixinglayer growth at moderate Reynolds number[END_REF] investigated the same phenomenon in a plane mixing layer: pairing was found to occur intermittently in these experiments, punctuated by occasional "shredding" events, i.e. the destruction of vortex cores by a subharmonic strain field, which interrupt the pairing process. The authors proposed a phenomenological model based on Stuart vortices. In the context of jets, [START_REF] Petersen | Influence of wave dispersion on vortex pairing in a jet[END_REF] examined the influence of higher Reynolds numbers, and they inferred an argument on the basis of wave dispersion intended to predict the location of vortex pairing.

More detailed experimental investigations of vortex pairing in jets were carried out by Zaman and Hussain (Zaman & Hussain 1980;Hussain & Zaman 1980). These authors considered high Reynolds number (Re = O(10 4 )) jets with a thin initial mixing layer, forced at a single frequency. Pairing was found to arise in two distinct frequency bands, one around St θ = f θ/U ≈ 0.012 and the other around St D = f D/U ≈ 0.85. These Strouhal numbers are formed with the jet exit velocity U , forcing frequency f , and either the initial shear layer thickness θ or the jet diameter D. According to their characteristic scaling, the two bands were identified with a "mixing layer mode" and a "jet column mode", respectively. In the former case, the vortices are very thin and dissipate quickly, whereas in the latter case, their radial extent is comparable to the jet radius, and their viscous dissipation takes place over a much longer travel distance. Both articles describe in much detail the vortex dynamics, their trajectories and velocities, as well as the transition to turbulence. When turbulence sets in close to the nozzle, the pairing becomes intermittent. High-quality flow visualisations of vortex pairing in jets, at Reynolds number 2300, are shown by [START_REF] Meynart | Speckle velocimetry study of vortex pairing in a low-Re unexcited jet[END_REF].

Vortex pairing is also a common event in plane shear layers. [START_REF] Ho | Subharmonics and vortex merging in mixing layers[END_REF] found that very low-amplitude subharmonic forcing in their shear layer experiments led to a vigorous flow response in the form of regular pairing, associated with a strongly increased spreading of the mean flow. Intermittently, simultaneous coalescence of several vortices occurred, and was described as "collective interaction". Similar observations had been made in jet experiments by [START_REF] Kibens | Discrete noise spectrum generated by acoustically excited jet[END_REF]. [START_REF] Ho | Subharmonics and vortex merging in mixing layers[END_REF] established experimentally that spatial growth of the subharmonic component only occurs in situations where its phase velocity is equal to that of the fundamental flow perturbation.

Interpretation in terms of wave interaction

Prior to theoretical explanations, [START_REF] References Arbey | Active cancellation of pure tones in an excited jet[END_REF] demonstrated experimentally the importance of the phase difference between fundamental and subharmonic perturbations for the onset of vortex pairing. A jet at Reynolds number 17500 was forced at moderate amplitude (about 2% of the centreline velocity) at two frequencies ω and ω/2. The spatial growth of perturbations at both frequencies was found to be strongly dependent on the relative phase of the applied forcing. Thanks to a numerical model, [START_REF] Mankbadi | On the interaction between fundamental and subharmonic instability waves in a turbulent round jet[END_REF] argued that pairing arises when the subharmonic component acquires sufficient energy, both from the fundamental wave and from the mean flow, to become the largest-amplitude perturbation in the jet. He observed that one or more stages of subsequent pairing can occur; in his framework, the number of stages and their spatial localisation depends on the Strouhal number.

A deeper theoretical understanding of the wave interaction involved in vortex pairing was reached by [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF], who formulated a weakly nonlinear model for the spatial development of fundamental and subharmonic instability waves in a parallel mixing layer. This model reflects the role of the phase shift between the two waves in triggering resonance, resulting in either pairing or "shredding" of vortices. It was predicted that the fundamental wave needs to reach a critical amplitude before subharmonic perturbations may phase-lock and grow. The same conclusions are supported by the more general theory of [START_REF] Cheng | Subharmonic instabilities of finite-amplitude monochromatic waves[END_REF].

The model of [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF] was confirmed experimentally by numerous studies: in a mixing layer, [START_REF] Husain | Subharmonic resonance in a shear layer[END_REF] showed that simultaneous forcing of a fundamental frequency and its subharmonic could either enhance or attenuate the pairing and shredding phenomena. Subsequently, detailed statistical analysis of experimental data was performed, first for uncontrolled "natural" perturbations in a mixing layer [START_REF] Hajj | Subharmonic growth by parametric resonance[END_REF], then for explicitly forced fundamental and subharmonic waves in the same setup [START_REF] Hajj | Fundamental-subharmonic interaction: of phase relation[END_REF]. These two studies gave clear evidence of a parametric resonance, determined by the phase difference between both waves. [START_REF] Husain | Experiments on subharmonic resonance in a shear layer[END_REF] carried out similar investigations of jets with a very thin shear layer, confirming that pairing was amplified for a large band of phase differences, but attenuated for a narrow band of phase differences. Moreover, these authors studied the influence of a slight frequency detuning in the subharmonic forcing, finding that the occurrence of vortex pairing depended on the instantaneous phase difference in the forcing. A parametric experimental study of pairing in jets was conducted by [START_REF] Raman | Axisymmetric jet forced by fundamental and subharmonic tones[END_REF], who varied the Strouhal number, phase difference and forcing amplitude for both the fundamental and the subharmonic instability wave. Consistent with the theoretical predictions, it was shown that a critical fundamental amplitude was necessary to trigger subharmonic resonance, the growth rate being controlled by the phase difference at small forcing amplitude. When the forcing was strong enough, however, the subharmonic growth became independent of the phase difference. [START_REF] Paschereit | Experimental investigation of subharmonic resonance in an axisymmetric jet[END_REF] confirmed these results, and further demonstrated that the subharmonic growth draws its energy from the mean flow, whereas the fundamental wave merely acts as a catalyst. In all these studies, explicit forcing of the subharmonic wave controlled the location of vortex pairing in turbulent jets.

The modal interaction framework was used by [START_REF] Bradley | Phase-locking in a jet forced with two frequencies[END_REF] to study interactions between more than two frequencies, or between frequencies different from the fundamental ω and its subharmonic ω/2. These authors experimentally studied a jet forced either at ω and ω/2 or at ω and ω/3 and studied the influence of frequency, amplitude ratio and phase shift. In the ω and ω/3 forcing case, they found more diverse behaviour, with collective interactions, or pairing between vortices of different sizes.

Further developments on pairing

Chaotic behaviour

In less controlled configurations, vortex pairing events are often observed to be irregular and intermittent. [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF], 1996) conducted jet experiments with single-frequency excitation; depending on the Strouhal number St D and on the forcing amplitude a f , different types of behaviour were reported, as summarised in figure 2. Regular dynamics were found to arise over large parameter regions, namely no pairing (FO regime in figure 2), stable pairing (SP) and stable double pairing (SDP), the latter referring to the occurrence of two successive stages of vortex pairing. Irregular dynamics were observed either in the mild form of "aperiodic modulations" of the first (AM) or the second pairing stage (SPMQ), or in more erratic ways, categorised as intermittency, chaos (QCA) and "nearly-periodic modulations" (NPMP).

The results of [START_REF] Broze | Transitions to chaos in a forced jet: intermittency, tangent bifurcations and hysteresis[END_REF] do not appear to depend significantly on the Reynolds number within their operating conditions of 11000 Re 90000. Drawing on chaos theory, the authors characterised the pairing dynamics in terms of attractors. For irregular scenarios, it was demonstrated that the occurrence of pairing of two vortices is strongly influenced by previous pairing events, implying that upstream-directed feedback is involved in the subharmonic growth. This observation is fully consistent with the waveinteraction model discussed in § 2.2. [START_REF] Narayanan | Chaos control in open flows-experiments in a circular jet[END_REF] attempted to stabilise the pairing dynamics in chaotic regimes.

Pairing-related jet noise

High-speed jets are potent sources of noise, and the role of vortex pairing as an aeroacoustic source mechanism has received much attention. [START_REF] Bridges | Roles of initial condition and vortex pairing in jet noise[END_REF] Figure 2: Pairing dynamics observed experimentally in turbulent jets at Re = O(10 4 ), forced at a Strouhal number St D with amplitude a f . No subharmonic forcing is applied. From [START_REF] Broze | Transitions to chaos in a forced jet: intermittency, tangent bifurcations and hysteresis[END_REF]. determined that the radiated noise is dominated by pairing events only in cases where the initial shear layer is laminar; in fully turbulent jets, vortex pairing appears to be acoustically unimportant. The dominant role of vortex pairing in laminar jet noise was further analysed numerically by [START_REF] Bogey | Influence of nozzle-exit boundary-layer conditions on the flow and acoustic fields of initially laminar jets[END_REF]. A detailed description of sound generation from vortex interaction in jets was given by [START_REF] Inoue | Sound generation by the leapfrogging between two coaxial vortex rings[END_REF], and model-based predictions of radiated sound levels were derived by [START_REF] Schram | Theoretical/empirical prediction and measurement of the sound produced by vortex pairing in a low Mach number jet[END_REF].

Kinematical modelling of vortex-ring interaction

Vortex interaction can be conveniently modelled by use of the Biot-Savart law. An account of early applications of such approaches to vortex rings is given by [START_REF] Shariff | Vortex rings[END_REF]; these studies were largely concerned with the initial roll-up or with the collective motion of a limited number of co-axial vortices. A common representation characterises a vortex ring by its core centre position and core radius, from which a semianalytical induction law can be derived [START_REF] Saffman | Vortex dynamics[END_REF]. Contrary to vortex elements in a two-dimensional plane, axisymmetric vortex rings move at a self-induced velocity, which becomes infinitely large as the core radius approaches zero. Some simplifications arise from the assumption that the core radius is much smaller than the ring radius.

Within the limiting assumptions of inviscid flow and compact cores, such a conceptual model then yields low-dimensional systems representing the mutually induced motion of a collection of vortex rings. The "leapfrogging" interaction between two rings, when one passes through the other, corresponds to the early nonlinear stage of pairing in a jet before actual merging occurs. A model of a plane jet, consisting of counter-rotating vortex pairs, was studied by [START_REF] Tophøj | Instability of vortex pair leapfrogging[END_REF], who formulated a linear stability criterion for the occurrence of leapfrogging. Similarly, [START_REF] Borisov | The dynamics of vortex rings: Leapfrogging, choreographies and the stability problem[END_REF][START_REF] Borisov | The dynamics of vortex rings: leapfrogging in an ideal and viscous fluid[END_REF] derived stability criteria for sets of two and three co-axial vortex rings, and validated these with respect to numerical simulations of viscous flow at high Reynolds number. Their work was completed by [START_REF] Cheng | Leapfrogging of multiple coaxial viscous vortex rings[END_REF], who explicitly portrayed parameter regions in which leapfrogging could occur, depending on the Reynolds number and the aspect ratio of the vortex arrangement. [START_REF] Meunier | A merging criterion for twodimensional co-rotating vortices[END_REF] discovered a merging criterion for two co-rotating vortices in a two-dimensional plane.

Pairing of helical vortices

In several industrial applications, such as helicopter rotors or wind turbines, a wake composed of several nested helical vortices is formed [START_REF] Vermeer | Wind turbine wake aerodynamics[END_REF]. Further downstream, these helical vortices diffuse and can interact together as well as with the hub vortex [START_REF] Delbende | Merging of two helical vortices[END_REF][START_REF] Felli | Mechanisms of evolution of the propeller wake in the transition and far fields[END_REF].

Formally, these vortices form a time-periodic three-dimensional flow, which can be regarded as a steady state in a frame of reference rotating with the blades. In addition, assuming no interaction between vortices, they diffuse slowly in the far wake, allowing for a quasi-static approximation (Selçuk et al. 2017b) when neglecting variations of the vortex structure in the axial direction. Therefore, Selçuk et al. (2017a) carried out a global stability analysis of these quasi-static states reduced to steady flows. At low pitch values, they found the existence of unstable modes that trigger leapfrogging, overtaking and eventually pairing when superposed onto the quasi-static base flow.

Vortex pairing in simulations at various Reynolds and Strouhal numbers

In this section, we give a general overlook of the vortex pairing phenomenon in axisymmetric laminar jets. After explaining the numerical simulation details, we show that, depending on the jet parameters and the forcing considered, vortex pairing can occur or not. When it does, we take a close look on a pairing sequence. Eventually, we finish with a parametric study to characterise its occurrence domain.

Setup of direct numerical simulations

Simulations were carried out using Nek5000 (Nek 5000), an incompressible spectral element code. A perfectly axisymmetric jet is described in cylindrical coordinates (z, r), z being the main flow direction and r being the radial distance from the jet axis. The flow is assumed to be governed by the incompressible axisymmetric Navier-Stokes equations with zero azimuthal velocity, written in dimensionless form as

∂v ∂t + v ∂v ∂r + u ∂v ∂z = - ∂p ∂r + 1 Re 1 r ∂ ∂r r ∂v ∂r - v r 2 + ∂ 2 v ∂z 2 , (3.1a) ∂u ∂t + v ∂u ∂r + u ∂u ∂z = - ∂p ∂z + 1 Re 1 r ∂ ∂r r ∂u ∂r + ∂ 2 u ∂z 2 , (3.1b) 1 r ∂(rv) ∂r + ∂u ∂z = 0, (3.1c)
where the velocity u has axial and radial components u and v, and p denotes pressure. The jet diameter D and the inlet centerline velocity U 0 are used to render the flow problem nondimensional, leading to a definition of the Reynolds number as Re = U 0 D/ν, with ν the kinematic viscosity. The computational domain extends, unless specified otherwise, over 40 × 5 diameters in the axial and radial directions, respectively, and it is discretised with 17600 spectral elements, each containing 64 mesh points.

Boundary conditions are specified as follows. In the inlet plane, z = 0, a hyperbolic tangent velocity profile is imposed. Its amplitude is modulated in time as

u(z = 0, r, t) = 1 2 1 -tanh 1 4θ 0 r - 1 4r (1 + A cos (ω f t)) e z , (3.2)
where A is the forcing amplitude, θ 0 = 0.025 is the initial non-dimensional shear layer thickness and ω f is the forcing frequency. The forcing period is given by T = 2π/ω f , and the Strouhal number is defined as St D = ω f D/(2πU 0 ). On the jet centreline, axisymmetric boundary conditions are imposed as

∂u ∂r = v = ∂p ∂r = 0 at r = 0. (3.3)
In the outlet plane, z max = 40, and on the lateral boundary, r max = 5, the standard outflow formulation provided by Nek5000 is employed. This prescribes a stress-free condition 1 Re ∂v ∂r p = ∂u ∂r = 0 at r = r max and 1 Re ∂u ∂z p = ∂v ∂z = 0 at z = z max .

(3.4)

The flow configuration is thus characterised by the Reynolds number Re, the Strouhal number St D , the initial shear layer thickness θ 0 and the forcing amplitude A. A single value θ 0 = 0.025 is used throughout this study.

Two distinct behaviours

Depending on flow parameters and forcing Strouhal number, rolled-up vortices may spontaneously undergo subsequent pairing. In the absence of free-stream turbulence, and if the harmonic forcing is well-controlled, this pairing takes place in a perfectly regular fashion.

In cases where pairing occurs, two neighbouring vortices merge into one, such that the passage frequency of vortices downstream of the pairing location is exactly half that of the imposed forcing. An example, obtained by direct numerical simulation with parameters Re = 2000, St D = 0.6 and A = 5%, is shown in figure 3a. If the forcing is characterised by the time period T , the "paired state" is globally 2T -periodic. The velocity field of a paired state will be denoted u p .

A different behaviour is observed when the Reynolds number is lowered to Re = 1300, as shown in figure 3b. Vortices roll up close to the nozzle and advect downstream, until they are dissipated by viscosity, but no spontaneous pairing is observed at this parameter setting. Such a flow state is globally T -periodic and will be called hereafter an "unpaired state". Its velocity field will be denoted u u .

When the natural time-asymptotic flow state for a given set of parameters involves pairing, it is still possible to recover an unpaired state as an alternative solution of the flow equations, by use of time-delayed feedback control (see § 4.1 and Shaabani-Ardali et al. (2017)). The unpaired state obtained in such a way for the previous configuration with Re = 2000 is shown in figure 3c. This solution is an exact solution of the Navier-Stokes equations without the time-delayed feedback.

The dynamics of vortex pairing

The process of the pairing of two vortices is inspected from snapshots of the vorticity, presented in figure 4, at four different phases of one pairing cycle. At t = 0 (figure 4a), two successive vortices, located around z = 1 and 2, have rolled up due to the primary shear instability. One half-cycle of the forcing period later (figure 4b), the leading vortex has slightly expanded radially, while the trailing one has contracted. This movement is accompanied by a deceleration of the expanding vortex, and inversely an acceleration of the contracting vortex, through the influence of the vortex ring radius on its self-induced propulsion. The same process continues at t = T (figure 4c), when the trailing vortex begins to pass through the interior of the leading vortex. At t = 3T /2 (figure 4d), both vortices are in the process of merging into one, which is largely achieved at the end of the cycle (figure 4a).

Conceptual arguments for the occurrence of vortex pairing in the literature are typically based on the interaction between fundamental and subharmonic fluctuations, and the possibility of energy transfer to the latter [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF]. Spatial energy variations of the fundamental ω f and the subharmonic ω f /2 Fourier modes are readily extracted from the present numerical simulations. These are presented in figure 5, for the configuration Re = 2000, St D = 0.6 and A = 5%, in the form of radially integrated kinetic energy. This plot allows the distinction of various stages in the pairing process. Immediately downstream of the inlet, the fundamental mode grows from its forced initial amplitude to its peak value at z = 1.4. This streamwise position may be identified with the shedding of a fully formed vortex (see figure 4). The subharmonic component experiences strong exponential growth over the same interval, starting from a much lower level, as the boundary condition imposes zero amplitude at the inlet. fundamental mode decreases over the distance 1.4 < z < 3.1, although a local maximum is found at the peak location of the subharmonic mode. As pairing is accomplished, at z > 3.1, both the fundamental and the subharmonic mode decay slowly in z, both at a similar rate, due to viscous dissipation of the convecting vortices.

Consistent with these observations, [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF] argued that a growth of subharmonic perturbations must be fed by energy transfer from the fundamental mode, which requires that both modes propagate at the same phase velocity. Phase velocities of fundamental and subharmonic fluctuations in the present simulation are compared in figure 6 as solid and dashed lines. The reported phase velocities are measured, for each streamwise location, at the radial distance where each respective Fourier mode has its maximum amplitude along r. The values for both modes match quite closely throughout the relevant interval upstream of the pairing location. They continue to match in the downstream region, but this is only the result of the fundamental mode being slaved to its subharmonic counterpart as a passive harmonic.

For comparison, energy and phase velocity results obtained for the fundamental mode of the corresponding stabilised (unpaired) flow are also presented in figures 5 and 6 (dotted lines), alongside the values found in the paired state. Upstream of the roll-up location, the fundamental modes show identical energy growth in both configurations. Between the roll-up and the pairing locations, the fundamental energy decay in the paired case is stronger than in the unpaired case, which again indicates that the growth of subharmonic perturbations feed on the energy of the fundamental mode. However, downstream of the pairing location, the fundamental mode in the stabilised case decays significantly faster, and at an increasing rate. This difference is also visible when one compares figures 3a and 3c: it appears that the lower frequency and the stronger circulation of the paired vortices hinders their diffusion and allow them to be sustained longer. When comparing the phase velocities in figure 6, it is seen that the unpaired fundamental mode propagates faster than its paired counterpart upstream of the pairing location, whereas their velocities are again equal in the downstream region. The discrepancy upstream of the pairing can be linked to different positions of the vortices: in the paired case, the vortices are radially more expanded than in the unpaired case, therefore moving at a slower speed. Vortex pairing does not only influence the flow downstream, but also upstream.

Parametric study

All simulations presented in this section were carried out on a domain with z max = 15 in order to save computational resources. Test runs with z max = 40 were performed for selected cases, showing no effect of domain truncation on the results presented here.

Effect of Reynolds number and Strouhal number

The effect of Reynolds and Strouhal numbers on the onset of vortex pairing is investigated first, for a fixed forcing amplitude A = 5%. Direct numerical simulations were run, first in an exploratory fashion for many (Re, St D ) combinations. After an initial Close to the instability thresholds, it becomes difficult to precisely characterise the flow behaviour, because of the long simulation times needed to achieve convergence. When the Reynolds number is about 1350, computations were run on full-length domains (40D), and even after several hundred forcing periods, the paired or unpaired nature of the final flow state cannot be determined. For instance, at St D = 0.6, the final state for Re = 1350 is unpaired, for Re = 1400 it is paired. For intermediate values Re = 1360, 1372 or 1375 (close to threshold Re c = 1371 predicted in § 4.4), the pairing location moves gradually downstream as time evolves, but with no indication whether it will eventually become stationary. These ambiguous data points are not displayed in figure 7.

As mentioned in § 2.1, a first parametric study of vortex pairing in terms of Reynolds and Strouhal numbers has been carried out experimentally by Zaman & Hussain (1980), for Reynolds numbers greater than 10 4 and thin initial mixing layer θ/D ≈ 0.25%). They found that pairing could occur in two frequency bands, one characterised by a Strouhal number St θ = f θ/U based on the shear layer thickness, the other by a Strouhal number St D = f D/U based on the jet diameter. Pairing has been reported for St θ around 0.012 and for 0.75 St D 1.0. With our choice of the initial shear layer thickness being 5% of the diameter, these two regimes are only weakly separated, which explains why our findings of a single band of instability at a given Reynolds is coherent; both scalings are indicated in figure 7. Again for Reynolds numbers greater than 10 4 , and for a 5% forcing amplitude, [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF], found several flow regimes depending on St D (see figure 2): no pairing for St D < 0.52, aperiodic modulations and coexistence of different states (stable or modulated pairing) for 0.52 < St D < 0.77, stable pairing for 0.77 < St D < 0.97, stable pairing with quarter-harmonic modulations for 0.97 < St D < 1.02 and stable double pairing for St D > 1.02, with some intermittent cases. The St θ -dependance is not documented in that study. In our computations, we do not encounter such a richness of scenarios, because our inlet condition is fully laminar and time-periodic, and because our Reynolds number is one decade lower; however, the occurrence of pairing, reported in figure 7, is consistent with the experiments of [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF]. Multiple stages of successive pairing are never observed in our computations, even at St D > 1. This difference with respect to the experiments is certainly owing to the lower Reynolds number, as [START_REF] Cheng | Leapfrogging of multiple coaxial viscous vortex rings[END_REF] demonstrated that viscosity inhibits pairing. However, the modulated states described by [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF] for St D ≈ 1 may be linked to our observation of fluctuations of the pairing position (triangle in figure 7). [START_REF] Cheng | Leapfrogging of multiple coaxial viscous vortex rings[END_REF] investigated numerically the leapfrogging of coaxial vortex rings. For a system of two adjacent vortices, they systematically documented the occurrence of leapfrogging as a function of Reynolds number and vortex separation. The Reynolds number, defined as the ratio of the ring vorticity over the viscosity and chosen of the order of 10 3 is related to the product between Reynolds and Strouhal numbers in our present notation. The ratio of the vortex spacing over the ring radius is related to the inverse of the Strouhal number as defined here. Cheng et al. found that leapfrogging occurs only above a critical Reynolds number, and for a narrow Strouhal band that increases with Reynolds number. This is consistent with our findings (figure 7), as pairing can be understood as an advanced stage of leapfrogging (figure 4).

Effect of forcing amplitude

The effect of forcing amplitude is investigated by including computations with two additional values, A = 1% and 10%. Instability maps in the (Re, St D ) plane for these configurations are displayed in figure 8a and8b. An increase in the forcing amplitude is seen to shift the onset of pairing towards lower Reynolds number values. Fundamental mode (i = 1) for A = 10% Fundamental mode (i = 1) for A = 5% Fundamental mode (i = 1) for A = 1% Subharmonic mode (i = 1/2) for A = 10% Subharmonic mode (i = 1/2) for A = 5% Subharmonic mode (i = 1/2) for A = 1% Another effect of increasing A is to move the roll-up and pairing locations upstream, as can be seen in figure 9. In the same way as discussed for figure 5, the energy of the fundamental mode in all cases grows from the inflow towards its maximum at the roll-up location, whereas the subharmonic mode peaks at the location of pairing. Increasing the fundamental amplitude at the inflow reduces the distance needed before roll-up, and it catalyses the subharmonic mode, inducing an earlier pairing.

However, by comparing the difference between the modes at A = 1% and 5%, and between the modes at A = 5% and 10% in figure 9, it is anticipated that a further increase in A will only marginally change the paired flow behaviour. This is consistent with [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF]: as shown in figure 2, increasing the forcing amplitude above 5% does not induce significant topological changes in the final state, but below 1% forcing, no stable pairing is observed in the experiments. [START_REF] Raman | Axisymmetric jet forced by fundamental and subharmonic tones[END_REF] also found that a critical minimal amplitude of the fundamental forcing was required to trigger pairing. Therefore, an expansion of the study to lower amplitude levels A < 1% could be of interest; however, this would require higher values of Re to be considered, rendering the assumption of laminar flow increasingly doubtful.

In their leapfrogging study of two vortex rings, [START_REF] Cheng | Leapfrogging of multiple coaxial viscous vortex rings[END_REF] varied the vortex thickness, and thereby the vortex concentration, which is similar to varying the amplitude of forcing. Consistent with the present study, they concluded that stronger vortices undergo pairing at lower Reynolds numbers.

Effect of inlet noise

In configurations where vortex pairing does not arise intrinsically, the flow may still be receptive to low-level subharmonic extrinsic perturbations, in the sense of a "slightly damped oscillator" [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF], and exhibit vortex pairing in their presence. In the following, this receptivity is probed by imposing a random noise in addition to the fundamental forcing at the inlet, such that the inflow condition is prescribed as Four distinct types of the flow response are observed: (i)At low Reynolds number, the noise barely impacts the flow behaviour. For instance, at Re = 500 and for both noise levels, no significant departure from the purely harmonically forced jet is observed. (ii)The noise induces a subharmonic modulation of the vortices, but the flow diffuses too quickly for pairing to occur, for instance in the case Re = 750 and A = 0.1%. (iii)Intermittent pairing is triggered, for instance at Re = 1000 and A = 0.1%, or at Re = 750 and A = 1.0% (figure 10). In the latter case, pairing occurs far downstream, where the vortices are indistinct due to diffusion. (iv)Continuous pairing is sustained by noise input, but its location fluctuates in time. This is observed for Re = 1000 with A = 1%, and for Re = 1300 with A = 0.1% (figure 11) as well as 1.0%. As described experimentally by [START_REF] Ho | Subharmonics and vortex merging in mixing layers[END_REF]; [START_REF] Husain | Subharmonic resonance in a shear layer[END_REF], other phenomena such as shredding, where one single vortex "escapes" between two successful pairing events, or collective interaction, where more than two vortices interact at once, can be observed, for instance in figure 12.

u(r, t) = 1 2 1 -tanh 1 4θ 0 r - 1 4r (1 + A cos (ω f t) + noise (t)) e z , (3.5)

Vortex pairing as an unstable global Floquet mode

In this section, it is investigated whether the onset of vortex pairing can be described as the manifestation of a global Floquet instability of the periodic unpaired state. This unpaired state must first be computed for a given combination of flow parameters ( § 4.1). After a short reminder of Floquet theory ( § 4.2), and a presentation of the numerical implementation ( § 4.3), the linear Floquet stability of the unpaired vortex street is analysed ( § 4.4), and the instability characteristics are compared to the observations documented in the previous section. A fully synchronised paired state can be decomposed into one component that is Tperiodic and another that is only 2T -periodic,

u p (x, t) = n u T n (x) exp (inω f t) + n u 2T n (x) exp i 2n + 1 2 ω f t , (4.1)
with n = 0, ±1, ±2, . . . , ±N.

In a T -periodic unpaired state, the second sum is zero. Time-delayed feedback control is applied by adding a forcing term of the form This solution is therefore a full solution of the Navier-Stokes equations. An example of such a stabilised unpaired state is shown in figure 3c.

f = -λ(u(t) -u(t -T )), ( 4
The value of λ in equation (4.2) affects the convergence of the stabilisation procedure. In a previous publication [START_REF] Shaabani-Ardali | Time-delayed feedback technique for suppressing instabilities in time-periodic flow[END_REF], it has been demonstrated that very small values of λ do not lead to convergence, as the resulting damping is insufficient to counteract the natural instability of the system. Similarly, too large values of λ result in overshooting of the damping force, which also inhibits convergence. An optimal value λ = 0.0432ω f was derived from a model problem, and has been used in the present calculations.

Floquet framework

The Floquet stability problem for a T -periodic unpaired base flow (U u (t), P u (t)) is set up by superposing small-amplitude perturbations (u , p ), which are governed to leading order by the linear equations

∂u ∂t + (U u (t) • ∇) u + (u • ∇) U u (t) = -∇p + 1 Re ∆u , ∇ • u = 0. (4.3)
The following boundary conditions are implemented. In the inlet plane, the flow is unperturbed, u (r, z = 0, t) = 0. We do not allow perturbations directly at the inlet, since we consider only the behaviour of intrinsic perturbations. On the centreline of the jet, r = 0, axisymmetric boundary conditions as in the nonlinear case are imposed, ∂ r u = v = ∂ r p = 0. In the outlet plane z = 40 and on the lateral boundary r = 5, stress-free conditions (3.4) are chosen.

The equations are written in compact form as ∂q ∂t = L(t)q , (4.4)

where q = (u , p ) represents the perturbation state vector. According to Floquet theory [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF], one may seek modal solutions of (4.4) in the form q (t) = P (t)e At q (0), (4.5) with P (t) a T -periodic and A a constant matrix. Noting that q (T ) = e AT q (0), (4.6) due to P (T ) = P (0) = Id, the time-shift operator Φ = e AT is introduced, such that q (nT ) = Φ n q (0). The eigenvalues µ i of Φ are known as Floquet multipliers, and the associated eigenmodes v i are the Floquet modes of the system (4.4). For a modal perturbation q (0) = (ṽ i , pi ), the Floquet multiplier µ i such that q ((n+1)T ) = µ i q (nT ) represents the complex amplitude gain over one cycle period. Therefore, the stability of the system is indicated by the modulus of µ i : if all Floquet multipliers have a modulus lower than unity, all perturbations decay at long time and the system is stable. Floquet modes with an associated |µ i | > 1 experience exponential temporal growth.

The complex phase of a Floquet multiplier, arg(µ i ), characterises the time-periodicity of its associated mode. At zero phase, the mode evolves with the same periodicity as the base flow, and may be qualified as being harmonic. A phase of arg(µ i ) = π, indicates that two base flow periods are needed to complete one perturbation cycle, and the mode therefore evolves as ω f /2. Such a mode is qualified as being subharmonic. As vortex pairing is a 2T -periodic phenomenon, subharmonic Floquet modes are expected to arise. Floquet multiplier phases that are not integer multiples of π characterise modes with periodicities unrelated to that of the base flow.

A standard result of Floquet theory applied to the linearised Navier-Stokes equations is that the time derivative of the base flow, (∂ t U u , ∂ t P u ), represents a neutral Floquet mode of the system, with µ = 1. However, such a mode does not exist in the present problem, because it is inconsistent with the boundary conditions. While the base flow is periodically forced at the inlet boundary, according to equation (3.2), linear perturbations are prescribed to be zero there.

Numerical implementation

The evolution of the linearised system (4.4) is calculated using a fully implicit finitedifference time-stepping scheme of second order implemented in FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. The mesh has the same size and resolution as the one used in the Nek5000 calculations. P 2 finite elements are used for the velocity perturbation whereas P 1 finite elements are used for the pressure.

In the numerical implementation of the Floquet mode calculation, only the velocity perturbation u is considered. This is possible in incompressible flow, since the full state q = (u , p ) is fully determined by u alone. Therefore, the standard projection operator P q→u from the q-space to the u-space can be defined, as well as its inverse P u→q . The operator Φ = P q→u ΦP u→q maps a given velocity perturbation to its value after one flow period. The modal stability properties of Φ are the same as those of Φ, and Φ will be considered in what follows.

By use of a block-Arnoldi method [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems: Revised Edition[END_REF], it is possible to construct a matrix representation of Φ in a reduced orthonormal basis, generated by power iterations. A N vec -vector block-Arnoldi is iterated over N stages, each stage consisting in timeintegration of the linear flow equations over one flow period. Contrary to the standard Arnoldi algorithm, where the image of only one vector is calculated in each iteration stage, N vec vectors are advanced simultaneously in the block-Arnoldi method. A value N vec = 30 was used in all calculations presented in this section. Eigenvalues of the resulting matrix, of reduced dimension N N vec × N N vec = 750 × 750, may then be obtained. Concretely, the algorithm involves the following steps: (i)An orthonormal basis of N vec initial velocity perturbation vectors û0 i is built, with i = 1, . . . , N vec . Orthonormalisation of these vectors with respect to the energy scalar product

u, v = Ω ru(r, z) • v(r, z)drdz (4.7) is enforced. (ii)At iteration n ∈ [1, N -1], the images Φû n-1 i
of each vector ûn-1 i after one flow period are computed simultaneously by time-stepping. (iii)The Gram-Schmidt algorithm is employed to extract and normalise the component ûn i of Φ ûn-1 i for 1 i N vec that is orthogonal to the already existing set of vectors {û j i }, with 1 i N vec when 0 j n -1 and 1 i < i when j = n. Thereby, the orthonormal basis {û j i } is augmented by dimension N vec in every iteration n. (iv)In the end, after N iterations over one flow period, a fully orthonormal Krylov basis {û 0 i , . .

or Φ = R T M (Φ R) , (4.9)
with M the mass matrix associated with the scalar product (4.7).

(v)By computing the eigenvalues µ k and eigenvectors α k of Φ, defined such as dim(α k ) = N N vec , the Floquet multipliers µ k are obtained directly, and the Floquet modes v k can be reconstructed as

v k = N -1 n=0 Nvec i=1 α k nNvec+i ûn i = Rα k (4.10)
As pointed out by [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems: Revised Edition[END_REF], the orthonormalisation step (iii) is essential for the recovery of non-dominant eigenmodes in the n th iteration amidst the numerical noise on the level of round-off error.

The above algorithm is designed to maximise numerical efficiency in combination with a linear flow solver based on implicit time-stepping. As the base flow is time-dependent, a linear operator is constructed and factorised at each time step. It would not be economical to use this factorised operator for the time advancement of one single state vector; by use of the block-Arnoldi method, N vec vectors can be advanced in time simultaneously, thus lowering significantly the numerical burden of constructing a high-dimensional Krylov subspace.

A minimum of N = 25 flow-period iterations has been used to generate the following results. This number was increased in steps of 5 as necessary in order for the dominant eigenvalue to converge to four significant digits.

list of all examined flow configurations is given in table 1; the influence of the Reynolds number, the Strouhal number and the forcing amplitude may thus be characterised. For high values of St D , a larger number of block-Arnoldi iterations are required in order to achieve convergence. Two competing time scales characterise the dynamics: the forcing period, which decreases with increasing St D , and the convection time, which is constant in all cases. Therefore, a constant number of iterations at high Strouhal number corresponds to a shorter convection time.

Floquet instability modes

As shown in table 1, nearly all calculations have been performed for constant forcing amplitude A = 5%, and with fixed values of either St D = 0.6 or Re = 2000, in order to track the isolated influence of Strouhal and Reynolds number on the instability behaviour. In nearly all cases, one strictly subharmonic Floquet mode is identified, characterised by a negative real Floquet multiplier µ. This mode is observed to be unstable over certain Figure 13: Absolute value of the dominant subharmonic Floquet multiplier for different (Re, St D ), obtained for A = 5%. When the system is unstable, the most unstable mode is always subharmonic (e.g. real negative), except for the two labeled cases in (b), where the most unstable mode is almost subharmonic. When the system is stable, the represented mode are the least stable of the subhamonic domain. Black and grey bars represent parameter regimes with and without vortex pairing, according to uncontrolled DNS (see figure 7).

ranges of St D and Re, where µ falls below -1. The absolute value |µ| = -µ for St D = 0.6 and A = 5% is plotted as a function of Re in figure 13a: by linear interpolation of the critical Reynolds at which µ c = -1, instability in this case is found to arise for Re 1371. This is to be compared to the critical band Re ∈ [1350; 1400], above (resp. below) which sustained pairing was found to occur (resp. not to occur) in the DNS, as discussed in § 3.4.1. The paired and unpaired regimes, as identified in the DNS, are indicated in figure 13a by black and grey bars, in order to highlight the agreement with the onset of subharmonic Floquet instability. Results for variations in St D , at fixed values Re = 2000 and A = 5%, are presented in the same manner in figure 13b. A finite band of subharmonically unstable Strouhal numbers is identified, again in agreement with the prevalence of vortex pairing as observed in the DNS. Two values of µ reported in this diagram are distinct from the others: at St D = 0.75, the Floquet multiplier of the most unstable mode appears as a complex conjugate pair with small imaginary parts, as indicated in the figure. This mode is therefore nearly subharmonic, but not strictly so, and further iterations of the block-Arnoldi procedure do not change this result. Higher deviations from the negative real axis are found in the stable case St D = 0.9. For St D = 0.75, the unstable mode is slightly detuned, but the DNS does not display any irregular behaviour.

The effect of the forcing amplitude A on the instability is demonstrated for a single setting Re = 2000 and St D = 0.6. As shown in table 2, a higher amplitude leads to stronger instability, consistent with the DNS observations discussed in Vortex pairing was found to occur in all three configurations.

The spatial shape of an unstable subharmonic Floquet mode is presented in figure 14, for parameters Re = 2000, St D = 0.6 and A = 5%, associated with µ = -1.17. A snapshot of perturbation vorticity is shown. Its axial wavelength corresponds to twice the spacing between vortices in the unpaired base flow, and its amplitude maximum occurs far downstream, around z = 20. The latter seems rather surprising, because the base flow vortices at this position are already quite diffuse, as can be seen in figure 3c. Furthermore, vortex pairing the unstabilised flow is observed around z = 3.

In to demonstrate the effect of this modal shape onto the unpaired base flow, the are superposed, with a perturbation amplitude that is chosen ad hoc. The resulting vorticity field is shown in figure 15. It is seen that the perturbation indeed displaces the vortices around z = 20 in a fashion that indicates pairing. This pattern was also found by Selçuk et al. (2017a): the superposition of their quasistatic helical base flow and of their most unstable mode shows that their global mode shifts the helical vortices to trigger pairing.

However, in the current problem, the result from the superposition is very distinct from that of the fully developed paired state shown in figure 3a. To explain this discrepancy, it might be speculated that nonlinear adjustments could occur, or that a different, nondominant Floquet mode could be responsible for the onset of pairing. Indeed, a second unstable mode exists at this parameter setting, characterised by a Floquet multiplier value -1.03 ± 0.14i, and its vorticity distribution is given in figure 16. It displays the same spatial pattern as the first Floquet mode, but shifted several diameters further downstream, and this shift presumably accounts for its weaker growth. As this mode is not strictly subharmonic, it does not evolve precisely at ω f /2; still its growth may trigger subharmonic pairing in the nonlinear regime. However, this mode does not provide a more plausible interpretation of the observed vortex pairing much further upstream.

Transient growth analysis

In this section, we show that the transient dynamics is essential to understand the bifurcation from an unpaired unstable flow to a paired flow. Indeed, it will be demonstrated that even though the stability analysis predicts in which parameter range pairing occurs, the transition from an unpaired state arises much and much closer to the inlet that what modal theory predicts. In contrast, transient dynamics predicts more accurately the transition rate and the perturbation structures, and it allows to explain the occurrence of intermittent phenomena in sub-critical but noisy jets shown in § 3.4.3.

Growth of random initial perturbations

The jet at Re = 2000, forced at St D = 0.6 with 5% amplitude is considered throughout this section. The natural state in this case is the paired one, shown in figure 7, consistent with an unstable Floquet multiplier µ = -1.17 as discussed in § 4.4. Direct numerical simulation results are presented here, which aim to show how pairing is triggered in the unstable unpaired flow.

A first simulation is performed starting from the stabilised unpaired state, displayed in figure 3c, as initial condition. Residual non-T -periodic components (u(t)u(t -T ))/2 in this flow state are of the order of 0.01% of the reference jet velocity.

Two additional simulations have been carried out with the same state, but with added white noise velocity perturbations, u(r, z, t = 0) = U u (r, z, t = 0) + (r, z). This noise exhibits zero spatial mean = 0, and two different standard deviations 2 = (10 -3 , 10 -3 ) T and (10 -4 , 10 -4 ) T are prescribed in the two simulations. Non-harmonic components of the flow state at any given time are measured by a norm defined as

e(t) = 1 2 u(t) -u(t -T ) 2 r dr dz. (5.1)
The time evolution of this norm is traced in figure 17 for the three different initial conditions. While all three cases evolve into the same paired attractor state, they arrive there along different trajectories.

Modal growth from very low perturbation amplitude

Starting from the unpaired state, as it has been obtained through flow stabilisation, without any added random noise, the initial perturbation that may give rise to a growing subharmonic component is given by the residual non-T -periodicity that remained when the stabilised calculation was halted. The total vorticity after one period T is shown in figure 18a, and the magnitude of the non-T -periodic residual at the same instant is presented in figure 19a. Note that this residual is computed by comparing the state at time t with that from time t -T ; therefore time-stepping over one period is required before the stabilisation effect can be evaluated.

10 0 • • • • • • • • • • t T e(t)
With initial perturbation of 10 -3 With initial perturbation of 10 -4 Without initial perturbation Modal growth rate (Floquet multiplier) Figure 17: Growth of non-harmonic components in simulations with and without added initial noise. The modal Floquet growth rate indicated for comparison.

Over the first forty periods, the subharmonic perturbation is dominated by a slow growth of the residual structure, but at such small amplitude that the total vorticity field in figures 18b and 18c is not noticeably affected. Then, as can be seen by comparing figures 19c and 19d with figure 14, the exponentially growing unstable Floquet mode becomes manifest. Its amplitude growth per flow period, between = 40T and 60T , is estimated from figure 17 as a factor 1.14, to be compared to the absolute value 1.17 of the computed Floquet multiplier. The spatial structure of this perturbation, displayed in figure 20, exhibits a similar structure as the one of the corresponding Floquet mode, shown in figure 14. Beyond = 60T , a nonlinear saturation of the perturbation sets in (figure 17), accompanied by a change in its spatial shape. At t = 90T , as the flow approaches the asymptotic periodic regime, the maximum perturbation amplitude has moved upstream to z ≈ 10 (figure 19e), where a pronounced pairing of vortices is observed in figure 18e. This pairing location still moves further upstream with time, until it will finally stabilise near z = 3, in the natural asymptotic paired state shown in figure 3a.

Non-modal growth from initial white noise

The simulations with added white noise in the initial condition (solid lines in figure 17) show a much faster convergence to the final paired flow state than the case discussed in the preceding section (dotted line in figure 17). In particular, the initial growth of e(t) in these two cases is markedly stronger than that of the unstable Floquet mode.

The time development from an initial condition with 2 = (10 -3 , 10 -3 ) T (thick solid line in figure 17) is visualised in figure 21 by successive snapshots of the total vorticity. It is seen that vortex pairing not only sets in faster than in the previous case of figure 18, but also much further upstream. Corresponding non-periodic perturbations are again displayed in figure 22: the smallest scales of the random initial condition are quickly dissipated (figure 22a), and a growing coherent perturbation structure is evident after a few period cycles (figure 22b). The maximum perturbation growth in this phase, as measured from figure 17, corresponds to a factor 2.01 per period, much stronger modal growth factor 1.17. Persistent vortex pairing is fully established at = 10T (figure 22c); subsequently, the pairing location slowly moves upstream, and stabilises around z = 3.

0 1 2 r (a) t = T 0 1 2 r (b) t = 20T 0 1 2 r (c) t = 40T 0 1 r (d) t = 60T
Vorticity perturbations after five flow periods, in these simulations with added initial noise, are represented in figure 23. At this early stage, their dynamics may still be regarded as linear; however, the perturbations are now located close to the inlet, and their spatial distribution bears no resemblance with the unstable Floquet mode (figure 14). Therefore, a non-modal mechanism is expected to underpin this growth.

Simulations with a lower initial noise level of 10 -4 show a similar behaviour, although slightly delayed (figure 17). The exponential phase is longer, since the amplitude takes more time to saturate, and the maximum growth rate is 2.16 per forcing period.

In both the "unperturbed" and "randomly perturbed cases", one should also note the absence of any sustained leapfrogging or overtaking events before pairing, contrary to what was observed by Selçuk et al. (2017a). This can be understood because the jet rings are much thicker compared their helical vortices, making leapfrogging and overtaking difficult.

Optimal linear perturbation growth over one cycle

In order to further analyse and understand the mechanism behind the non-modal onset of vortex pairing, as observed in § 5.1.2, the optimal growth of subharmonic perturbations is now investigated. One cycle period T is chosen as the time horizon over which optimisation is performed. As discussed in the context of (4.6), perturbations are propagated over one cycle by the time-shift operator Φ . The optimal perturbation is then found as the solution to the maximisation problem

u opt (T ) = max u (t=0) u (t = T ) u (t = 0) = max u (t=0) Φ u (t = 0) u (t = 0) , (5.2)
The norm used in the following is derived from the standard real u-scalar product in cylindrical coordinates (equation (4.7)). This scalar product defines a full norm for u, but only a semi-norm for q, because the separation condition is not fulfilled. The solution of the maximisation problem (5.2) is given by the norm of the operator Φ .

To evaluate this norm, the orthonormal basis {û n i } calculated in the context of modal analysis ( § 4.3) is once more exploited. Instead of maximising the norm of Φ , a infinitedimensional operator, we maximise the norm of its projection Φ, of finite dimension, onto this basis.

A perturbation state u is projected onto {û n i } as L. Shaabani-Ardali, D. Sipp and L. Lesshafft

β nNvec+i = u , ûn i such that β = R T M u (5.4)
The residual r is orthogonal to the basis {û n i }, such that u 2 = β 2 + r 2 . Therefore,

Φ u = N -1 n=0 Nvec i=1 β nNvec+i (Φ ûn i ) + Φ r (5.5) = N -1 n=0 Nvec i=1 Φβ nNvec+i ûn i + Φ r (5.6)
and

Φ u 2 = Φβ 2 + Φ r 2 + 2 N -1 n=0 Nvec i=1 Φβ nNvec+i ûn i , Φ r , (5.7)
with the norm of Φβ calculated in a finite-dimensional space of dimension N N vec . Therefore, maximising Φ u is equivalent to maximising the right-hand side of the previous equation. However, the ûn i -basis has not been chosen randomly: being constructed from the successive iterations of a single group of random vectors, it selects numerically the fastest-growing modes of the full system, in a similar way as power iterations [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems: Revised Edition[END_REF]). These modes are then gradually excluded from the residual space of r , and the norm of the image of r through Φ decreases as the number N of Krylov subspace iterations is increased. On the right-hand side of equation (5.7), the first term becomes dominant as N increases ; the second and third terms are bounded by Φ r 2 and Φβ Φ r , respectively. Therefore, the approximation Φ u ≈ Φβ is valid for large N . This explains why, for a given value of N N vec , a trade-off needs to be found between N and N vec : N must be large enough to capture the salient flow dynamics, whereas N vec must be sufficiently large to make the block-Arnoldi calculations computationally efficient.

The norm Φβ is evaluated by use of the singular value decomposition (SVD) (5.8) with Σ a real positive diagonal matrix, and Ũ and Ṽ real unitary matrices. Columns ṽk and ũk represent forcing and response pairs in the orthonormal basis of the {û n i }. Σ contains the singular values ordered in descending order, such that Φṽ k = σ k ũk with σ k σ k+1 0.

Φ = Ũ Σ Ṽ T ,
Note that the Krylov base constructed in § 4.3 is orthonormal with respect to the scalar product (4.7). Therefore, the optimal initial perturbations or responses in the reduced space correspond to the optimal initial perturbations or responses in the fullstate space, using the R matrix to change basis. For a given optimal initial perturbation (resp. response) ṽk (resp. ũk ) in this reduced basis, the corresponding full-state initial perturbation velocity (resp. response velocity) v k (resp. u k ) is obtained as v k = Rṽ k (resp. u k = Rũ k ). The corresponding matrices of optimal initial perturbations and responses are V = R Ṽ and U = R Ũ . Because of this preserved optimality, the maximum linear gain achievable over one period in the full-state is still σ 1 , obtained by perturbing the velocity field at t = 0 with v 1 . After one period, the perturbation velocity field is given by u 1 .

This method is computationally efficient, as it is entirely based on the results already available from the modal analysis. This is in contrast with the additional computations required in a direct-adjoint approach. Moreover, with the direct-adjoint approach, the optimisation is carried out for a single time horizon T f . In this section, only T f = T has been chosen ; however, as will be seen in the next section, this can be extended with little additional effort to any time-horizon of integer periods T f = nT . Time-horizons that are Figure 24: Norm of the time-shift operator for different Reynolds and Strouhal numbers, with a 5% forcing level and a normalised momentum thickness of 5%. The grey or black domains represent the parameters values for which DNS calculations show that the "natural" solution is unpaired or paired, as shown in figure 7. In figure 24b, the leading amplification rates are also depicted rescaled over one dimensionless time unit. Additional singular values for lower Reynolds number flows are discussed in § 5.4.

Forcing amplitude 1% 5% 10% uopt(T ) 3.96 3.40 3.32

Table 3: Evolution of the optimal perturbation gain over one period with the forcing amplitude for Re = 2000 and St D = 0.60.

not integer multiples of T , corresponding to transient phenomena within a forcing period, are not considered here. Singular value decomposition of Φ has been carried out for all flow configurations listed in table 1. The maximum gain values σ 1 are shown in figure 24 for fixed values of St D and Re. Variations of σ 1 with the forcing amplitude A are given in table 3 for one setting Re = 2000 and St D = 0.6. In all cases, the achievable transient amplification over one period is significantly higher than the maximum modal growth rate. A comparison with results of perturbed DNS flows (figure 17) is discussed in § 5.3.

At fixed Strouhal number (figure 24a), the leading singular value changes weakly as a function of Reynolds number, even as the system goes from stable to unstable. Variations of σ 1 are within 20% as the value of Re is doubled.

At fixed Reynolds number (figure 24b), variations of the Strouhal number also do not affect σ 1 in a strong way. Doubling St D from 0.45 to 0.9 is accompanied by a 10% decrease of the maximum gain over one flow period. However, when σ 1 is rescaled to give the mean amplification over a constant time unit, as σ1 = σ 1 St D , this rescaled gain increases by more than 30% over the investigated interval of St D . This is consistent with [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF], 1996) : in their parametric study, as shown in figure 2, they found that when the Strouhal number is increased, the jet is more prone to experience non-T -or non-2T -periodic behaviour (periodic or not modulations, intermittency or chaos). Our study shows that the larger the Strouhal number, the more amplified random perturbations can be over a given time unit. With increasing forcing amplitude, the maximum gain is moderately diminished (table 3). This may indicate that the receptivity to subharmonic perturbations is larger when the shear layer is not yet fully rolled up, since increasing forcing amplitude results in a faster roll-up of vortices, as shown in figure 9. Again, this is coherent with [START_REF] Broze | Nonlinear dynamics of forced transitional jets: periodic and chaotic attractors[END_REF], 1996): as shown in figure 2, with larger forcing levels, except at Strouhal numbers greater than those considered in our study, the flow is less prone to non-periodic behaviour, which may be brought about by transient amplification of small disturbances.

The leading singular mode may not be the only relevant way to trigger transient growth. For the case (Re = 2000, St D = 0.6), the 20 largest singular values are represented in figure 25. Indeed, a strong dominance of the optimal perturbation cannot be affirmed: the first value is 15% larger than the second, 27% larger than the third, and only 87% larger than the twentieth.

The shapes of the optimal initial perturbation and response structures are shown for the case (Re = 2000, St D = 0.6) in figures 26a and 26b, respectively. In comparison with the Floquet mode discussed in § 4.4 (see figure 14), these perturbation structures are localised much closer to the inlet. Over the course of one period, the perturbation shape is largely conserved, while it is convected at the same pace as the rolled-up vortices. The manifestation of these perturbations in the total vorticity field is visualised in figure 27, where the linear optimal initial perturbation and response structures are superposed onto the periodic base state with a small amplitude of 5%. It can be observed that the optimal perturbation structures displace the vortices towards and away from the axis, in an alternating fashion. This displacement is furthermore oriented at an oblique angle, such that the distance between neighbouring vortices is modulated, thus initiating the pairing interaction.

The second singular mode pair in the same configuration, corresponding to σ 2 , is displayed in figure 28. It exhibits a very similar structure as the first singular mode, the main difference being that the maximum amplitude is shifted downstream to the next vortex pair.

Optimal linear perturbation growth over many cycles

technique that has been described in the previous section for optimisation over a single flow period T is easily extended in order to construct optimal perturbations for time horizons of multiple periods nT . Let f n denote the optimal perturbation at t = 0 that leads to the largest possible flow response r n at t = nT . The associated amplitude gain is then given by the norm u opt (nT ) = r n / f n , which is identical to the largest singular value of the n th power of the operator Φ. Substituting Φ with its approximation Φ, the SVD of Φn is straightforward to compute. The ensuing variation of the optimal gain with nT is traced in figure 29, alongside the perturbation growth associated with the most unstable Floquet mode ( § 4.4). Optimal perturbation at t = 0 enables a vigorous transient growth over short time horizons; at longer times, however, the optimal growth rate approaches that of the modal solution, given by the modulus of the dominant Floquet multiplier. This behaviour is identical to that of perturbation growth in steady base flows, which has been discussed theoretically by [START_REF] Trefethen | Spectra and pseudospectra: the behavior of nonnormal matrices and operators[END_REF]. In the present case, the transient non-modal growth provides a significant boost, on the order of 10 5 , of the overall long-time amplitude gain. This additional factor corresponds to 73 cycle periods of modal growth. The evolution and convergence of the optimal initial perturbation and response structures are examined in figure 30 in terms of projections. The scalar products of the initial perturbation and response structures obtained for any number of cycles with those at n = 1 and n = 50 are represented. It is observed that the optimal initial perturbation structure does not evolve much with time, and that the shape determined for one single period is close to optimality for all time horizons. In contrast, the shape of the optimal response changes significantly through time; the optimal response at t = 50T is indeed localised at a different location than its counterpart at t = T . While the optimal initial perturbation structure always retains a nearly identical shape, the optimal response structure gradually shifts downstream with increasing time horizon, and eventually it

product | f n , f 1 | | f n , f 50 | | r n , r 1 | | r n , r 50 |
Figure 30: Projection of the optimal initial perturbation and response at time nT over the optimal initial perturbation and responses at times T and 50T for a jet at Re = 2000 and St D = 0.60 forced with a 5% amplitude. While the optimal initial perturbation experiences little change, the optimal response changes drastically. All vectors considered have unit norm.

| f1, ṽ1 | | f50, ṽ1 | | r1, ṽ1 | | r50, ṽ1 |
1.24 × 10 -3 1.10 × -5 3.60 × 10 -3 9.9990 × 10 -1

Table 4: Scalar products (4.7) between the Floquet mode ṽ1 and the optimal finite-time initial perturbations f n and responses r n . All modes and perturbations have unit norm. Consistent with the notations introduced in § 5.2, f 1 is identical to v 1 , and r 1 to u 1 .

converges towards the structure of the Floquet mode (see figure 14). The projection values in table 4 demonstrate this convergence. It is expected on theoretical grounds [START_REF] Trefethen | Spectra and pseudospectra: the behavior of nonnormal matrices and operators[END_REF]) that, at long time horizons, the optimal response structure corresponds to the dominant Floquet mode, whereas the optimal forcing is given by the associated adjoint mode of the time-shift operator.

The different evolution from different initial perturbations, documented in § 5.1.1 and 5.1.2 can now be interpreted in terms of the optimal perturbation results.

In the two cases that were initialised with white noise perturbations, this initial condition contained significant components of strongly amplified optimal initial perturbation modes. After approximately five periods, these structures emerged from the background noise, as may be inferred from a qualitative comparison between figures 23, 26b and 28b. However, the maximum growth observed in figure 17 for 2 = (10 -3 , 10 -3 ) T corresponds to a factor 2.16 per period, whereas the leading singular value was determined to be 3.40 in § 5.2 for the time horizon t = T . Yet, the optimal growth rate decreases with time; this is shown in figure 31, which represents the marginal gain σ m (nT ) at t = nT of the nT -optimal perturbation: Table 5: For St D = 0.6 and A = 0.05, evolution of the least stable subharmonic Floquet multiplier and of the maximum growth rate over one period as a function of Re in the stable case.

served in figure 17, is consistent with the marginal gains observed after 5 to 10 periods in figure 31.

In the case discussed in § 5.1.1, initialised without added noise, the exponential growth observed between 40T and 60T has been associated with the unstable Floquet mode. The initial perturbation (figure 19a) is spatially disjoint from optimal perturbations near the inflow. At t = 20T , a non-modal pairing perturbation is faintly visible in figure 19b, still in the process of transitionning towards the Floquet mode in figure 19c, the optimal linear response at long times.

Whether or not transient processes lead to a rapid onset of nonlinear vortex pairing depends on the energy level of oscillatory perturbations near the inflow, with appropriate wavelengths for subharmonic vortex modulation.

Importance of transient dynamics at low Re

In § 3.4.1, a Reynolds number threshold was determined below which pairing did no longer occur naturally. However, in § 3.4.3, it has been demonstrated that forcing a jet harmonically with a small noise level was able to trigger pairing for Reynolds numbers where pairing does not normally occur. It is now examined how this behaviour can be explained in the light of the transient growth analysis developed in the previous section. At St D = 0.6 and A = 5%, four Reynolds number values 1300, 1000, 750 and 500 are considered here. All these cases are modally stable, as seen in table 5. The transient growth of these cases is shown in figure 32.

Although in all cases, the gain decreases in the long term, the maximum achieved during the transient growth varies largely, between 4.5 after 6 periods for Re = 500, and 8.3 × 10 3 after 26 periods for Re = 1300. This change explains the different behaviour observed in § 3.4.3. At Re = 500, whatever the forcing level, the noise is not energetic enough to trigger pairing. On the other hand, at Re = 1300, even very small perturbations eventually initiate pairing via transient growth. At intermediate Reynolds numbers, the behaviour observed depends on the perturbation level.

Conclusion

The onset of sustained vortex pairing in the street of axisymmetric vortices in a laminar jet has been investigated by means of global instability analysis.

The numerical framework developed for this study is applicable to the instability analysis of any spatially developing time-periodic flow. It includes the computation of a strictly T -periodic base flow by means of DNS with added time-delay control (Shaabani-Ardali et al. 2017), the identification of dominant linear Floquet modes, and a singular value decomposition of the propagator in order to characterise non-modal transient growth phenomena. The modal and non-modal analysis are both achieved via timestepping of the linearised flow equations, without the need for adjoint calculations, in combination with a block-Arnoldi algorithm [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems: Revised Edition[END_REF].

It has been demonstrated that self-sustained vortex pairing in a jet is the manifestation of a subharmonic linear Floquet instability of the underlying periodic vortex street. Direct numerical simulations of axisymmetric laminar jets, forced harmonically at Strouhal number St D at the inlet, have been perfomed in order to delineate the region in the Re-St D parameter plane where pairing in the ensuing vortex street arises spontaneously. This parameter regime has then been shown to be characterised by subharmonic Floquet instability. Furthermore, the phase velocities of fundamental and subharmonic fluctuations in the simulation results have been found to be consistent with the resonance model proposed by [START_REF] Monkewitz | Subharmonic resonance, pairing and shredding in the mixing layer[END_REF]. DNS results also show a strong influence of inflow forcing amplitude as well as random ambient noise on the onset of vortex pairing.

The spatial structure of the unstable linear Floquet mode, throughout the parameter regime considered in this study, reaches its maximum amplitude far downstream of the inlet, whereas the subharmonic perturbation amplitude in the nonlinear saturated paired flow state peaks a few diameters away from the inlet, where stationary vortex pairing occurs. When DNS calculations in the unstable regime are initiated with extremely low subharmonic perturbations, the Floquet mode structure and growth rate are indeed observed over a short time interval in the bifurcation process. However, when additional random noise is added to the initial condition, considerably faster subharmonic perturbation growth is observed, with spatial amplitudes concentrated near the inlet. This behaviour motivated a linear optimal perturbation analysis. It has been found that linear transient mechanisms may induce, for a standard configuration with Re = 2000 and St D = 0.6, an additional amplitude gain of about five orders of magnitude with respect to purely modal growth. Furthermore, the spatial structure of the linear optimal perturbation response reflects the DNS observations at early times, when subharmonic perturbations may be assumed to be governed by linear dynamics.

While the asymptotic occurrence of vortex pairing is determined by modal Floquet instability, its emergence is dominated by non-modal transient growth mechanisms.

Figure 3 :

 3 Figure 3: Vorticity snapshots of the periodic paired and unpaired states, obtained naturally for two different parameter settings. Forcing amplitude, Reynolds and Strouhal numbers are defined in § 3.1.

Figure 4 :Figure 5 :

 45 Figure 4: Vorticity snapshots of a paired case at Re = 2000 forced harmonically at St D = 0.60 and A = 5%. Only the region near the inlet is shown. The colour coding is the same as in figure 3c.

Figure 6 :

 6 Figure 6: Phase speed along the line of maximum subharmonic mode of the subharmonic and fundamental modes of a paired jet at (Re = 2000, St D = 0.60), and comparison with the phase speed of the fundamental mode along the lime of maximum fundamental mode of the corresponding stabilised (unpaired) jet. The phase of the Fourier mode is denoted as ϕ.

Figure 7 :

 7 Figure 7: Occurrence of vortex pairing as a function of Reynolds number and Strouhal number for configurations with forcing amplitude A = 5% and initial shear layer thickness θ 0 = 0.025. Only points near the pairing boundary are represented.

Figure 8 :

 8 Figure 8: Occurrence of vortex pairing as a function of the Reynolds and the Strouhal number for two different forcing amplitudes A = 1% and 10%. The same symbols as in figure 7 are used.

Figure 9 :

 9 Figure 9: Total energy of the fundamental and subharmonic modes in each plane z = const as a function of z at St D = 0.6 and Re = 2000. Three values of the inlet forcing amplitude are considered, A = 1%, 5% and 10%.

Figure 10 :

 10 Figure 10: Vorticity snapshots at Re = 750, St D = 0.60 and A = 5% with A = 1.0% taken at two different times, and both at the same phase. Pairing is intermittent, occurring in (a) and not in (b). The colorbar is rescaled compared to figure 3b, to magnify the behaviour of downstream vortices.

Figure 11 :

 11 Figure 11: Vorticity snapshots at Re = 1300 and St D = 0.60 with A = 0.1% taken at two different times, and both at the same phase. This shows the pairing location fluctuation, emphasised in grey, in this setup. The colorbar is in figure 10.

Figure 12 :

 12 Figure 12: Vorticity snapshots at Re = 1000 and St D = 0.60 with A = 1.0%: (b) is taken two forcing periods after (a). The three vortices located in the shaded are on (a) merge altogether, as an example of collective interaction. The colorbar is in figure 10.

  .2) to the right-hand side of the governing equations (3.1a)-(3.1c). Such forcing attenuates all 2T -periodic fluctuations, but leaves T -periodic components unaffected. The simulations converge towards a purely T -periodic state, in which the artificial forcing term vanishes.

Figure 14 :Figure 15 :

 1415 Figure 14: Vorticity component of the most unstable Floquet mode at Re = 2000 and St D = 0.60. It peaks far downstream.

Figure 16 :

 16 Figure 16: Real part of the vorticity component of the second most unstable Floquet mode at Re = 2000 and St D = 0.6. Its peak is located far downstream.

Figure 18 :

 18 Figure 18: Vorticity snapshots at five different instances without added initial noise. The colour coding is the same as in figure 3b.

Figure 20 :Figure 21 :

 2021 Figure 20: Vorticity field of the perturbation without any forcing applied at t = 40T .

Figure 22 :Figure 23 :

 2223 Figure 22: Perturbation magnitude u(t)u(t -T ) at five different instances evolving out of initial white noise level 10 -3 . 0 5 10 15 20 25 30 35 40 0 1 2 z

Figure 25 :

 25 Figure 25: The twenty first singular values of Φ for Re = 2000, St D = 0.6 and A = 5%.

Figure 26 :Figure 27 :

 2627 Figure 26: For Re = 2000 and St D = 0.6, vorticity of (a) the optimal perturbation and (b) its response after one flow period. The maxima are located close to the inlet.

Figure 28 :

 28 Figure 28: For Re = 2000 and St D = 0.6, vorticity of (a) the second singular mode: (a) perturbation and (b) its response after one flow period. The maxima are located further downstream than in the optimal setting shown in figure 26.

Figure 29 :

 29 Figure 29: Transient growth gain achieved by a (Re = 2000, St D = 0.60) jet forced with 5% amplitude and comparison with modal growth.

Figure 31 :

 31 Figure 31: Marginal gain σ m (nT ) of the optimal perturbation at the horizon time nT for Re = 2000, St D = 0.60 and A = 0.05. The marginal growth rate asymptotically converges towards the modulus of the most unstable Floquet multiplier.

Figure

  Figure Transient growth gain for a jet forced at St D = 0.60 with 5% amplitude at different subcritical Reynolds number values.

Table 1 :

 1 Parameter combinations for which Floquet analysis is performed.

	represented by the matrix
	ΦnNvec+i,n

. , ûN-1 i } and their images {Φ û0 i , . . . , Φ ûN-1 i } after one period of time-stepping are obtained. Let R denote the matrix of this basis, R = û0 1 , . . . , û0 Nvec , . . . , ûN-1 1 , . . . , ûN-1 Nvec , (4.8) which is of dimension N dof × N N vec , with N dof the number of degrees of freedom of the initial velocity perturbation. It is then possible to construct the projection Φ of the infinite-dimensional operator Φ onto the finite-dimensional space spanned by R. Φ is Nvec+i = Φ ûn i , ûn i

Acknowledgements

We thank Holly Johnson for her help with Nek5000 and Robin Yegavian for his help with the parallelisation of FreeFem++.

This work was granted access to the HPC resources of TGCC under the allocation Grant No. 2016-2a6451 made by GENCI, and it was financially supported by the Délégation Générale de l'Armement under Grant No. 2015.60.0004.00.470.75.01.