Samir F Matar
email: matar@icmcb-bordeaux.cnrs.fr

Adel F Al Alam

Naim Ouaini

Hydrogen Insertion Effects on the Electronic Structure of Equiatomic MgNi Traced by ab initio Calculations

Keywords: Hydrogen, Equiatomic MgNi, Iono-covalent Character, Bader Charge, Chemical Bonding

For equiatomic MgNi which can be hydrogenated up to the composition MgNiH 1.6 at an absorption/desorption temperature of 200 °C, the effects of hydrogen are approached with model structures MgNiH, MgNiH 2 and MgNiH 3 . From full geometry optimization and calculated cohesive energies obtained within DFT, the MgNiH 2 composition close to the experimental limit is identified as most stable. Charge density analysis shows increasingly covalent character of hydrogen: MgNiH (H -0.67 ) → MgNiH 2 (H -0.63 ) → MgNiH 3 (H -0.55 ).

While Mg-Ni bonding prevails in MgNi and hydrogenated model phases, extra itinerant low energy lying Ni states appear when hydrogen is introduced signaling Ni-H bonding which prevails over Mg-H as evidenced from total energy calculations and chemical bonding analyses.

Introduction

In spite of the large hydrogen sorption capacity (7.6 wt.-%), the strong metal-hydrogen bonding in ionic like MgH 2 prevents its use, and other Mg based compounds were studied as candidates for hydrogen storage. It was found early that introducing nickel leads to materials with reversible adsorption ability such as in Mg 2 NiH 4 [1].

The Mg-Ni system comprises Mg 2 Ni and MgNi 2 , but no report on the equiatomic composition has been published [2]. Recently, Kamata et al. [3] obtained equiatomic MgNi under HP/HT (5 GPa / 300 °C) conditions. The structure is tetragonal with CuTi-type and small distortion with respect to cubic CsCl. The equiatomic was readily hydrogenated into MgNiH 1.6 and found to absorb and desorb hydrogen reversibly at 200 °C [3]. This result characterizes it as a hydride with 1.9 wt.-% H capacity, smaller than the rate of Mg 2 NiH 4 (~3.6 wt.-% H).

However, from earlier work on amorphous nanostructured Mg-Ni a higher H content of 2.2 wt.-% in the equiatomic region was reported, i. e. with the composition MgNiH 1.84 , close to MgNiH 2 [4]. Yet this amount of maximum hydrogen uptake is the lowest among nickel based binary equiatomics such as YNi which can absorb up to 4 H per formula unit (FU) [5].

In the context of our continuing investigations of candidate compositions for hydrogen storage, focusing on the electronic structure of hydrogenated nickel based intermetallics (cf.

refs. [5,6] and works cited therein), we address the bonding properties in the equiatomic MgNi and its hydrogenated compound based on computations within the well established quantum theoretical density functional (DFT) framework [7].

Computation Methods

Two computational methods within the DFT were used in a complementary manner. The Vienna ab initio simulation package (VASP) code [8,9] allows geometry optimization and total energy calculations. For this we use the projector augmented wave (PAW) method [9,10], built within the generalized gradient approximation (GGA) scheme following Perdew, Burke and Ernzerhof (PBE) [11]. Also semi-core 2p-Mg states were accounted for beside 3s valence states in the PAW potential construction. Preliminary calculations with local density approximation LDA [12] led to a largely underestimated volume versus the experiment. The conjugate-gradient algorithm [START_REF] Press | Numerical Recipes[END_REF] is used in this computational scheme to relax the atoms.

The tetrahedron method with Blöchl corrections [10] as well as a Methfessel-Paxton [START_REF] Methfessel | [END_REF] scheme were applied for both geometry relaxation and total energy calculations. Brillouinzone (BZ) integrals were approximated using special k-point sampling. The optimization of the structural parameters was performed until the forces on the atoms were less than 0.02 eV/Å and all stress components less than 0.003 eV/Å 3 . The scalar relativistic calculations were converged at an energy cut-off of 400 eV for the plane-wave basis set with respect to the k-point integration with a starting mesh of 4  4  4 up to 8  8  8 for best convergence and relaxation to zero strains.

