R V K Chakravarthy

L Lesshafft

P Huerre

Global stability of buoyant jets and plumes

The linear global stability of laminar buoyant jets and plumes is investigated under the low Mach number approximation. For Richardson numbers in the range 10 -4 Ri 10 3 and density ratios S = ρ ∞ /ρ jet between 1.05 and 7, only axisymmetric perturbations are found to exhibit global instability, consistent with experimental observations in helium jets. By varying the Richardson number over seven decades, the effects of buoyancy on the base flow and on the instability dynamics are characterised, and distinct behaviour is observed in the low-Ri (jet) and in the high-Ri (plume) regimes. A sensitivity analysis indicates that different physical mechanisms are responsible for the global instability dynamics in both regimes. In buoyant jets at low Richardson number, the baroclinic torque enhances the basic shear instability, whereas buoyancy provides the dominant instability mechanism in plumes at high Richardson number. The onset of axisymmetric global instability in both regimes is consistent with the presence of absolute instability. While absolute instability also occurs for helical perturbations, it appears to be too weak or too localised in order to give rise to a global instability.

Introduction

Vertical injection of light fluid into a denser ambient creates a flow that either bears the characteristics of a jet, if the injected momentum is dominant over the buoyant forces, or those of a plume, if the momentum that is generated by buoyancy is dominant over the momentum that is imparted at the orifice.

The instability behaviour of jets is known to be strongly affected by density variations, even if buoyancy is not taken into account. [START_REF] Monkewitz | Absolute instability in hot jets[END_REF] found that jets at a jet-to-ambient density ratio below 0.72 in zero gravity display absolute instability, which leads to the self-sustained formation of ring vortices at a well-defined frequency. This phenomenon has been observed experimentally [START_REF] Sreenivasan | Absolute instability in variable density round jets[END_REF][START_REF] Monkewitz | Self-excited oscillations and mixing in a heated round jet[END_REF][START_REF] Boujemaa | Analyse spatio-temporelle de jets axisymétriques d'air et d'hélium[END_REF][START_REF] Hallberg | On the universality of global modes in low-density axisymmetric jets[END_REF] and numerically [START_REF] Lesshafft | Nonlinear global modes in hot jets[END_REF][START_REF] Nichols | Self-sustained oscillations in variabledensity round jets[END_REF]. Lesshafft & Huerre (2007) established that the absolute instability arises from non-buoyant baroclinic effects. [START_REF] Mollendorf | An experimental and numerical study of the viscous stability of a round laminar vertical jet with and without thermal buoyancy for symmetric and asymmetric perturbations[END_REF] included the action of buoyancy in the form of weak forcing terms in a local instability analysis. Recently, [START_REF] Coenen | Global instability of lowdensity jets[END_REF] performed a global instability analysis for light jets in the zero Mach number limit, achieving good agreement with the helium jet experiments of [START_REF] Hallberg | On the universality of global modes in low-density axisymmetric jets[END_REF]. A small Richardson number in these experiments characterises buoyant effects as being small, and the global analysis confirms that their impact on the instability behaviour is negligible in this regime.

The instability of plumes, at high Richardson numbers, has received far less attention from a theoretical perspective. [START_REF] Wakitani | The stability of a natural convection flow above a point heat source[END_REF] and [START_REF] Riley | On the stability of an axisymmetric plume in a uniform flow[END_REF] investigated the local instability characteristics, both temporal and spatial, of selfsimilar plumes within the limits of the Boussinesq approximation. Under the same assumption, [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] considered the convective/absolute nature of local instability in plumes, both in the self-similar regime far removed from the buoyancy source, and in the near-source region for one particular configuration. It was established that helical perturbations of azimuthal wavenumber m = 1 undergo a transition to absolute instability, due to a saddle point in the dispersion relation that is conditioned on the presence of buoyancy. However, the associated growth rates seem to be small, and their relevance for global and non-Boussinesq dynamics remains to be proven. The axisymmetric mode was found to be at most convectively unstable.

The instability of internal plumes in a confined domain appears to be a separate subject. In direct numerical simulations performed in the Boussinesq limit, [START_REF] Lopez | Instability of plumes driven by localized heating[END_REF] document a sequence of global state bifurcations in such closed flows, occurring at successive critical Rayleigh numbers. A linear global instability analysis of the same configuration [START_REF] Lesshafft | Linear global stability of a confined plume[END_REF] suggests that at least the first of these bifurcations arises due to pressure feedback between the top and bottom solid walls.

Numerous experiments have been performed on plumes with large density differences, where the Boussinesq approximation is not justified. [START_REF] Subbarao | Investigation of a co-flowing buoyant jet: experiments on the effect of reynolds number and richardson number[END_REF] conducted helium-air experiments, and they reported periodic axisymmetric puffing at Reynolds and Richardson numbers, Re and Ri, above critical values. Similar observations were made by [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF] for a larger range of Ri. A power law was obtained in the latter study that relates the puffing Strouhal number to Re and Ri. These experimental findings were corroborated by numerical simulations [START_REF] Jiang | Direct numerical simulation of the puffing phenomenon of an axisymmetric thermal plume[END_REF][START_REF] Satti | Computational analysis of gravitational effects in lowdensity gas jets[END_REF] and in additional recent experiments by [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF]. Through systematic variation of the gravity parameter, [START_REF] Satti | Computational analysis of gravitational effects in lowdensity gas jets[END_REF] demonstrated that the onset of self-sustained oscillations in their setting is contingent on the presence of gravity. The large majority of experiments and simulations suggest a dominant role of axisymmetric instability structures, contrary to the conclusions drawn from local instability analysis of self-similar Boussinesq plumes by [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF].

The present investigation addresses the linear instability of buoyant jets and plumes in a global and non-Boussinesq framework. The low Mach number approximation of [START_REF] Mcmurtry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF] is used in a form where density variations arise from heating of a single-species fluid. This formulation allows to examine the stability of buoyant jets and plumes on a continuous scale provided by the Richardson number, while the density ratio as an independent parameter characterises the departure from the Boussinesq condition. Special attention will be given to the physical origin of flow instability, by means of sensitivity analysis.

A similar approach has been pursued by [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF] in their analysis of helium plumes. That study demonstrated close agreement between the occurrence of self-excited puffing in experiments and the onset of global linear instability. Furthermore, the linear analysis was shown to accurately predict the puffing frequency, even far from the instability threshold.

The paper is organised in the following manner. Section 2 introduces the governing equations and the numerical procedures employed for the computation of base flows and their instability characteristics in global and local frameworks. Global instability results are presented in §3, followed by a discussion of the relevant physical mechanisms in §4. The global results are complemented by a local absolute/convective analysis in §5, which provides a link with the study by [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] of local instability in Boussinesq settings. Conclusions are offered in §6.

Problem formulation

Governing equations

A calorically perfect fluid is injected into an unstratified quiescent ambient of the same fluid at lower temperature, from a circular orifice in an adiabatic wall. In order to model a flow with strong density variations but negligible compressibility, a low-Mach-number approximation of the compressible Navier-Stokes equation is used. This approximation, which retains all effects of variable density in the convective terms, but discards the compressible dependency of density on pressure, was introduced by [START_REF] Mcmurtry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF] for a study of non-buoyant jets in the limit of zero Mach number. It was then extended to include a buoyancy term by [START_REF] Nichols | Self-sustained oscillations in variabledensity round jets[END_REF] and [START_REF] Chandler | Sensitivity analysis of low density jets and flames[END_REF], and their formulation is used in the present investigation. The dimensional governing equations in this approximation are given by

∂ ρ ∂t + div (ρũ) = 0, (2.1a) ρ Dũ Dt = -grad p + µ ∆ũ + 1 3 grad (div ũ) + g(ρ ∞ -ρ)e z , (2.1b) ρC p D T Dt = α∆ T , (2.1c) ρR T = p 0 , (2.1d )
where ρ, ũ, p, T denote the dimensional density, velocity, pressure and temperature, ρ∞ is the ambient density, g is the acceleration due to gravity, α is the thermal conductivity, C p is the specific heat, µ is the dynamic viscosity, and R is the specific gas constant. Note that the pressure in this formulation is split into a thermodynamic component p 0 , which is constant throughout the flow, and a fluctuating hydrodynamic component p. While the continuity and momentum equations (2.1a,b) are of the same form as in the fully compressible case, the energy equation (2.1c) simplifies to a simple advection-diffusion equation for temperature.

In dimensionless form, scaled with the nozzle radius R, the inlet centreline velocity ũj , the temperature difference Tj -T∞ between inflowing and ambient fluid, and the ambient density ρ∞ , equations (2.1) become (2.6)