Then all-electron calculations with GGA were carried out for a full description of the electronic structure and the properties of chemical bonding, using full potential scalarrelativistic augmented spherical wave (ASW) method [15,[START_REF] Eyert | The Augmented Spherical Wave Method -A Comprehensive Treatment[END_REF]. In the minimal ASW basis set, we chose the outermost shells to represent the valence states and the matrix elements were constructed using partial waves up to l max + 1 = 3 for Mg and Ni, and l max + 1 = 1 for H. Selfconsistency was achieved when charge transfers and energy changes between two successive cycles were below 10 -8 and 10 -6 eV, respectively. BZ integrations were performed using the linear tetrahedron method within the irreducible wedge. In order to optimize the basis set, additional augmented spherical waves were placed at carefully selected interstitial sites (IS).

Besides the site projected density of states, we discuss qualitatively the pair interactions based on the overlap population analysis with the crystal orbital overlap population (COOP) [17]. In the plots, positive, negative, and zero COOP indicate bonding, anti-bonding, and non-bonding interactions, respectively.

Geometry Optimization, Cohesive Energies and Charge Transfers

MgNi is found to crystallize in the tetragonal CuTi type with space group P4/mmm and a = 2.997, c = 3.166 Å and V = 28.44 Å 3 per FU which can be considered as a tetragonally distorted CsCl cubic structure [3]. Upon hydrogenation the structure becomes cubic, Pm3

m with a = 3.215 Å and V = 33.23 Å 3 per FU and a volume increase of V/V ~ 17 %.

Inserting hydrogen in the cubic structure leads to a random occupation of the threefold 3c positions within Pm3 m. However, this is difficult to model in computational situations in which one can assume the occupation of only two out of the three (3c) positions in order to approach the experimental stoichiometry MgNiH 1.6 . Due to symmetry constraints the (3c)

positions split into (1c) ½ ½ 0 (model structure MgNiH with H2 in Table 1 and Fig. 1) and (2e) 0 ½ ½; ½ 0 ½ (MgNiH 2 with H1 in Table 1 and Fig. 1). The structure is then described within P4/mmm (subgroup of Pm3 m). On the opposite, the full occupancy of the 3c sites such as in the MgNiH 3 model structure keeps the cubic symmetry. In order to appreciate the distortion effects and to follow the progressive hydrogen introduction, the model compositions MgNiH, MgNiH 2 and MgNiH 3 besides MgNi have been examined. In the intermetallic CsCl-related structure the larger atom, i. e. Mg (r = 1.60 Å), is expected to be found at the origin (0 0 0), and the smaller one, i. e. Ni (r = 1.24 Å) at body center (½ ½ ½) position. The two positions are equivalent as long as no neighbors are introduced, because the CsCl structure is composed of two interpenetrating simple cubic substructures. In fact the calculations with straight and inverted structures in both cubic and tetragonal MgNi led to identical volumes and energies. The results of the geometry optimization are given in Table 1, 2 nd column; they show that the calculated volume of 28.32 Å 3 per FU is in the range of the experimental value (V exp. = 28.44 Å 3 per FU). Close total energy magnitudes between the tetragonal and cubic structures are found, with the former identified as the ground state structure due to slightly lower energy, in agreement with experiment.

The introduction of hydrogen is modeled with the structures sketched in The cohesive energies are obtained from energy differences between the total energy of the compound and those of the atomic constituents Mg and Ni and dihydrogen which have the following energies from PAW-GGA calculations (in eV units per atom): Mg (-1.543), Ni (-5.539) and H 2 (-6.577). The resulting cohesive energies in eV/atom are given in Table 1:

MgNi (-0.440) → MgNiH (-0.226) → MgNiH 2 (-0.247) → MgNiH 3 (-0.223)

This trend can be interpreted as a destabilization of the intermetallic upon the introduction of hydrogen, but the least unstable composition is the one with 2H. This resembles the experimental data where the observed MgNiH 1.6 composition [3] is the closest to our model MgNiH 2 which has the largest cohesive energy among the three model structures.

Focusing on MgNiH 2 the other twofold insertion sites 2f and 2g were examined for the sake of completeness. Full geometry optimization led respectively to total electronic energies:

-12.22 eV and -11.62 eV (with calculated z = 0.605), i. e. to less stable structures than MgNiH 2 with H at (2e) positions (Table 1).

The analysis of the charge density issued from the self consistent calculations can be done using the AIM (atoms in molecules theory) approach [18] developed by Bader who devised an intuitive way of splitting molecules into atoms as based purely on the electronic charge density. Typically in chemical compounds, the charge density reaches a minimum between atoms and this is a natural region to separate them from each other. Such an analysis can be useful when trends between similar compounds are examined; it does not constitute a tool for evaluating absolute ionizations. Bader's analysis is done using a fast algorithm operating on a charge density grid. The program [19] yields charge densities obtained from high precision VASP calculations and generates the total charge associated with each atom.

The results of computed charge changes (Q) are such that they lead to neutrality when the respective multiplicities are accounted for; the values are given for each atomic species in Table 1. The covalent character brought into the electronic structure is assigned mainly to nickel owing to its binding with hydrogen, as illustrated below.

Electronic Density of States and Chemical Bonding Properties

Starting from the calculated crystal structure data in Table 1, the electronic density of states and the chemical bonding were analyzed for the four studied compounds. At self consistent convergence, charge transfers follow the above discussed trends with small residues ~0.15 electron into the interstitial spheres IS.

Fig. 1 .

 1 They are characterized by Ni-H distances at ~a/2 shorter than Mg-H at ~a2/2 and Mg-Ni at ~a3/2, considering a cube with an a edge for simplicity. Also from the positions of the hydrogen in the model MgNiH and MgNiH 2 the tetragonality ratios are expected as: c/a > 1 (elongation) and c/a < 1 (flattening), respectively. The results shown in Table 1 confirm the trends for the distances (d Mg-H > d Ni-H ) and the tetragonality c/a ratios close to 1, i. e. to a cubic structure. The volume increases following the series with additional inserted hydrogen atoms. V calc. (MgNiH 2 ) = 34.4 Å 3 is slightly larger than the experimental value of 33.23 Å 3 due to 2 H instead of 1.6, but the agreement allows to assume a correct modeling of the experiment. In spite of starting with the assumption tetragonal structure with H1 and H2 at different hydrogen sites in model MgNiH 3 , the calculations result in the cubic structure for the ground state which resembles the perovskite structure.Nevertheless, for the sake of completeness, we further hypothesized the interchange of Mg and Ni positions in hydrogenated MgNi. This resulted in a substantial rise (destabilization) of the total energies, shown with 'INV' labels in Table1, as a function of the amount of introduced hydrogen. The energy differences per atom in eV/FU units are then: E(MgNiH) = 0.13, E(MgNiH 2 ) = 0.58 and E(MgNiH 3 ) = 0.79, i. e. the higher the hydrogen content, the larger the destabilization. This confirms the stable configuration of the Mg (0 0 0) and Ni (½ ½ ½) positions with shorter Ni-H versus Mg-H distances, and the trend of hydrogen to bind with nickel preferentially.

  Due to the large electronegativity difference between Mg ( = 1.31) and Ni ( = 1.91), a charge transfer is expected from Mg to Ni as observed for MgNi with a charge transfer of ±1.61. In view of such large charge transfer featuring a partially ionic character, one should rather assume that s-like charges departing from Mg are close to free electron as in normal metals and delocalized over the whole MgNi crystal. The trends for the hydrides are such that in MgNiH the charge transfer is from Ni to H resulting in H -0.67 . This iono-covalent (in between -1, ionic, and ~ < 0, covalent) character of hydrogen is shifted to more covalent in MgNiH 2 with H -0.63 and then to MgNiH 3 with H -0.55 . While in all three models the charge transfer is between Ni and H, the charge on Mg becomes slightly more ionic, changing from +1.61 to +1.66 and +1.73 along the series. At this point it becomes interesting to compare the relative charge results for cubic Mg 2 NiH 4 on the one hand and rutile-type MgH 2 on the other hand calculated with the same conditions. In MgH 2 the obtained charges of Q (Mg) = +1.64; Q (H) = -0.82 show a much more ionic character, close to -1 for H. By contrast, Mg 2 NiH 4 shows results closer to the presently studied model compounds: Q (Mg) = +1.63; Q (Ni) = -0.71 and Q (H) = -0.64. The charge carried by hydrogen is close to the one calculated for MgNiH 2 .

Fig. 2

 2 Fig.2shows in four panels the site projected density of states (PDOS) for tetragonal

Fig. 1 .

 1 Fig. 1. MgNiH x : the three model structures used in the calculations. a) MgNiH, b) MgNiH 2 , c) MgNiH 3 . The c tetragonal axis is along the paper sheet.

Fig. 2 .

 2 Fig. 2. (color) Site projected density of states: a) MgNi, b) MgNiH, c) MgNiH 2 and d) MgNiH 3 .

Fig. 3 .

 3 Fig. 3. (color) Chemical bonding for pair interactions (atom-to-atom) of a) MgNi, b) MgNiH, c) MgNiH 2 and d) MgNiH 3 .

Table 1 )

 1 leading to less overlap between the valence states of Mg and Ni. The consequence is the progressive destabilization of the metal substructure. Ni-H is the second bonding intensity interaction with much larger intensity than Mg-H which is vanishingly small because of the larger Mg-H distance and because Mg-H bonds are non-directional in as far as only s-like orbitals are involved.For equiatomic MgNi, the effects of hydrogen uptake on the electronic structure and the ionocovalent character of bonding have been approached using model structures. Hydrogen insertion leads to a destabilization of the metal substructure through its binding with nickel. In agreement with experiments, the most stabilized hydrogenated model structure is MgNiH 2 , close to the experimental value (MgNiH 1.6 ) in which hydrogen is found to carry a charge of H -

	Conclusion
	0.65	with substantial covalent character. This is concomitant with a medium
	absorption/desorption temperature of 200 °C whereas in ionic hydrides like MgH 2 the
	corresponding temperature is ~350 °C.

Table 1 .

 1 Calculation results for geometry optimizations and energy quantities for the studied compounds.

	Symmetry		MgNi	MgNiH	MgNiH 2	MgNiH 3
	Pm3 -m (cubic)	a , Å V, Å 3	3.05 (exp. 28.44) 28.32	-	-	3.32 37.01
		d Mg-Ni , Å	2.64	-	-	2.88
		d Ni-H1/H2 , Å				1.66
		d Mg-H1/H2 , Å	-			2.34
		Total energy	-7.45	-	-	-18.07
		per FU (eV)				INV (-
						14.65)
		Cohesive	-0.435	-	-	-0.223
		energy per
		atom (eV)
	P4/mmm	a, Å	2.987	3.09	3.26	-
	(tetragonal)	c/a	1.047	1.055	0.994
		V, Å 3	28.30 (exp.	31.08	34.4 (exp.
			28.44)		33.23)
		d Mg-Ni , Å	2.63	2.72	2.82	-
		d Ni-H1/H2 , Å		1.63	1.63
		d Mg-H1/H2 , Å	-	2.18	2.30
		Total energy	-7.46	-11.09	-14.65	-
		per FU (eV)		INV (-10.71)	INV (-12.31)
		Cohesive	-0.437	-0.226	-0.247	-
		energy per
		atom (eV)
	Charge trends	Q(Mg)	+1.61	+1.61	+1.66	+1.70
	(in electrons)	Q(Ni)	-1.61	-0.93	-0.41	-0.06
		Q(H1) 2	-	-	-0.63	-0.55
		Q(H2) 1	-	-0.67	-	-0.54

Acknowledgements

Parts of the calculations benefited from the facilities of MCIA-Université de Bordeaux.

We thank the CSR-USEK and URA-GREVE for financial support.


J F O N T A N E

P B R A N C H E R A N

Stochastic forcing of the Lamb-Oseen vortex