∂ρ ∂t + div (ρu) = 0, (2.2a) ρ Du Dt = -grad p + 1 ReS ∆u + 1 3 grad (div u) + Ri S -1 (1 -ρ)e z , (
In the homogeneous limit S → 1, the state equation (2.2d) prescribes ρ → 1, and the system (2.2) is then identical to the Boussinesq equations used in [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF], provided ρ -1 is taken to be of order S -1; the present formulation is therefore consistent with our earlier study. A rigorous derivation of the Boussinesq equations and a discussion on the underlying assumptions may be found in Tritton (2012).

Base flow

In a cylindrical coordinate system (r, θ, z), the flow variables q = (ρ, u, p, T ) T are split into steady and unsteady components as q(r, θ, z, t) = q(r, θ, z) + q (r, θ, z, t).

(2.7)

The numerical tools used for the construction of the base flow, as well as for the linear perturbation analysis described in §2.3, are very similar to those employed by [START_REF] Coenen | Global instability of lowdensity jets[END_REF], except that density variations are modelled as an effect of heating, rather than species mixing. Equations (2.2) are discretised with finite elements in FreeFEM++, and a steady axisymmetric solution q is obtained from Newton-Raphson iterations [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF]. The numerical domain is 80 radii long in the streamwise direction and 30 radii in the transverse direction. Iterations are performed until all flow quantities are converged to within 10 -9 times their maximum values. At the inlet, z = 0, boundary conditions

u z = 1 2 + 1 2 tanh 5 2 1 r -r , u r = 0 and ρ = 1 -1 - 1 S u z (2.8)
are prescribed for the axial velocity u z , the radial velocity u r and the density ρ. The initial shear layer momentum thickness is 10% of the nozzle radius. All other boundary conditions are specified as

1 Re ∂u ∂r -pe r = 0, ρ = 1 at r = r max , (2.9a) 1 Re ∂u ∂z -pe z = 0, ∂ρ ∂z = 0 at z = z max , (2.9b) ∂u z ∂r = u r = ∂ρ ∂r = ∂p ∂r = 0 at r = 0.
(2.9c)

These boundary conditions are obtained from the kinematic constraints on the axis [START_REF] Khorrami | Application of spectral collocation techniques to the stability of swirling flows[END_REF], and by imposing zero normal stresses at r max and z max , together with Dirichlet and Neumann conditions for the density.

As the objective is to characterise the role of buoyancy in the instability dynamics, the main parameters to be varied are the Richardson number and the density ratio. The ranges of parameters 10 -4

Ri 10 3 and 1.05 S 7 will be investigated. The effect of the Reynolds number above a value of 100 is found to be weak, and a standard value of Re = 200 (in some cases Re = 500) will be used, while P r = 0.7 is maintained throughout.

For a strong density ratio S = 7, and the two extreme values Ri = 10 -4 and Ri = 10 3 , steady base flow profiles of axial velocity and density are shown in figure 1. The flow at low Ri is dominated by the momentum of the injected fluid, which diffuses radially with axial distance under the effect of viscosity. This is clearly the case of a jet, in a configuration where buoyancy has no noticeable impact on the base flow dynamics, despite the strong density variations. The flow at high Ri, on the contrary, is principally driven by the buoyancy force, as the injected momentum is negligibly weak in comparison. The fluid in this case is not pushed out of the orifice, but rather pulled out by buoyancy, forming a slender rising column around the axis (note the different radial scales in figure 1). This flow is characteristic of a plume, and it is often called a 'lazy' plume, as the momentum at its base is much lower than that of a self-similar profile, where momentum and buoyancy are in balance.

Between the two extremes shown in figure 1, the aspect of the base flow at different Richardson numbers changes gradually. Cases with Ri < 1 will be denoted 'buoyant jets' in the following, as opposed to 'plumes' with Ri > 1. Vertical variations of the centreline velocity are shown in figure 2 for the same two configurations as in figure 1. The jet, at Ri = 10 -4 , exhibits a short potential core region, where the centreline velocity is constant, followed by hyperbolic decay. The plume flow, at Ri = 10 3 , accelerates progressively with vertical distance from the inlet, and the centreline velocity approaches asymptotically a limit value in the self-similar regime [START_REF] Yih | Fluid mechanics[END_REF]. Buoyant jets at low but finite Richardson number behave as plumes at large distances from the nozzle, when their excess momentum is sufficiently diffused. While a jet entrains ambient fluid only through momentum diffusion, entrainment into a plume tends to be much stronger due to its continuous production of axial momentum. The plume base flow presented in figures 1(c,d ) is particularly marked by radial entrainment close to z = 0.

Linear stability problem

Infinitesimal perturbations of a steady base flow are sought with a global ansatz [ρ , u , p , T ] T = ρ(r, z), û(r, z), p(r, z), T (r, z) T e i(mθ-ωt) + c.c.

(2.10)

The integer m denotes the azimuthal wavenumber and ω = ω r + iω i is a complex frequency. Upon linearising the governing equations (2.2), and substitution of (2.10), the linear perturbation equations are found as

-iω ρ + div m (ρ u + ρ û) = 0, (2.11a) -iωρû + ρ (grad m u) • û + (grad m û) • u + ρ (grad m u) • u = (2.11b) -grad m p - Ri S -1 ρe z + 1 ReS [∆ m û + 1 3 grad m (div m û) , -iωρ T + ρ (grad m T ) • û + (grad m T ) • u + ρ (grad m T ) • u = 1 P rReS ∆ m T , (2.11c) ρ + ρ 2 T (S -1) = 0. (2.11d )
Differential operators in the above equations are written with a subscript m in order to indicate that azimuthal derivatives are replaced with a factor im; these operators are documented in the appendix. While the base flow is taken to be swirl free, u θ = 0, the azimuthal perturbation velocity u θ may in general be non-zero. Homogeneous Dirichlet conditions are imposed on û and ρ at the inlet, z = 0, and a homogeneous Neumann condition is prescribed for p. On the axis, depending on the azimuthal mode considered, appropriate boundary conditions as detailed in [START_REF] Khorrami | Application of spectral collocation techniques to the stability of swirling flows[END_REF] and [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] are enforced. On the outer boundaries at r max and z max , no-stress conditions consistent with (2.9) are used. The system (2.11) is cast into the form of an eigenvalue problem ωB q = Lq.

(2.12)

According to the ansatz (2.10), the real part of the eigenvalue, ω r , represents the oscillation frequency while the imaginary part, ω i , is the growth rate of the perturbation.

As before, the variable T is eliminated through the equation of state (2.11d). So-called global eigenmodes are computed by resolving ρ, û and p in both r and z, such that spatial variations of the base flow and the perturbation quantities are accounted for without further limiting assumptions [START_REF] Theofilis | Advances in global linear instability analysis of nonparallel and threedimensional flows[END_REF]. The eigenvalue problem (2.12) is then solved with an iterative shift-invert IRAM algorithm, in the same way as in [START_REF] Garnaud | Modal and transient dynamics of jet flows[END_REF], with an accuracy close to machine precision. In addition, a local analysis is performed in §5, in order to identify the absolute mode in a parallel base flow [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF]. Perturbations (2.10) are then Fouriertransformed in z, leading to the standard ansatz ρ(r, z), û(r, z), p(r, z), T (r, z)

T = ρ(r), ǔ(r), p(r), Ť (r) T e ikz .
(2.13)

Global spectra and eigenfunctions

Instability results obtained from the global formulation (2.10) are presented first for axisymmetric modes, m = 0, since experimental and numerical evidence suggests their leading role in the self-sustained dynamics of plumes, as discussed in §1. A brief discussion of helical perturbations, m = 1, follows in §3.2.

Axisymmetric perturbations

Eigenvalues ω obtained for two configurations, Ri = 10 -4 and 10 3 , with otherwise identical parameters S = 7 and Re = 200, are presented in figure 3 as black symbols. These are the two extreme jet and plume cases discussed in §2.2.

The jet, at Ri = 10 -4 (figure 3a), exhibits one unstable mode with ω = 0.558 + 0.025i. When buoyancy effects are eliminated, by setting Ri = 0 in the perturbation equations but still using the same base flow, eigenvalues shown as red crosses are obtained; visibly, the buoyancy term in the perturbation equations has no significant impact on the instability dynamics. This observation, as well as the overall appearance of the spectrum, is fully consistent with the non-buoyant and slightly buoyant results of [START_REF] Coenen | Global instability of lowdensity jets[END_REF], and many more details on the unstable eigenmode of pure jets are provided in their study.

The plume, at Ri = 10 3 (figure 3b), possesses five unstable modes, one being strongly dominant with ω = 29.4+11.3i. When the perturbation Richardson number is set to zero in this configuration, all unstable modes vanish from the spectrum. It can be concluded that buoyancy plays a determining role for the instability of this plume, not only by setting up the base flow, but also by the coupling of density and velocity perturbations. It is also noted that the high-Ri plume base flow is not subject to the non-buoyant instability that affects the low-Ri jet.

Unstable eigenvalues of the plume take on significantly higher values, both in their real and in their imaginary parts, than those of the jet. This is clearly a result of the scaling with the inflow centreline velocity, which appropriately characterises a jet, but is less pertinent for 'lazy' plumes. A common scaling is employed here for the sake of consistency across the full range of Richardson number values; if only high-Ri plumes were considered, a buoyancy-based velocity scaling would be more suitable.

The distinct nature of the instability modes of the jet and the plume flows is also apparent in the shape of their eigenfunctions. Figure 4 shows the spatial distribution of the axial velocity amplitude in the dominant modes for the two respective cases. At low Richardson number, the maximum perturbation amplitude is found on the centreline, 12 radii downstream of the nozzle. Perturbations are confined inside the jet column, as documented in more detail by [START_REF] Coenen | Global instability of lowdensity jets[END_REF]. In the high-Ri plume, the spatial eigenmode structure is very different: the maximum amplitude is located inside the mixing layer very close to the inflow, in the region where the density gradient in the base flow is maximal.

In the reference experiments by [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF] and [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF], a slightly different definition is chosen for the Richardson number, which corresponds to 2Ri/S in our nomenclature. Variations of the dominant eigenvalue in the present analysis are therefore presented as functions of Ri/S in figure 5 in order to facilitate a comparison. The baseline case, with Re = 200, S = 7 and inflow conditions (2.8), is represented by solid circles. Eigenvalues of this configuration display continuous variations both in the Strouhal number St = ω r /(2π) and in the growth rate ω i . Strouhal number values are asymptotically constant in the low-Ri/S regime, whereas they follow a power law at values Ri/S > 0.1. A regression fit yields the dependence St = 0.55(Ri/S) 0.43 , which is in good agreement with experimental results in the range 1

Ri/S 250: after conversion to the present definition of the Richardson number, the power law determined by [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF] in this regime is given by St = 0.52(Ri/S) 0.38 (shown as a line in figure 5a), and the corresponding measurements of [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF] for plumes from an orifice convert to St = 0.51(Ri/S) 0.39 . A different scaling, measured as St ∝ (Ri/S) 0.28 for Ri/S > 250 by [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF], is outside the parameter range considered here. The growth rate in the baseline configuration, shown in figure 5b, increases monotonically with Ri/S, and it is positive throughout.

Eigenvalues from three other flow configurations are included in figure 5 in order to assess the sensitivity of the instability with respect to the Reynolds number, to the density ratio and to the inlet velocity profile. With the standard profile (2.8), parameter combinations Re = 500, S = 7 (white circles) and Re = 200, S = 4.5 (triangles) are chosen. The Strouhal number values in subfigure a are barely affected by these changes, and the growth rates in subfigure b follow a similar trend as in the baseline case. The less heated configuration (triangles) is stable at Ri < 0.1, consistent with the analysis by [START_REF] Coenen | Global instability of lowdensity jets[END_REF].

Square symbols indicate results for a special case where Re = 200 and S = 7 are maintained, but the velocity inlet profile is changed to a parabolic pipe flow, while the density profile is still given by (2.8). This flow case is introduced in order to better approach the experimental conditions of [START_REF] Subbarao | Investigation of a co-flowing buoyant jet: experiments on the effect of reynolds number and richardson number[END_REF] and [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF], where the fluid exits from a nozzle as a developed laminar pipe flow. This change in the velocity profile barely has any effect on the Strouhal number across all Ri/S values, but it does inhibit the global instability in the low-Ri/S regime. The latter is consistent with the absence of self-excited behaviour at low Ri in the experiments of [START_REF] Subbarao | Investigation of a co-flowing buoyant jet: experiments on the effect of reynolds number and richardson number[END_REF] and [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF], and in the simulations of [START_REF] Satti | Computational analysis of gravitational effects in lowdensity gas jets[END_REF].

As the effect of heating enters the problem both through the density ratio S and through the Richardson number Ri, the onset of instability is examined for independent variations of these two parameters. The main results of the local (see §5) and global instability analyses are summarized in the state diagram in the (Ri,S) plane shown in figure 6. The thin line delineates the neutral boundary separating locally convectively unstable inlet base flows (in white below the curve) from locally absolutely unstable ones (in blue and red above the curve). The thick neutral line separates the globally stable states (in white and blue below the curve) from the globally unstable states (in red above the curve). In the white area, convective instability prevails and in the blue and red areas, absolute instability prevails. Along a diagonal line in the state diagram, the flow changes from convectively unstable to absolutely unstable to globally unstable, once the absolutely unstable region is large enough. For buoyant jets (low Ri) and plumes (large Ri) the same sequence takes place as S increases at a given Ri. Note that the globally unstable domain is reached much 'sooner' for plumes than for buoyant jets. The dip in the global stability neutral curve for Ri of order unity signals a gradual shift from a shear dominated instability to a buoyancy dominated instability, as further discussed in §4.

Helical perturbations

The local analysis of plumes in the Boussinesq whereas axisymmetric perturbations in that setting were found to be only convectively unstable. Although the global instability of axisymmetric eigenmodes in non-Boussinesq situations, as documented above, appears to be fully consistent with experimental and numerical observations of self-excited behaviour, the possibility of helical global instabilities remains to be explored.

Eigenvalues pertaining to helical instability modes are displayed in figure 7 for two different calculations, both for the same physical parameter setting S = 7, Ri = 1 and Re = 200. One case, represented by blue unfilled circles, was computed with the same boundary conditions as all previous results. A branch of regularly spaced modes is seen to be unstable over the interval 0.8 ω r 4. The features of this branch are very similar to several jet cases discussed by [START_REF] Coenen | Global instability of lowdensity jets[END_REF], as well as observations made in many different flow cases, especially in the boundary layer calculations by Åkervik et al. (2008). In a recent study [START_REF] Lesshafft | Artificial eigenmodes in truncated flow domains[END_REF], the occurrence of such 'arc branches' is ascribed to the presence of spurious pressure feedback from the downstream end of a truncated flow domain, and it is predicted that artificial damping near the outflow should be effective in reducing the growth rate of such unphysical modes.

Indeed, if artificial damping is applied in an 'absorbing layer' [START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF]) at z > 30, the growth rates of the arc branch modes decrease, and all modes recede to the stable half plane if the damping is sufficiently strong. Such a case is represented by red filled circles in figure 7, where a damping term -λ(z)q has been added to the right-hand side of equation (2.12). The damping coefficient λ(z) ramps up from zero to 16, over the interval 30 < z < 50, according to equation (2.4) of [START_REF] Chomaz | Fully nonlinear dynamics of parallel wakes[END_REF]. No unstable helical modes are found with this boundary treatment. In contrast, the same artificial damping has no discernible effect on the unstable eigenvalues for axisymmetric perturbations shown in figure 3. This behaviour is in full agreement with the arguments of [START_REF] Lesshafft | Artificial eigenmodes in truncated flow domains[END_REF], as m = 0 perturbations are locally stable in the downstream flow region, whereas m = 1 perturbations remain convectively unstable, as will be shown in §5.

Influence of buoyant, baroclinic and shear-related mechanisms

The results discussed in §3.1 suggest that different mechanisms drive the global instability dynamics in the low-and the high-Ri regimes. These mechanisms are investigated in the present section, on the basis of the formalism proposed by [START_REF] Marquet | Identifying the active flow regions that drive linear and nonlinear instabilities[END_REF]. This formalism is introduced here in a slightly different and weaker manner, which is sufficient for the present purpose.

A sensitivity analysis is to be performed, in order to quantify the influence of the various forces in the momentum equation onto the unstable growth of velocity perturbations. The latter are governed by the equation

-iωû = C + S + P + B + V, (4.1)
with the right-hand side terms

C = -(grad û) • u, (4.2a) S = -(grad u) • û, (4.2b) P = - grad p ρ + ρ grad p ρ 2 , (4.2c) B = - Ri S -1 ρ ρ 2 e z , (4.2d ) V = 1 ReS 1 ρ ∆û + grad(div û) 3 - ρ ρ 2 ∆u + grad(div u) 3 . (4.2e)
These individual terms represent the mechanisms of base flow convection C, base flow shear S, pressure force P, buoyancy B, and viscous diffusion V. As only axisymmetric m = 0 perturbations are considered in this chapter, it is not necessary to append a subscript m to the differential operators. For a better physical discussion, the pressure force can be split into a a divergence-free (baroclinic) and a curl-free (barotropic) component; the former is linked to the baroclinic torque in the vorticity equation, after application of the curl operator to (4.1), whereas the latter does not affect the evolution of perturbation vorticity.

A Helmholtz decomposition is performed on the pressure force P, such that

P = curl A -grad φ, (4.3)
where A and φ are found from

A = 1 4π (curl P) ⊗ (1/r), (4.4a) φ = 1 4π (div P) ⊗ (1/r). (4.4b)
The operator ⊗ denotes a convolution over the entire volume V of the numerical domain, and r represents any position in V . This decomposition is performed numerically, such that P = P 1 +P 2 is explicitly obtained, with a baroclinic component P 1 and a barotropic component P 2 .

Further analysis is restricted to the action of shear, baroclinic and buoyant forces, because all other contributions are found to be strictly stabilising at all Richardson numbers. The dependence of an eigenvalue ω on these three components is obtained by introducing small variations into (4.1),

-iωû = C + (1 + ε S )S + (1 + ε P )P 1 + P 2 + (1 + ε B )B + V, (4.5)
from where sensitivities can be defined as

∂ S ω = ∂ω ∂ε S = q † , S q † , B q , ∂ P ω = ∂ω ∂ε P = q † , P 1 q † , B q , ∂ B ω = ∂ω ∂ε B = q † , B q † , B q . (4.6)
Note that the terms S, P 1 and B contain components of the eigenvector q, and that q † is the associated adjoint eigenvector, defined with respect to the scalar product •, • . A standard non-weighted discrete scalar product has been chosen in the present calculations, but the scalar quantities ∂ω in (4.6) are independent of this choice, as demonstrated by [START_REF] Marquet | Identifying the active flow regions that drive linear and nonlinear instabilities[END_REF]. The sensitivities (4.6) are interpreted in the following way. An infinitesimally small positive value ε S proportionally increases the strength of the shear-related force term, resulting in an eigenvalue variation δω = ε S ∂ S ω. If the imaginary part of ∂ S ω is positive, then S has a destabilising effect; if it is negative, then S acts in a stabilising way. The same reasoning applies to ∂ P ω and ∂ B ω. The three sensitivities are commensurate, so that a larger absolute value of one compared to another indicates a stronger stabilising or destabilising effect.

Results from this analysis are presented in figure 8 over the full range of Ri values, for the standard setting S = 7 and Re = 200. Imaginary values of ∂ S ω, ∂ P ω and ∂ B ω are shown in two separate diagrams for low and high Richardson numbers, for better readability. In the jet regime Ri < 1, the strongest destabilising force is due to the base flow shear. At very low Ri, the effect of buoyancy vanishes, while the baroclinic force provides a small additional destabilisation. The local analysis of Lesshafft & Huerre (2007) demonstrated that the baroclinic torque is the determining ingredient that renders a non-buoyant heated jet absolutely unstable, through co-operation with the basic shear instability. The present global results are consistent with that conclusion. In the plume regime Ri > 1, the buoyancy force becomes strongly destabilising, dominating all other contributions for Ri > 5. Shear and baroclinic effects are negligible in comparison at very high Ri, the baroclinic force even becomes stabilising above Ri = 100.

It is concluded that the observed global instability in the jet and plume regimes indeed

10 -4 10 -3 10 -2 10 -1 10 0 0 0.2 0.4 0.6 Ri ∂ S ω i , ∂ P ω i , ∂ B ω i a) 10 0 10 1 10 2 10 3 0 10 20 Ri ∂ S ω i , ∂ P ω i , ∂ B ω i b)
Figure 8: of the growth rate with respect to ( ) shear, (•) baroclinic and ( ) buoyancy terms, as functions of the Richardson number, with parameters S = 7 and Re = 200.