The aim of the present paper is to analyse the dynamics of the Lamb-Oseen vortex when continuously forced by a random excitation. Stochastic forcing is classically used to mimic external perturbations in realistic configurations, such as variations of atmospheric conditions, weak compressibility effects, wing-generated turbulence injected in aircraft wake, or free-stream turbulence in wind tunnel experiments. The linear response of the Lamb-Oseen vortex to stochastic forcing can be decomposed in relation to the azimuthal symmetry of the perturbation given by the azimuthal wavenumber m. In the axisymmetric case m = 0, we find that the response is characterised by the generation of vortex rings at the outer periphery of the vortex core. This result is consistent with recurrent observations of such dynamics in the study of vortex-turbulence interaction. When considering helical perturbations m = 1, the response at large axial wavelengths consists of a global translation of the vortex, a feature very similar to the phenomenon of vortex meandering (or wandering) observed experimentally, corresponding to an erratic displacement of the vortex core. At smaller wavelengths, we find that stochastic forcing can excite specific oscillating modes of the Lamb-Oseen vortex. More precisely, damped critical-layer modes can emerge via a resonance mechanism. For perturbations with higher azimuthal wavenumber m 2, we find no structure that clearly dominates the response of the vortex.

Introduction

From the large hurricanes developing in the atmosphere to the well-known Kelvin-Helmholtz billows in shear layers, vortices are ubiquitous in fluid flows. They are notably major actors of turbulence as they are involved in the energy cascade, entrainment and mixing. The understanding of their dynamics is then of considerable interest. In the context of aeronautics, the will to reduce the aircraft wake and the associated hazards to forthcoming planes has motivated the study of the stability of columnar vortices. Since the early works of [START_REF] Crow | Stability theory for a pair of trailing vortices[END_REF], [START_REF] Moore | The instability of a straigth vortex filament in a strain field[END_REF] and [START_REF] Tsai | The stability of short waves on a straight vortex filament in a weak externally imposed starin field[END_REF], many studies have shown the sensitivity of the wing-tip vortex pair to both long-and short-wave cooperative instabilities, see [START_REF] Leweke | Cooperative elliptic instability of a vortex pair[END_REF] or [START_REF] Billant | Three-dimensional stability of a vortex pair[END_REF] amongst others. Conversely, a single vortex is asymptotically stable but it supports various families of oscillating and damped modes amongst which are the so-called Kelvin waves [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF]. Although commonly used, the standard modal stability approach fails to fully predict the vortex linear dynamics. Indeed, transient growth can occur when specific perturbations are introduced in the flow. [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF] have calculated such disturbances for the Lamb-Oseen vortex and evidenced a core contamination mechanism combining Orr (1907a,b) and induction effects. The associated energy amplification can reach levels high enough to activate the nonlinearities and eventually lead to another equilibrium state or trigger a so-called "bypass" transition to turbulence. In the axisymmetric case, [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] (cited as PH06 in the following) recently encountered similar transient amplifications. [START_REF] Antkowiak | On vortex rings around vortices : an optimal mechanism[END_REF] (hereinafter referred to as AB07) completed the picture with the identification of a physical mechanism specific to vortices leading to the generation of vortex rings in the potential region around the vortex core.

The inability of modal analysis to prefigure such transient energy amplification was previously keyed out in wall-bounded shear flows, see [START_REF] Butler | Three-dimensional optimal perturbations in viscous shear flow[END_REF], 1993), Reddy & Henningson (1993) and Farrell & Ioannou (1993a). Thus both Couette and Poiseuille plane flows experience energy amplification even in their modal stability domain when aroused with the adequate perturbation. The physical mechanism involved consists in the emergence of strong streamwise velocity streaks emanating from streamwise vortices in the flow.

This general occurrence of transient growth in asymptotically stable flows is intimately related to the non-normality of the associated linear dynamical operator [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF][START_REF] Farrell | Variance maintained by stochasting forcing of nonnormal dynamical systems associated with linearly stable shear flows[END_REF]. Let us consider a flow represented by the following dynamical system for the state vector x

dx dt = Ax , (1.1)
where A is the operator corresponding to the Navier-Stokes equations linearised around a given basic state. The base flow is said to be asymptotically stable if all the eigenvalues of A have a negative real part. When the dynamical operator is non-normali.e. A H A = AA H with (.) H denoting the Hermitian transpose (Farrell & Ioannou 1993b[START_REF] Farrell | Variance maintained by stochasting forcing of nonnormal dynamical systems associated with linearly stable shear flows[END_REF])this approach does not address the issue of the flow energetics for finite times. Indeed, the stability analysis of a non-normal operator results in a set of modes that decay individually but that do not form an orthogonal basis. As a result, one can construct a perturbation on this basis with expansion coefficients that are large but with modes cancelling each other to give an initial energy of order one. Since each eigenmode evolves independently, the initial cancellation may not persist. The energy of the disturbance can thus increase substantially before decaying ultimately to zero. [START_REF] Schmid | Nonmodal stability theory[END_REF] gives an illustrative two-dimensional geometric example of this scenario in his recent review on nonmodal stability analysis. This mathematical property of non-normal operators is revealing of physical mechanisms that lead to transient energy amplification. The lack of orthogonality corresponds to the potential for energy extraction from the basic flow by a subspace of perturbations leading to transient growth despite the absence of modal (i.e. exponential) instability, a result already pointed out in the seminal work of Orr (1907a,b). This specificity of non-normal operators can be explained in a more formal way by considering the equation governing the instantaneous rate of energy change

dE dt = d dt x H x = x H A H + A x, (1.2)
where the energy of the system is defined by E = x H x in the H 2 norm. Energy growth occurs when the right-hand side of equation (1.2) is positive. This is mathematically equivalent to requiring that a portion of the numerical range of A, defined by N (A) = z ∈ C, z = x H Ax, x H x = 1 , lies in the right half-plane (Reddy et al. 1993). This can be turned into a condition on the spectrum of the energy operator K = A H + A, since the largest eigenvalue of 1 2 K, referred to as the numerical abscissa [START_REF] Schmid | Nonmodal stability theory[END_REF], is the supremum of Re [N (A)]. Transient growth of energy is likely to occur if K has at least one eigenvalue with a positive real part.

For flows experiencing such transient growth, it is of interest to find the "most dangerous initial condition" through a process of optimisation. It consists in identifying the initial condition that maximises its energy growth at a fixed time τ . Such a disturbance is referred to as an optimal perturbation. A second optimisation can be conducted with respect to time, leading to the optimal time τ opt that defines the most amplified optimal perturbation, often called the global optimal perturbation. Besides, in some cases, it is also relevant to look at the other disturbances experiencing transient energy amplification. They are termed sub-optimal perturbations. Among them, short-term optimals correspond to optimal perturbations for fixed times τ < τ opt [START_REF] Corbett | Optimal linear growth in swept bounday layers[END_REF].

Given the existence of such optimal disturbances in columnar vortices [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF][START_REF] Antkowiak | On vortex rings around vortices : an optimal mechanism[END_REF], the point is to know if they can naturally emerge from uncontrolled perturbations as diverse as atmospheric turbulence, background noise in wind tunnel experiments or turbulence generated by the aircraft wings. Indeed, while the potential for substantial transient growth of properly defined initial perturbations certainly exists, recurrent critics against optimal perturbation analyses concern the particular structure of these disturbances. These can be quite intricate and unlikely to occur spontaneously in real conditions since there is no apparent mechanism to excite such specific perturbations. This issue is theoretically addressed in the present paper to some extent. The general technique is to linearise the Navier-Stokes equations of small perturbations to a particular mean flow and then to augment these linear dynamics with stochastic forcing, which is uncorrelated in time (i.e., "white noise") and also possibly uncorrelated in space. This general maintained forcing can be interpreted as a crude way to mimic perturbations arising incessantly in real transitioning flows due to background turbulence or any kind of uncontrolled (in space and time) ambient fluctuations. But it should be kept in mind that the main objective of the present work is to understand the role played by the transiently growing disturbances found in the previous optimal perturbation analyses when the forcing lacks the bias of any specific forcing function. In that context, while this analysis may not necessarily apply directly to, say, real-life turbulent vortical flows where only certain types of disturbances may be introduced, it will nevertheless lend insight into the importance of the physical mechanisms of transient growth uncovered by the previous optimal perturbation analyses.

The associated dynamical equations can be thought of as a system where background noise is regarded as an "input" and the resulting random velocity field representing the response of the flow as the "output". The ratio of the output energy or variance to that of the input noise gives the energy amplification or gain of the system [START_REF] Schmid | Nonmodal stability theory[END_REF]. This analysis was successfully applied to wall-bounded shear flows (Farrell & Ioannou 1993b;[START_REF] Bamieh | Energy amplification in channel flows with stochastic excitation[END_REF] and to two-dimensional vortices with radial inflow derived from geophysical applications [START_REF] Nolan | The intensification of two-dimensional swirling flows by stochastic asymmetric forcing[END_REF]. In atmospheric sciences this approach is indeed a classical tool of investigation for the prediction of the statistics of meteorological flows such as the long-and short-term deviations of the hurricane tracks from that prescribed by the surrounding flow. For instance, [START_REF] Whitaker | A linear theory of extratropical synoptic eddy statistics[END_REF] used such a stochastic forcing analysis to recover successfully the observed variances of the winter Northern Hemisphere flow (in particular the location and structure of the storm tracks). Not only are these previous works in good agreement with the results issued from optimal perturbation analyses, pointing out the robustness of the growth mechanisms uncovered by these latter, but they also show to be complementary to them.

The paper is organised as follows. The stochastic forcing formulation is introduced for any general linear dynamical system in §2. Then its derivation for the Navier-Stokes equations is considered. The results are presented in §3, classically ordered with increasing azimuthal wavenumbers. The results are discussed and interpreted in §4. The paper ends with the conclusions and perspectives to this work in the last section §5.

Numerical formulation

Stochastically driven linear dynamical systems

The formalism employed in this section is classical in control theory and we only give a brief synthesis of the main steps for the sake of completeness. We consider an asymptotically stable linear dynamical system of the form (1.1) under the influence of an external forcing ξ. The system of governing equations can be conveniently written in the following state-space form classically used in the control literature [START_REF] Jovanovic | Componentwise energy amplification in channel flows[END_REF])

dx dt = Ax + Bξ, (2.1a) y = Cx, (2.1b)
where x is the state vector with initial condition x 0 and y is the output vector. The matrix B and C denote the input and output operators. The forcing considered is stochastic in nature and is assumed to be a temporal Gaussian white noise process with zero mean

ξ = 0, (2.2a) ξ(t)ξ H (t ′ ) = Rδ(t -t ′ ), (2.2b)
where is the ensemble averaging operator and R the spatial covariance matrix. For purpose of representativeness of a specific configuration, the matrix R can be implemented with respect to experimental data. In the absence of such information it can be set equal to the identity matrix, i.e. , R ij = δ ij , leading to a spatio-temporal Gaussian white noise. Some refinements may be included in the forcing term in order to represent more specific perturbation fields. For instance, one can let the forcing amplitude peak near the walls to take into account their influence in wind tunnel experiments as done by [START_REF] Jovanovic | Componentwise energy amplification in channel flows[END_REF] for wall-bounded shear flows. We limit here the statistical properties of the forcing terms to the above-mentioned so as to mimic the most generic free stream disturbances occurring in real conditions without favouring any particular region of the flow. As claimed by Farrell & Ioannou (1993b), this forcing is not intended to reproduce the full complexity of turbulence observed in experiments. Its aim is to retain the essential physics underlying the maintained variance from any external continuous perturbation field. In that sense, the analysis will give an insight into the receptivity of the flow without introducing any a priori bias through the characteristics of the forcing in the physical or spectral space.

As a response to this forcing, the emerging state y is a stochastic process with secondorder statistics given by the covariance matrix

y(t)y H (t) = C t 0 e A(t-s) BB H e A H (t-s) ds C H = CQ(t)C H , (2.3)
where Q(t) is referred to as the controllability Gramian. If the dynamical operator A does not vary in time, the system will reach a statistical steady state where the matrix

Q ∞ = lim t→∞ Q(t) is solution of a Lyapunov equation AQ ∞ + Q ∞ A H = -BB H .
(2.4)

The mean energy corresponds to the variance of the output stochastic process and can be extracted from the covariance, i.e. E(t) = y H (t)y(t) = trace B H Q(t)B . As classically done in the control literature [START_REF] Zhou | Robust and Optimal Control[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF], this quantity can be interpreted as the H 2 -norm of the transfer function H(ω) = C(iωI -A) -1 B associated to the linear system (2.1). To obtain the energy amplification denoted as G ∞ = E ∞ / E inp , one has to determine the input energy introduced by the stochastic forcing. It comes from the variance equation which reads

d dt E(t) = x H (A H C H C + C H CA)x + x H C H CBξ + ξ H B H C H Cx .
(2.5) Using (2.1b) and (2.2b), the input energy appearing on the left-hand side of (2.5) becomes

E inp = x H C H CBξ + ξ H B H C H Cx = 2 trace B H (C H C) 2 B . (2.6)
The extraction of coherent structures from the random flow field resulting from the forcing can be obtained through the computation of the eigenmodes of the controllability Gramian. The eigenvalue decomposition of Q ∞ , known as the Karhunen-Loève (KL) or proper orthogonal decomposition (POD), will provide with the flow patterns that participate to the response, ordered according to their contribution to the variance of the statistically steady state [START_REF] Schmid | Nonmodal stability theory[END_REF]) p) .