involve distinct physical mechanisms. In buoyant jets at low Richardson number, the dynamics are driven by a shear instability, which is strengthened by a baroclinic force. In high-Ri plumes, the instability arises principally from buoyancy effects. These conclusions are fully consistent with the results of [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF], who observed that the leading eigenmode could be stabilised through artificial compensation of the baroclinic torque at low Ri, and through suppression of the buoyancy force at high Ri.

Local analysis

The results so far have shown a dominance of axisymmetric global instabilities, which is in stark contrast to our earlier local analysis in the S → 1 Boussinesq limit [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF], where absolute instability was found only for helical perturbations. The absolute/convective nature of local instability in non-Boussinesq settings is now investigated.

The same algorithm as in [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] is used for tracking saddle points of the local dispersion relation in the complex k plane. Again, the group-velocity rootfinding procedure of [START_REF] Lesshafft | Optimal velocity and density profiles for the onset of absolute instability in jets[END_REF] has been found to be more efficient and robust than the classical Briggs or the cusp map methods (see [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF].

Axisymmetric perturbations, m = 0

For the standard setting S = 7 and Re = 200, and for Richardson numbers of 10 -4 and 10 2 , variations of the absolute frequency ω 0 and of the absolute wavenumber k 0 along the streamwise z direction are shown in figure 9. In the more extreme case Ri = 10 3 , numerical difficulties led to unreliable results. In both configurations, the flow is found to be absolutely unstable (ω 0,i > 0) over a streamwise interval of 6 or 7 radii downstream of the inlet. In the weakly non-parallel case of Ri = 10 -4 , the real part ω 0,r shows moderate variations around a value of 0.5, in reasonable agreement with the global frequency ω r = 0.56 as given in figure 5a. In the strongly non-parallel setting Ri = 10 2 , ω 0,r displays a variation between 1.3 and 80 over the absolutely unstable interval. This is not inconsistent with the global frequency ω r = 10.99, but it does not provide a means of predicting ω r at leading order. Nonetheless, the link between local absolute and global instability of axisymmetric perturbations is very plausibly confirmed by these results, both in the low-and in the high-Richardson number regime. The neutral curve for the onset of absolute instability in the (Ri,S) plane is reported in figure 6 for axisymmetric perturbations at z = 0. It is found that the transition from convective to absolute local instability at the inlet, with increasing Ri and S, occurs before global instability sets in. This is consistent with the common observation, both in model problems [START_REF] Chomaz | A frequency selection criterion in spatially developing flows[END_REF] and in slowly varying flow (e.g. Lesshafft et al. 2007), that absolute instability must prevail over a sufficiently long streamwise region with sufficient growth rate in order to prompt a global instability.

Helical perturbations, m = 1

The saddle point in the complex k-plane that is associated with helical absolute instability in the study of [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] is also recovered in the analysis of the present inflow profiles. The absolute growth rate ω 0,i for m = 1 perturbations at z = 0 is displayed in figure 10a density ratio S characterises the departure from the Boussinesq limit, and it is seen to have a very weak effect on the growth rate ω 0,i of the absolute helical mode.

The spatial variation of ω 0,i , over a short interval of z adjacent to the inlet, is shown in figure 10b for the highly non-Boussinesq setting S = 7; values beyond this streamwise region could not be obtained with sufficient confidence, due to numerical difficulties. The absolute helical growth rate at S = 7 displays the same qualitative characteristics as the one described in [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] for the Boussinesq limit: the growth rate declines quickly downstream of the inlet, perhaps asymptotically tending towards zero. It is not surprising that under these conditions the weak absolute instability of m = 1 perturbations does not lead to a global flow destabilisation, as was found in §3.2.

Conclusions

The linear instability dynamics in spatially developing buoyant jets and plumes have been investigated for a wide range of values of the Richardson number and the fluid density ratio. In the limit of zero Mach number, all variable-density effects have been taken into account in the mathematical formulation, so that configurations outside the realm of validity of the Boussinesq approximation could be considered.

Axisymmetric global instability modes have been found and documented over the entire investigated range of parameters, whereas no helical global instability could be identified. Some doubts remain only in the very high Richardson number regime, Ri > 100, where spurious helical modes could not be entirely stabilised due to numerical limitations. The preponderance of axisymmetric instability is in agreement with experimental observations by [START_REF] Subbarao | Investigation of a co-flowing buoyant jet: experiments on the effect of reynolds number and richardson number[END_REF], [START_REF] Cetegen | Experiments on the oscillatory behavior of buoyant plumes of helium and helium-air mixtures[END_REF] and [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF], who reported axisymmetric puffing in free plumes. The present global analysis furthermore reproduces the experimentally measured frequencies with satisfactory accuracy. In particular, the experimental power law ω r ∝ (Ri/S) 0.38 has been retrieved as ω r ∝ (Ri/S) 0.43 in the present calculations. It is noted that [START_REF] Bharadwaj | Global instability analysis and experiments of buoyant plumes[END_REF] found even closer agreement from their linear analysis, which was designed to specifically model helium plumes, as opposed to thermal plumes in the present study.

The physical mechanisms behind global instability have been characterised by means of a sensitivity analysis. As proposed by [START_REF] Marquet | Identifying the active flow regions that drive linear and nonlinear instabilities[END_REF], the sensitivity of the perturbation growth rate with respect to individual terms in the linear equations has been evaluated, which provides a quantitative measure for the stabilising or destabilising role of mechanisms represented by these terms. The main conclusion is that instability in the low-Ri jet regime is caused primarily by a shear mechanism, aided by a baroclinic force that arises from density variations, whereas the instability in the high-Ri plume regime is brought about principally by way of the buoyancy force, with a small contribution from shear. Nothing in the results indicates an abrupt switching between two distinct instability modes; the most unstable mode appears instead to vary smoothly as a function of Ri, with a gradual shift from the dominance of one mechanism to a dominance of the other.

These global results contrast with the conclusions of [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF], on the basis of local analysis in the Boussinesq limit, that the instability dynamics of selfsimilar plumes are dominated by helical perturbations. In particular, that earlier study showed axisymmetric perturbations to be only convectively unstable, whereas helical perturbations exhibit absolute instability, as well as larger temporal growth at real wavenumbers than axisymmetric modes. The global analysis in the present study was performed on base flows that develop from an orifice with prescribed inlet profiles, and that only relax asymptotically in the streamwise direction towards a self-similar flow solution. Close to the orifice, these base flows are markedly different from self-similar conditions, and this is the flow region where unstable axisymmetric perturbations reside in high-Ri plumes, according to the present results (see figure 4b). Global instability in low-Ri buoyant jets has been shown to depend on baroclinic effects, which are absent in the Boussinesq approximation. The Boussinesq framework used by [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF] is therefore inappropriate for an instability analysis in this regime.

It has finally been demonstrated, for selected configurations, that the global instability characteristics are consistent with the absolute or convective nature of local instability. All globally unstable settings in the (Ri, S) plane, with Re = 200, feature an absolutely unstable flow region in the vicinity of the inflow. In the examined cases, the absolute growth rate of axisymmetric perturbations is significantly larger than that of helical perturbations, and it remains positive over a longer streamwise region. Although absolute instability also arises for helical perturbations, it appears to be too weak to set off a global instability.

It must be cautioned that the conclusions drawn from the present results may not be easily extendable to generic plume and jet flows. In particular, the instability characteristics seem to be rather sensitive to details of the inflow profiles: with similar but not identical inflow profiles, axisymmetric perturbations are absolutely unstable in the present settings, but convectively unstable in the configuration of [START_REF] Chakravarthy | Local linear stability of laminar axisymmetric plumes[END_REF]. Test calculations, documented in [START_REF] Chakravarthy | Local and global instabilities in buoyant jets and plumes[END_REF], indicate that the functional shape of the density profile has a marked influence on the local stability characteristics, even when the mixing layer thickness is matched. [START_REF] Subbarao | Investigation of a co-flowing buoyant jet: experiments on the effect of reynolds number and richardson number[END_REF] point out, for instance, that helium release and diffusion flames create plumes with very distinct density variations, which therefore may present quite different instability dynamics. It can also not be ruled out that nonlinear effects alter the threshold of global instability [START_REF] Couairon | Absolute and convective instabilities, front velocities and global modes in nonlinear systems[END_REF]. The influence of the Prandtl number has not been investigated in this study, but it has been shown by [START_REF] Lakkaraju | Effects of Prandtl number and a new instability mode in a plane thermal plume[END_REF] that the instability behaviour of planar plumes undergoes qualitative changes as P r is varied far from unity.

The authors are grateful to Wilfried Coenen for his aid with the adjoint implementation in FreeFEM++, to Olivier Marquet for his contribution to the sensitivity analysis, and

Figure 1 :Figure 2 :

 12 Figure 1: Axisymmetric steady base flows obtained for Ri = 10 -4 (a,b) and for Ri = 10 3 (c,d ). Axial velocity (a,c) and density (b,d ) are shown in the (r, z) plane.

Figure 3 :

 3 Figure 3: Global spectra of axisymmetric perturbations (m = 0), for the two configurations shown in figure 1. a) Ri = 10 -4 , jet case; b) Ri = 10 3 , plume case. Re = 200 and S = 7 are set identically for both cases. True eigenvalues (•) are compared to their counterparts (+) that are obtained when the buoyancy term in (2.11b) is removed.