Q ∞ z (p) = γ (p) z (
(2.7)

These distributions will be referred to as output structures hereinafter. It is also of interest to know which coherent structures from the input noise participate the most to the excitation of the system. This issue can be addressed by considering the adjoint dynamical system † (A H , B H , C H ) forced with a similar Gaussian white noise. The second-order statistics of the adjoint output stochastic process are contained in the covariance matrix B H P (t)B where P (t) is the observability Gramian defined as

P (t) = t 0
e A H (t-s) C H Ce A(t-s) ds .

(2.8)

Its long-time value is also solution of a Lyapunov equation

A H P ∞ + P ∞ A = -C H C.
(2.9)

Finally, the eigenvalue decomposition of P ∞ , coined the back Karhunen-Loève decomposition by Farrell & Ioannou (1993b), will provide with the coherent forcing structures (referred to as input structures in the sequel) hierarchised according to their contribution to the system excitation

P ∞ ζ (p) = β (p) ζ (p) .
(2.10)

Stochastic forcing applied to Navier-Stokes equations

The base flow considered in this paper is the Lamb-Oseen vortex. Its non-dimensional azimuthal velocity field is

V (r) = 1 r 1 -e -r 2 , (2.11)
where the characteristics scales used are the vortex dispersion radius r 0 and the angular velocity at the axis Ω 0 . In the following, Ω(r) = V (r)/r = (1 -e -r 2 )/r 2 and Z(r) = (1/r)∂ r (rV (r)) = 2e -r 2 will represent respectively the angular velocity and the axial vorticity of the columnar vortex. We now consider infinitesimal disturbances expanded in a modal form

[û r , ûθ , ûz , p ](r, θ, z, t) = [u, v, w, p ](r, t)e imθ+ikz + c.c. , (2.12)
where c.c. stands for the complex conjugate, m is the azimuthal wavenumber and k the axial wavenumber. The linearisation of the incompressible Navier-Stokes equations gives the following set of equations

1 r ∂ r (ru) + im r v + ikw = 0 , (2.13a) ∂ t u + imΩu -2Ωv = -∂ r p + 1 Re ∆ m,k - 1 r 2 u - 2im r 2 v , (2.13b) ∂ t v + imΩv + Zu = - im r p + 1 Re ∆ m,k - 1 r 2 v + 2im r 2 u , (2.13c) ∂ t w + imΩw = -ikp + 1 Re ∆ m,k w , (2.13d )
where the Reynolds number is Re = Ω 0 r 2 0 /ν = Γ/(2πν) and

∆ m,k ≡ ∂ rr + 1 r ∂ r - m 2 r 2 -k 2 . (2.14)
The set of variables is reduced to the three velocity components by eliminating the pressure with the Poisson equation. The state variables correspond to the primitive ones, i.e. x = [u, v, w], and both the input and output operators B and C of §2.1 simply reduce to the identity operator. A Chebyshev spectral collocation method with an algebraic mapping identical to that described by [START_REF] Fabre | Viscous instabilities in trailing vortices at large swirl numbers[END_REF] and [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF] is used for the spatial discretisation of the problem. The energy of the perturbation, which defines a weighted inner product, is defined classically by

E = (x, x) E = 1 2 ∞ 0 (u * u + v * v + w * w) rdr, (2.15)
where * denotes the complex conjugate. The use of the energy-based inner product requires a coordinate transform to convert the statistics measures into the more standard Hermitian H 2 norm [START_REF] Hoepffner | Stability and control of shear flows subject to stochastic excitations[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF]. From the Choleski factorisation of the symmetric definite positive weight operator W = M H M , the relation between the two norms is ||x|| E = ||M x|| 2 and ||A|| E = M AM -1 2 . The calculations are carried out with Matlab using the DMSuite package developed by [START_REF] Weideman | A matlab differentiation matrix suite[END_REF]. The convergence of the numerical procedure depends on both the truncation level N and the radial domain extent r max . We tested the results sensitivity to both parameters and converged calculations were obtained for N ∈ [150, 300] and r max ∈ [5000, 20000], according cases.

Results

The axisymmetric case (m = 0)

The presentation of the results starts with the axisymmetric case. Figure 1 shows the energy amplification as a function of the axial wavenumber for Reynolds numbers varying from 500 to 10000. As in optimal perturbation analyses (PH06, AB07), high levels of amplification are observed. The plot also confirms that the largest amplifications occur for small axial wavenumbers. The amplification levels reach higher values than those found by optimal perturbation analyses, which can be easily understood since the system is excited by a constant input and not solely by an initial condition [START_REF] Bamieh | Energy amplification in channel flows with stochastic excitation[END_REF]. This difference in the level of amplification between these two approaches is systematic.

For the k = 1 case, we consider the hierarchy of both the input and output structures according to their contribution to the sustained variance of the steady state. The spectra resulting from the eigenvalue problems (2.7) and (2.10) are plotted in figure 2. The energy amplification essentially results from the contribution of the few first structures as the eigenvalues rapidly decrease to negligible values. The first input structure is responsible for 62% of the vortex excitation and dominates the forcing field. In the same manner, the first ouput structure prevails in the vortex response by contributing to 61% of the variance sustained by the flow. Figure 3 gives the spatial distribution of these two dominant structures. They are very similar to the optimal perturbation at initial and optimal times found in the optimal perturbation analyses of PH06 and AB07. The forcing structure consists of azimuthal velocity streaks located in the quasi-potential region in the outer periphery of the vortex. The physical mechanism leading to the azimuthal vortex rings of the figure 3(b) has been previously explained by AB07. It has been called "anti-lift-up" in reference to the so-called "lift-up" mechanism occuring in planar shear flows, though it is radically different in nature. Briefly, the azimuthal velocity streaks induce a local Coriolis force field yielding to a radial displacement of the fluid particules. Its potential part is balanced by the pressure gradient to ensure the flow incompressibility. Its rotational part feeds the tori of azimuthal vorticity that constitute the vortex response displayed in figure 3(b). AB07 established this original mechanism with the assumption of zero base vorticity. This is a good approximation since the velocity streaks are localised in the quasi-potential region of the base-flow. The essential trends are correctly captured by this scenario. Nevertheless, the residual base vorticity in the outer periphery of the Lamb-Oseen vortex is at the origin of a small deviation from this idealized view when approaching the vortex core. Indeed, one can identify weak secondary rolls in figure 3(b) located closer to the vortex core. [START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF] showed that they are the signature of waves generated in the region of non-zero epicyclic frequency of the base-flow κ defined by κ(r) 2 = 2Ω(r)Z(r) (Rayleigh discriminant). Mathematically, this comes from the coupling term associated with the base vorticity Z(r) in the equation for the azimuthal velocity, see equation (2.13c).

When varying the axial wavenumber, the radial location of both the azimuthal velocity streaks and the resulting azimuthal vortex rings increases with increasing wavelengths (data not shown). This is in agreement with the results of PH06 who related the increasing radial position of the optimal perturbation (and consequently the associated optimal time) with increasing k to the radial distribution of the vorticity-to-strain ratio in the base flow. The anti-lift-up mechanism is observed to dominate the flow dynamics even more as the axial wavenumber decreases. Indeed, contribution of the dominant structures to the variance of the response is found to shrink at large k. For k = 0.1, the output (resp. input) structure accounts for 78% (resp. 79%) of the energy amplification, whereas for k = 2.5 the output (resp. input) structure only participates for 32% (resp. 32%) of the gain.

The next point to be discussed is the influence of the Reynolds number. From figure 1, it is obvious that larger Reynolds numbers result in larger energy amplifications by the flow. This is not surprising because of the energy balance between the extraction of energy from the background flow by the stochastic forcing, on the one hand, and the viscous dissipation on the other hand. When looking more precisely at the Reynolds number dependence for the sustained variance, one can obtain a scaling law from a sufficient collection of calculations, see figure 4(a) for k = 1. For the range of axial wavenumbers explored here, the scaling was always found to be of the form G ∞ = O(Re α ) with α decreasing from 3 for k = 0.1 to 2.5 for k = 3. [START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF] showed that a O(Re 2 ) scaling for the transient energy growth is an upper bound. In the present study, the system is being continuously excited and the energy accumulates before it dissipates. The characteristic diffusion time is O(Re). Hence, the upper bound for the energy growth is O(Re 3 ) here, which is consistent with the results obtained. This scaling has been derived with the hypothesis of zero base vorticity. The deviation from the predicted scaling law is explained by the existence of the radially propagating waves mentioned previously, which are responsible of an outward radiation of energy [START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF]. This point also enlightens on the variation of the exponent α with k in the scaling law found in this study. As k is decreased, the input structures are located further away from the core where the zero base vorticity assumption is more accurate (PH06, AB07). Hence for small wavenumbers, there is no wave generation and the scaling law in O(Re 3 ) becomes exact.

In addition, variation of the viscous diffusion does not affect the spatial structure of forcing and response. The velocity and vorticity distributions displayed in figure 3 remain unchanged when the Reynolds number is varied. The only difference resides in the peaks of both azimuthal velocity and vorticity that rise to higher values when the Reynolds number is increased. Finally, we also obtained an empirical scaling law for the energy amplification when k goes to zero. We found numerically that the gain roughly scales as O(Re 3 /k 2 ) as can be seen in figure 4(b). 3.2. The helical case (m = 1)

Figure 5 displays the energy gain for the m = 1 helical waves as a function of the axial wavenumber for different Reynolds numbers. Compared to the optimal perturbation analyses of [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF] and PH06, similar global trends are retrieved and large amplification levels are found. The large increase of the gain for small k is observed as well as the peak around k = 1.5 and an emerging one at k = 2.5 although less pronounced than in the optimal perturbation analysis of [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF]. The latter point suggests that the dominant structure does not emerge conspicuously in the vortex response for this range of wavenumber. This is confirmed quantitatively since the first output structure only contributes to 32% of the sustained variance for k = 1.35 and Re = 1000. This is the consequence of the coexistence of several perturbations participating in the energy amplification when the system is stochastically forced. For k = 1.35, the second emerging structure contributes for 11% to the gain. Hence, the vortex response to the forcing will be dominated by the first output structure. Its axial vorticity is plotted in figure 6 as well as the axial vorticity of the input structure that essentially excites it. This forcing structure is composed of a pair of left-handed spiraling vorticity sheets similar to the initial shape of the optimal perturbation found by [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF]. These two folded vorticity layers located in the quasi-potential region of the flow are of alternate sign. Their respective velocity induction on the vortex core initially cancel each other. As time evolves, they progressively uncoil via an Orr mechanism induced by the base flow differential rotation. Through this process, the reorganisation of the vorticity sheets promotes an increasing velocity induction in the vortex core that leads to the emergence of the dominant output structure of figure 6(b). This response corresponds to the first critical layer mode of [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF] as can be seen from the comparison of the velocity radial profiles in figure 7(a). This phenomenon can be interpreted as a transient resonance mechanism: a Kelvin wave of the Lamb-Oseen vortex is excited by the perturbation field induced by the uncoiling of the initial vorticity spirals. This conjecture was first confirmed by Antkowiak (2005) who made a successful comparison with a forced harmonic oscillator elementary model. Choosing an adequate forcing term so as to correctly reproduce the influence of the spiraling vorticity arms, the frequency of the forcing term is found to fit very well those of the selected waves. PH06 performed a similar analysis on the top-hat (Rankine) vortex. Their results supported the hypothesis of a resonance-driven excitation of vortex waves. From the simple harmonic oscillator model, one can obtain an expression for the energy gain at large time when the forcing term has a constant amplitude: G ∞ = 1/(s 2 + (ω f -ω) 2 ) where ω f is the pulsation of the forcing term, ω the pulsation of an eigenmode of the system and s its damping rate. From this first estimation of the energy growth, one can see that the selection process is based on two aspects. The first criterion consists in selecting one of the least damped waves -minimising the s term -since it allows a maximal energy amplification. For the present set of parameters, the first critical layer mode -referred to as the first L1 mode according to the nomenclature of [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF] -has the minimum damping rate, see their figure 6. As a consequence it is found to be the dominant output structure. Secondly, the disturbance may be in phase with the wave in order to impose a continuous excitation -minimizing the (ω f -ω) term. This implies that the rotation rate of the spirals be close to the pulsation of the mode. In the quasi-potential region of the flow, the rotation rate of the vorticity sheets is Ω(r) ≈ 1/r 2 . Hence, the radial position of the spiraling vorticity sheets is imposed by the mode frequency, i.e. r ≈ 1/

-4 -2 0 2 4 -4 -2 0 2 4 (a) -3 -2 -1 0 1 3 -3 -2 -1 0 1 2 3 (b) -4 -2 0 2 4 -4 -2 0 2 4 (c) -3 -2 -1 0 1 3 -3 -2 -1 0 1 2 3 (d)
√ ω where ω is the wave pulsation, a criterion previously established by PH06 with the top-hat vortex model (see their figure 21c). For the k = 1.35 case, the first L1 wave frequency is 0.121, which gives a mean radial position of r = 2.9. This is consistent with the location of the spiral arms in figure 6(a). Keeping in mind these two criteria for the resonance phenomenon and considering the stability results of [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF], it is expected that the mode selection depends on the axial wavenumber. Indeed for larger wavenumbers, we found the second L1 mode as the first output structure, see figure 6(d) for k = 2.5. When k is decreased, the emerging structure displayed in figure 6(f ) is also of a different nature and corresponds to the displacement wave of the Lamb-Oseen vortex (referred to as the D wave in the following). According to the criteria previously mentioned, the emergence of the D mode is not linked to an exact resonance. This wave being countergrade, i.e. its pulsation is negative, no disturbance can have a similar rotation rate. In this case, the D wave emerges only because its damping rate is very small compared to those of the waves of the L family. It is noteworthy that the spiraling input structure is located far from the vortex core, at a mean radial position r ≈ 15. This represents a slow time-rotating disturbance with an associated rotation rate Ω(r) ≈ 0.0044. This point is linked to the fact that the displacement wave emergence corresponds to an exact resonance in the two dimensional limit for a steady forcing. For k = 0, the D wave is stationary and the exciting perturbation can satisfy the in-phase condition if the spiral arms of vorticity are pushed outwards to an infinite radial position. Obviously, this could not be verified numerically, but we found that the radial location of the input structure increased when reducing the axial wavenumber from k = 0.5 to k = 0.1. This exact resonance comes with an infinite energy growth as revealed by the large increase of the gain when k goes to zero in figure 5.

The emergence of the displacement wave under continuous external perturbations is a good candidate for explaining the vortex meandering (or wandering) phenomenon. The nature of the wave, namely a long wavelength bending wave, is in agreement with the characteristics of the meandering observed in wind tunnel experiments [START_REF] Baker | Laser anemometer measurements of trailing vortices in water[END_REF][START_REF] Devenport | The structure and development of a wing-tip vortex[END_REF]. The theoretical approach used here, which takes into account the external disturbances in the form of a continuous random forcing, is consistent with real experimental conditions where background noise is constantly exciting the wing-tip vortex. This point has to be confirmed experimentally by controlling, or at least, quantifying, the perturbation field and by correlating the data with measurements of the vortex response. Here we only give a physical mechanism that could explain the occurrence of this phenomenon. But we are still far from a comprehensive and fully predictive description that could help experimentalists to filter out this erratic core displacement. To our knowledge, the best way to proceed is proposed by [START_REF] Devenport | The structure and development of a wing-tip vortex[END_REF], who used a Gaussian model for the fluctuations of the vortex core position.

Inspection of the effect of viscosity has been explored numerically and a general power law has been obtained for the energy amplification, see figure 7(b) for k = 1. For the range of axial wavenumbers considered here, we found the scaling G ∞ = O(Re 1.8 ). This exponent is considered as valid for any k since very small variations were observed, i.e 1.8 ± 0.02. [START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF] showed that the transient energy growth is O(Re) for the resonance phenomenon while it is O(Re 2/3 ) when only the Orr mechanism is active. Following the same reasoning as in the axisymmetric case, the theoretical upper bound for energy amplification under stochastic forcing is O(Re 2 ). The results obtained here are in good agreement with this scaling law. A decrease in viscous diffusion also affects the structures. While the vorticity distribution of the Kelvin waves remains practically unchanged, the spiral arms of the dominant disturbance get thinner as the Reynolds number increases, a point already noticed by PH06. The peaks of vorticity also reach larger values, which leads to a stronger velocity induction in the core when the spirals unfold and a higher level of energy amplification.

The double-helix case (m = 2) and higher azimuthal wavenumbers

The energy gain of the vortex in response to double-helix perturbations is plotted versus the axial wavenumber in figure 8 for different Reynolds numbers. The levels of amplification are lower than in the axisymmetric case but with similar tendencies. For a given Reynolds number, the larger the wavelength the larger the amplification. These results are noticeably different from those of the optimal perturbation analysis (PH06). Not only are the levels much higher but the dependence on k is completely different. The amplification maximum occurs for small wavenumbers and the curve does not present the flat-shaped aspect displayed in figure 19(b) of PH06. The difference in the results can be explained by two complementary arguments. Firstly, the development of the optimal perturbation is completed within few rotation periods of the vortex [START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF]. In this case, the transient energy growth of the optimal perturbation may occurs on too short a timescale, compared to the statistical mean time between two consecutive excitations of the mode by the external forcing, for being efficiently amplified. After an excitation of the vortex with the optimal disturbance extracted from the random noise, the process of energy amplification and decay is likely to be over before another excitation happens. Hence, the energy gain resulting from successive excitations is statistically unlikely to cumulate. As a result, the optimal perturbation will be less amplified under a random forcing than a disturbance associated with a less efficient transient mechanism but occurring on a longer timescale. This point is related to the second argument. Indeed, strong sub-optimal perturbations competing with the optimal one may exist. They can eventually dominate the response of the vortex when submitted to stochastic maintained forcing if they develop on a larger timescale than the optimal one. Evidence to this conjecture can be obtained by looking at both the input and output structures.

The structures that dominate the input forcing consist of the same kind of spiralling structures as the ones identified in the m = 1 helical case (fig. 9). The resonance mechanism leading to the selection of a Kelvin wave is active and the first output structure for k = 2.5 displayed in figure 9(b) corresponds to the flattening wave of [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF]. This mode being by far the least damped, it is selected on frequency arguments according to the criteria mentioned in the previous section. The mean radial position of the entangled vorticity spirals is imposed by the pulsation of the mode. However, while active, the resonance no longer dominates the flow dynamics since the first input structure account for only 8% of the sustained variance in this case. The sub-optimals participate equally to the excitation of the vortex. For smaller wavenumbers, transient resonance is not even the dominant mechanism leading to energy amplification. Figure 9(d) shows the first output structure for k = 1. The vortex response is only composed of the four initial vorticity sheets that have been uncoiled. In this case, the Orr mechanism gives the largest energy growth and accounts for 10% of the total amplification. The spiralling vorticity arms of the input structure are located far from the vortex core so that the timescale of the Orr mechanism is long enough to be efficiently amplified by the stochastic forcing. Its timescale is given by the duration of the uncoiling of the vorticity arms which is roughly evaluated by

τ Orr ≈ - 2π e∂ r Ω(r h ) , (3.1)
where e and r h are respectively the mean width and the radial position of the spiral arm. For the perturbation plotted in figure 9(c), τ Orr is found to be about 50 rotation periods of the vortex, which is one order of magnitude larger than the optimal time. As previously argued, this energy growth mechanism is activated on a longer timescale than the transient resonance, and it is more efficiently amplified by the maintained forcing. Inspection of the less dominant structures reveals that the flattening wave is the fourth preferred response of the flow accounting only for 4% of the variance. For the m = 2 case, there is no structure that clearly dominates the flow response, by contrast to what has been found for the m = 0 and m = 1 cases. Whatever the axial wavenumber, the contribution of the most dominant output structure never exceeds 15%. Finally, our survey ends with the description of the vortex response for higher azimuthal wavenumbers. We explored the flow behaviour under a continuous stochastic forcing for m as large as 10 and always found energy amplification. The observed gain levels decrease with increasing m and are far below those of the three previous cases. The amplification curves present the same shape as for the m = 2 case with a maximum when k tends to zero. The dynamics leading to energy growth always consists of the Orr mechanism. The dominant input structure lies in the spiralling vorticity sheets that uncoil due to the differential rotation. No transient resonance mechanism has been found as could have been expected considering the large damping rate of the Kelvin waves [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF]). This phenomenon is valid for all m > 2. While energy growth is observed, these cases are of limited interest as they are exceedingly unlikely to prevail over the responses corresponding to the first three azimuthal wavenumbers.

Discussion

The results obtained in the present study show a noticeable similarity between the amplification curves for the different azimuthal wavenumbers considered here (m ∈ [0, 10]): the energy gain is maximum for k = 0 and decreases as k increases. This theoretical lack of intrinsic axial wavelength selection raises questions regarding the use of the present results in order to predict the response of the vortex in real-life conditions.

At that stage, we wish to underline that the interpretation of the amplification curves must be undertaken carefully. If one wants to predict the selection of a particular wavelength, a direct and hasty extrapolation of the results presented in this paper can prove completely erroneous when compared to experiments without any caution. Thus, considering that the most amplified (k, m) modes correspond to m = 0 and k tending to zero, it could be tempting to predict the systematic and dominant occurrence of very large axisymmetric vortex rings at the outer periphery of the vortex, as presented in section 3.1. But this conclusion, if applied to realistic experiments or observations, must be balanced by taking into account the departure of the experimental conditions from the hypothesis that define the framework of the present analysis. It concerns the unbounded character of the flow, the spectral content of the external forcing and the influence of the nonlinearities. These points are discussed further in this section.

Bounded vs unbounded flows

First, the largest amplifications are obtained in the limit of infinite wavelengths (k → 0) for every azimuthal wavenumber m, suggesting that the response of the vortex should be systematically quasi-two-dimensional, at least for m > 0, and involve very large structures. But it must be kept in mind that the input structures triggering the vortex response (azimuthal streaks for m = 0 and axial vorticity spirals for m > 0) are localised radially further away from the vortex axis when k decreases, eventually extending to infinity for k → 0. This singular behaviour questions the validity of the unbounded flow hypothesis that is implicitly made in the present formulation. It is expected that the predictions concerning the long-wave response of the vortex put forward in the previous sections will be distorted by the presence of physical boundaries at finite distance such as the side walls of the wind tunnel, or by taking into account flow features ignored in the present model such as other vortices for instance. The present study found no intrinsic wavelength selection by the vortex and the process of wavelength selection is extrinsic and therefore case-dependent.