Figure 4 :

 4 Figure 4: Spatial distributions of axial velocity eigenfunctions ûz , corresponding to the most unstable modes of the two respective cases shown in figure 3. a) Jet at Ri = 0.0001, S = 7 and Re = 200; b) plume at Ri = 1000, S = 7 and Re = 200.

Figure 5 :

 5 Figure 5: a) number and b) growth rate of the global eigenvalue ω, as a function of Ri/S. Legend: (•) Re = 200, S = 7; (•) Re = 500, S = 7; ( ) S = 4.5, Re = 200; ( ) Re = 200, S = 7, with parabolic inlet velocity profile. All other configurations take the inlet velocity profile (2.8). ( ) Power law from the Cetegen & Kasper (1996) experiments, rescaled to match the present definition of Ri.

Figure 6 :

 6 Figure 6: Instability regimes in the S-Ri plane for the axisymmetric mode at Re = 200. The thick line denotes the global stability boundary, and the thin line denotes the boundary between local convective and absolute instability of the inlet profile.

Figure 7 :

 7 Figure 7: Global spectrum of helical perturbations (m = 1), for the configuration S = 7, Ri = 1 and Re = 200. Results from two calculations are shown: (•) without absorbing layer; (•) with absorbing layer at z > 30.

Figure 9 :

 9 Figure 9: Variation of the absolute frequency ω 0 with streamwise location z for the axisymmetric mode m = 0 and Ri = 10 -4 (left) and Ri = 10 2 (right) at S = 7, Re = 200.

Figure 10 :

 10 Figure 10: growth rate ω 0,i of helical perturbations, for Re = 200 and Ri = 1, (a) at the inlet, as a function of the density ratio S; (b) as a function of z, for S = 7.

R. V. K. Chakravarthy, L. Lesshafft and P. Huerre to Prof. Debopam Das for precious comments on the comparison with experiments. This work has been funded by the DGA under grant number LADX2331 and by the Agence Nationale de la Recherche under the "Cool Jazz" project. RVKC was supported by a PhD fellowship from École polytechnique.

Appendix A. Differential operator definitions for azimuthally decomposed field quantities

The operators for divergence, gradient, Laplacian and advection in cylindrical coordinates in equations (2.11) are written with a subscript m. This is meant to express that azimuthal derivatives of perturbations (2.10) are included in these operators in the form of terms in m, such that formally

and accordingly for all other flow variables. All relevant terms from the equations (2.11) are written out below for the sake of completeness.