Nevertheless, the present results show that the physical mechanisms of growth involved in the response of the vortex systematically favour the structures with the largest axial wavelength admissible as long as their radial extent does not exceed the limit of representativeness of the vortex flow model used here. We then expect the present results to be relevant above a critical axial wavenumber k c ∼ 2π/r c corresponding to a characteristic radial extent r c above which the vortex flow model used here significantly departs from the real-life conditions.

This discussion can also be supplemented, to a lesser extent, with the validity of the infinite time limit required in the present formulation to reach theoretically a statistical steady state. The input structures that are favoured in the k → 0 limit correspond to infinitely slowly growing disturbances [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF][START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF]. These perturbations thus may not have enough time to grow significantly in finite-time experiments. This is a second source of distortion for the large-wavelength results presented here. If the flow develops on a short or medium timescale, one can expect the selection of finite-time optimal perturbations which can be radically different from the structures predicted here for an infinite time horizon.

Initial turbulence vs. continuous white noise

The linear evolution of perturbations of a vortex flow is formally given by the general solution of the governing equation (2.1a) which has two components: the homogeneous solution which describes the evolution of initial conditions and the particular solution which represents the long-time response of the flow to a continuous external forcing. As mentioned by [START_REF] Schmid | Nonmodal stability theory[END_REF], both parts are complementary and fully describe the general dynamics of small perturbations.

The studies of perturbed vortices generally focus on the former point, the evolution of initial conditions, such as temporal modal stability analyses (see [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF] for the Lamb-Oseen vortex), optimal perturbation analyses [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF][START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] or theoretical and numerical investigations of the response of vortices to initially injected turbulence [START_REF] Melander | Coupling beetwen a coherent structure and fine-scale turbulence[END_REF][START_REF] Risso | Direct numerical simulations of wake vortices in intense homogeneous turbulence[END_REF][START_REF] Miyazaki | Linear and non-linear interactions between a columnar vortex and external turbulence[END_REF][START_REF] Takahashi | The influence of turbulence on a columnar vortex[END_REF][START_REF] Marshall | External turbulence interaction with a columnar vortex[END_REF]. More particularly in the latter case, the objective is to understand how the vortex immersed in an initial turbulent field responds to this perturbation and how in return the initial turbulence is affected by the presence of the vortex and the associated shear and rotation which are known to drastically alter the statistics of turbulence on a short time scale.

By contrast with these studies, the approach adopted in the present paper is quite different, as it focuses on the long-time response of the vortex when submitted to a continuous external forcing. This corresponds formally to the study of the particular solution of equation (2.1a) as stated by [START_REF] Schmid | Nonmodal stability theory[END_REF]. Physically this particular solution provides a model for receptivity process, where the external forcing may represent either free-stream turbulence, wall roughness or other non smooth geometries, body forces or even neglected terms such as nonlinearities. This forcing can also be linked to the deviation from the model base flow, like the faraway presence of other vortices or the boundary-layer perturbations generated at the side walls of the wind tunnel. In that context of receptivity analysis, the continuous forcing used in the study is chosen as generic and unbiased as possible in the form of white noise. Being equally distributed in space allows to excite all possible regions of the flow and not to favour any specific wavelength or frequency. Thus, if the wavelength selection is case-dependent (on the dominant wavelengths or frequencies of the external forcing in real-life conditions for instance †, see preceding section), the quantification of the amplification is nevertheless a measure of the intrinsic transfer function of the Lamb-Oseen vortex, and reveals the intrinsic mechanisms that are favoured by the flow. Here it is shown that the largest response is observed in forcing scenarios that convert perturbations in the form of azimuthal velocity streaks into intense azimuthal vortex rings (for the axisymmetric part of the flow), and perturbations in the form of spiraling axial vorticity sheets outside the vortex into large bending (m = 1) and deformation waves (m > 1) within the vortex. It is noteworthy that these mechanisms are similar to the ones uncovered by nonmodal studies of the initial-condition problem such as optimal perturbation analyses (PH06, AB07) and are consistent with the numerical simulations of the interaction of a vortex with an initial turbulent field. This suggests that such mechanisms are robust and fundamental to the dynamics of perturbed vortices. It is expected that they are potentially active whatever the details of the vortex flow and the perturbations considered.

Validity of the linear approach

This work has been conducted within the linear approximation and it is important to figure out how the results obtained here are adjusted by the nonlinearities. The physical mechanisms of transient energy amplification considered here and in the studies of PH06 and AB07 are linear. As argued by [START_REF] Miyazaki | Linear and non-linear interactions between a columnar vortex and external turbulence[END_REF], linear analysis correctly describes these processes or at least their initial development, while the nonlinearities impact the further evolution of individual structures. This point is currently under investigation via direct numerical simulations of the nonlinear evolution of the optimal perturbations. Preliminary results for the m = 0 case show that the azimuthal vorticity rings of the output structure are self-advected away from the vortex axis, carrying along the streaks of high azimuthal velocity of the forcing structure. Due to the conservation of angular momentum, the intensity of the streaks decreases as the structures move radially outwards. Thus, the axisymmetric mechanism of energy growth is progressively damped by the nonlinearities. In the case of a maintained forcing, this nonlinear damping effect is expected to be of minor importance as the streaks of azimuthal velocity would be statistically regenerated continuously.

Concerning the external forcing from a more general perturbation field, the nonlinear development of the flow has to be considered for all wavenumbers simultaneously. Numerical simulations including a continuous forcing would be of great interest in order to know how the amplification factors obtained here are affected by the nonlinear effects (i.e. , saturation). In the case of a transient forcing induced by an initially homogeneous isotropic turbulence, [START_REF] Takahashi | The influence of turbulence on a columnar vortex[END_REF] performed such analysis. The departure from the linear regime happens when the bending wave growing in the core reaches a finite amplitude. They measured quantitatively the impact on nonlinearities through the time dependence of the axisymmetric axial correlation function. Its temporal growth was observed to be proportional to t contrary to linear RDT analyses where a quadratic time dependence was found [START_REF] Miyazaki | Linear and non-linear interactions between a columnar vortex and external turbulence[END_REF]. Thus, the nonlinear terms tend to limit the energy growth and a similar saturation is expected for any generic random forcing.

Conclusions

The goal of the present study is to analyse the dynamics of the Lamb-Oseen vortex when continuously forced by a random excitation. Considering the existence of transient growth in vortices when submitted to specific perturbations, this work aims at determining if the optimal perturbations found by [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] could naturally emerge from background noise. For this purpose, the linear Navier-Stokes equations are continuously forced with Gaussian white noise so as to mimic any perturbations occurring in real transitioning flows. We then look at the large-time statistical response of the vortex and at the coherent structures participating to its excitation. Not only this method allows the energy gain of the system to be quantified but also to find and order the coherent input (resp. output) structures according to their contribution to the excitation (resp. variance or response) of the flow.

For all azimuthal wavenumbers investigated, energy amplification is always observed but when m > 2, the levels reached are too small compared to those obtained from smaller values of m to be significant. Compared to the optimal perturbation analysis, the levels of amplification obtained here are always higher since it corresponds to energy growth resulting from a continuous noise input and not from a single initial condition. This difference in the value of the gain is particularly marked when the optimal perturbation evolves on a short timescale. In such cases, the sub-optimal perturbations have to be considered as they participate equally in the vortex response. This point has to be recognised as one advantage of the stochastic forcing approach. Indeed, the optimal perturbation analysis classically focuses on the most amplified disturbance and can sometimes miss some transitional scenarios. This is what we observe for the m = 2 case where no prevailing mechanism has been found to dominate the vortex dynamics.

Focusing on the physical mechanisms leading to energy growth, the scenarios leading to the vortex excitation identified by [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] are recovered. This study confirms that the optimal perturbations can naturally be activated by the background noise present in uncontrolled conditions.

In the axisymmetric case, the mechanism called "anti-lift-up" by AB07 consists in the emergence at the vortex periphery of strong tori of azimuthal vorticity fed by the azimuthal velocity streaks of the forcing structure. This scenario gives a theoretical counterpart to the observation of the recurrent development of vortex rings at the periphery of vortices when submerged in an ambient turbulence [START_REF] Melander | Coupling beetwen a coherent structure and fine-scale turbulence[END_REF][START_REF] Risso | Direct numerical simulations of wake vortices in intense homogeneous turbulence[END_REF][START_REF] Takahashi | The influence of turbulence on a columnar vortex[END_REF].

For helical perturbations, the excitation of the vortex is explained by a transient resonance phenomenon. The dominant forcing structure is composed of two left-handed entangled vorticity sheets located in the quasi-potential region of the flow. As time evolves, they progressively uncoil under the base flow differential rotation and trigger the appearance of a Kelvin wave through a process combining induction effects and the Orr mechanism. The emerging Kelvin mode depends on the axial wavenumber as the selection is based both on a minimisation of the wave damping rate and a concordance between the pulsations of both the mode and the vorticity spirals. At large wavelengths where the gain is maximal, the displacement mode of [START_REF] Fabre | Kelvin waves and the singular modes of the Lamb-Oseen vortex[END_REF] is preferentially excited. The emergence of the displacement wave under the influence of a maintained background noise is thus an interesting candidate for the vortex meandering observed in the experiments [START_REF] Baker | Laser anemometer measurements of trailing vortices in water[END_REF][START_REF] Devenport | The structure and development of a wing-tip vortex[END_REF]. This mechanism of resonance originating from a localised spiraling vorticity disturbance also occurs for m = 2 perturbations with the emergence of the flattening wave. But it is in competition with other disturbances experiencing transient growth. The sole disentanglement of the vorticity spirals via the base flow differential rotation leads to energy amplifications as large as those deriving from the flattening mode resonance. The existence of sub-optimal perturbations living on a longer timescale than the optimal one is put forward to explain the discrepancy between the results obtained here and those of PH06. This point has to be confirmed. It will be investigated in the near future with the computation of the sub-optimal perturbations of the Lamb-Oseen vortex for m = 2. For larger wavenumbers, the Orr mechanism becomes the dominant process of energy amplification from a continuous random excitation of the vortex.

Finally, the approach conducted here fulfills the interrogations brought by the results of the optimal perturbation studies of [START_REF] Antkowiak | Transient energy growth for the Lamb-Oseen vortex[END_REF][START_REF] Schmid | Nonmodal stability theory[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF]. This work demonstrates both analyses to be complementary. Such tools may be applied to asymptotically stable flows whose dynamical operator is non-normal in order to find out an eventual transient mechanism. For example, [START_REF] Joly | The Rayleigh-Taylor instability of twodimensional high-density vortices[END_REF] recently showed that low-density vortices were insensitive to the Rayleigh-Taylor instability. Since the corresponding dynamical operator does not commute with its Hermitian transpose, nonmodal stability analysis (both initial value and stochastic forcing formulations) should be conducted.

Figure 1 .

 1 Figure 1. Energy amplification reached for the statistical steady-state in the axisymmetric case (m = 0) as a function of the axial wavenumber k for various Reynolds numbers.

Figure 2 .Figure 3 .

 23 Figure 2. (a) Contribution of the output structures to the response of the system. (b) Contribution of the input structures to the excitation of the system. The case shown corresponds to m = 0, k = 1 and Re = 1000.

Figure 4 .

 4 Figure 4. (a) Influence of the Reynolds number on the energy amplification factor for k = 1 and m = 0. The cross markers represent calculations and the dashed-dotted line is the theoretical scaling derived from[START_REF] Antkowiak | Dynamique aux temps courts d'un tourbillon isolé[END_REF] analysis. For this axial wavenumber, the linear regression from the numerical data gives a gain scaling as G∞ = O(Re 2.77 ). (b) Scaling law for small k.

Figure 5 .

 5 Figure 5. Energy amplification in the helical case (m = 1) as a function of the axial wavenumber k for various Reynolds numbers.

Figure 6 .Figure 7 .

 67 Figure 6. Isocontours of axial vorticity for the dominant structures for m = 1 and Re = 1000. The same convention as in figure 3 is used except for (d) where ten equally spaced levels have been used. The dotted circle corresponds to the location of the maximum azimuthal velocity of the Lamb-Oseen vortex at r = 1.1209r0. First input structure for (a) k = 1.35 accounting for 31% of the energy amplification, (c) k = 2.5 and 17%, (e) k = 0.5 and 50%. First output structure for (b) k = 1.35 accounting for 32% of the flow excitation, (d) k = 2.5 and 16%, (f) k = 0.5 and 52%.

Figure 8 .

 8 Figure 8. Energy amplification reached for the statistical steady-state in the double-helix case (m = 2) as a function of the axial wavenumber k for various Reynolds numbers.

Figure 9 .

 9 Figure 9. Isocontours of axial vorticity for the first functions for m = 2 and Re = 1000. The same convention as in figure 3 is used. First input structure for (a) k = 2.5 accounting for 8% of the energy amplification, (c) k = 1 and 10%. First output structure for (b) k = 2.5 accounting for 8% of the flow excitation, (d) k = 1 and 10%.

† A specific notation, say A + , should have been used for the adjoint operator as its definition strictly depends on the chosen inner product. In the present case the adjoint operator A + coincides with the Hermitian transpose A H since the standard inner product of the Hilbert space H with identity weighting is used[START_REF] Hoepffner | Stability and control of shear flows subject to stochastic excitations[END_REF]).

† More precisely, as mentioned by one referee, a forcing of the form of a realistic turbulent signal would inject very low amplitude perturbations at large scales k → 0, and therefore a peak is expected at some finite wavelength in the response variance for each m. This wavelength selection will depend on the spectral content of the azimuthal Fourier decomposition of the turbulent signal, that might be highly sensitive to the particularities of the experimental environment in the long wave limit.

This work was supported by the European Community in the framework of the FAR-WAKE project under grant number AST4-CT-2005-012238 and the French National Research Agency as part of the VORTEX project.