A Rabinovich

O M Umurhan

N Harnik

F Lott

Eyal Heifetz
email: eyalh@post.tau.ac.il

Vorticity inversion and action-at-a-distance instability in stably stratified shear flow

The somewhat counter-intuitive effect of how stratification destabilizes shear flows is reexamined, in what we believe to be the simplest example, in terms of action-at-a-distance interaction between 'buoyancy-vorticity gravity wave kernels'. The setup consists of an infinite uniform shear Couette flow in which the Rayleigh-Fjortoft necessary conditions for shear flow instability are not satisfied. When two stably stratified density jumps are being added, the flow however becomes unstable. At each density jump the perturbation can be decomposed into two coherent gravity waves propagating horizontally in opposite directions. We show, in detail, how the instability results from a phase locking actionat-a-distance interaction between the four waves (two at each jump), but can as well be reasonably approximated only by the interaction between the two counter-propagating waves (one at each jump). From this perspective the nature of the instability mechanism is similar to the barotropic and baroclinic ones.

Next we add a small ambient stratification to examine how the critical level dynamics alters our conclusions. We find that strong vorticity anomaly is generated at the critical level due to the persistent vertical velocity induction by the edge waves at the jumps. This critical level anomaly acts in turn at-a-distance to decay the edge waves at the jumps.

When the ambient stratification is increased, so that the Richardson number exceeds the value of a quarter, this destructive interaction overwhelms the constructive interaction between the edge waves and consequently the flow becomes stable. This effect is manifested when considering the different action-at-a-distance contributions to the energy flux divergence at the critical level. The edge wave interaction is found to contribute toward positive divergence, that is, toward instability, whereas the critical level -edge waves interaction contributes toward an energy flux convergence, that is, toward stability.

Introduction

Recently Harnik et al. (2008) developed a 'buoyancy-vorticity wave interaction approach' to describe the linear dynamics of stably stratified shear flow. The essence of their view is that in stably stratified flows the buoyancy acts as a restoring force and therefore a horizontal (say zonal) gradient of material displacement generates vorticity in the meridional direction. By inversion such an initial localized vorticity anomaly at some particular layer induces a non-local vertical velocity field which deforms the flow at remote layers, generating vorticity in the far field. The latter, in turn, induces a far field velocity which deforms the initially perturbed layer and consequently alters the initial vorticity anomaly (Fig. 1). Thus, even in the absence of a mean vorticity gradient vorticity inversion and action-at-a-distance are at the heart of stratified shear flow dynamics.

In barotropic and baroclinic shear flows the basic (potential) vorticity building blocks which interact at a distance are the Counter-propagating Rossby Waves -CRWs (Bretherton (1966), Hoskins et al. (1985)). In stratified shear flows the basic building blocks at each level are mixed gravity-Rossby wave kernels that, in the absence of shear curvature, are reduced to two oppositely propagating internal gravity wave vortex sheets. The fact that at each given level two building blocks exist (rather than one CRW kernel in the barotropic/baroclinic case) complicates the dynamics, however this is a direct consequence of the interplay between the buoyancy and vorticity fields.

The two necessary conditions for modal instability of barotropic and baroclinic shear flows are straightforwardly rationalized from the CRW perspective. Mutual amplification between two CRWs are possible only if the mean vorticity gradient at their vicinities changes sign (Charney and Stern (1962) condition). The Fjortoft (1950) condition states that the mean vorticity gradient should be positively correlated with the mean velocity profile, as measured with respect to the point of zero vorticity gradient. When this condition is satisfied the CRWs are able to counter-propagate against the shear, which allows then to phase lock in their constructive interaction configuration. It is thus easy to understand why these conditions are irrelevant to stratified shear flow. The mechanism of vorticity anomaly generation by horizontal buoyancy gradients is independent of the presence of mean vorticity gradients. Furthermore, at every layer one of the two gravity wave kernels (propagating in opposite directions in the mean flow frame of reference) is always counter-propagating against the mean flow. The Charney -Stern and Fjortoft conditions result, respectively, from the global conservation of pseudo-momentum and pseudo-energy (Shepherd 1990). The equivalent Miles (1961) -Howard (1961) (M-H) necessary condition for stratified shear flow instability (that the Richardson number should be smaller than a quarter somewhere within the domain) does not seem to result from an eddy associated global conservation law, and therefore it is more challenging to interpret.

One of the main aims of this study is to gain a deeper insight on this condition from the Kernel Gravity Wave (KGW) interaction standpoint.

The mechanistic picture of the simple barotropic Rayleigh (1880) and the baroclinic Eady (1949) models, in terms of a two-CRW interaction (Heifetz et al., 1999;Davies and Bishop, 1994), served as a guiding line in much more complex shear flows (Heifetz et al. 2004), including the nonlinear eddy life cycles on realistic jets (Methven et al. 2005). Here we follow the same route and examine in detail an equivalent simple model for stratified shear flow of two density jumps in the presence of constant shear1 (Fig. 2a). After formulating the problem and briefly reviewing the buoyancy-vorticity KGW interaction approach (Section 2), we analyze the four KGW interaction of this problem (Section 3). As expected the instability is dominated by the two counter-propagating kernels, however for long wavelengths the other two pro-propagating waves, which must be hindered, contribute as well to the instability.

The common interpretation for the M-H condition stems from the over-reflection perspective, which is based on a cross-shear wave propagation view. Lindzen and Barker (1985) found that the wave geometry which allows over-reflection, and hence instability, must include a critical level (hereafter CL) which separates wave propagating regions from wave evanescence ones. To obtain an evanescence region near the vicinity of a CL the local Richardson number there must be less than a quarter (Booker and Bretherton 1967). When analyzing the index of refraction of wave propagation across the shear in the Taylor-Goldstein equation Lindzen and Barker (1985) found that the wave geometry which allows over-reflection, and hence instability, must include a critical level (hereafter CL) which separates between wave propagation and wave evanescence regions. To obtain an evanescence region at the vicinity of a CL one must have local values of the Richardson number which are smaller than a quarter (Booker and Bretherton 1967). Alternatively, Baines and Mitsudera (1994, see also Baines text book chapter 4, 1998) used a wave interaction approach, similar to our KGW interaction, but they interpret the suppressing of the instability in terms of absorption of the two waves at the CL between them (a cross shear wave propagation view). Here we interpret this stabilizing effect of the CL purely in terms of counter-propagating kernel-wave interaction.

The M-H condition is obviously satisfied in the simple setup of the two density jumps, since the Richardson number is concentrated in two delta functions there. When increasing slightly the stratification elsewhere, so that the ambient Richardson number increases from zero to a quarter the instability is suppressed, in agreement with the M-H condition.

Nonetheless, it seems counter-intuitive that such weak ambient stratification (a value of quarter compared to delta functions) will dominate over the dynamics of the density jumps. In Section 4 we analyze the destructive effect of the CL consistently from the buoyancy-vorticity action-at-a-distance approach. It is shown how the KGWs at the density jumps generate vorticity at the CL which acts both to destroy itself as well as to mask the constructive interference between the KGWs in the far field. We discuss these results in Section 5.

KGW formulation for stratified shear flow

We consider an inviscid, incompressible, Boussinesq, 2D flow in the zonal-vertical (x-z) plane, with a zonally uniform basic state that varies only with height and is in hydrostatic balance. For this setup the linearized momentum and continuity equations become:

Du Dt = -U z w - 1 ρ 0 ∂p ∂x (1a) Dw Dt = b - 1 ρ 0 ∂p ∂z (1b) Db Dt = -wN 2 (1c) ∂u ∂x + ∂w ∂z = 0 (1d)
U is the mean zonal velocity where the mean shear, U z , is taken to be constant in this study. The linearized material derivative is D Dt = ∂ ∂t + U ∂ ∂x , (u, w) are the perturbation velocities in the x and z directions and p is the perturbation pressure. ρ 0 is a constant reference density, b = -( ρ ρ 0 )g is the perturbation buoyancy, where ρ is the density perturbation, and g being the gravity. The Brunt-Väisälä frequency N 2 = -g ρ 0 dρ dz = b z is the mean vertical buoyancy gradient where ρ is the mean density profile.

Equation set (1) yields, after some algebra, the vorticity equation:

Dq Dt = ∂b ∂x = -b z ∂ζ ∂x (2)
where q = ∂w ∂x -∂u ∂z is (minus) the perturbation component of the vorticity perpendicular to the zonal-vertical plane and ζ is the vertical displacement. Equation (2) simply states that since the buoyancy is a restoring force its tendency to flatten back the deformed isopycnals generates vorticity to the right of the perturbation buoyancy gradient (Fig. 1, see also Fig. 3). By inversion, q induces a far field vertical velocity that generates vertical displacement which, in turn, generates vorticity in the far field. To formulate this mathematically we look at a zonal Fourier component of the form of e ikx , with wavenumber k, and introduce its streamfunction ψ(k, z) satisfying: u = -∂ψ ∂z , w = ikψ, q = -k 2 ψ + ∂ 2 ψ ∂z 2 . Then the vertical velocity at some height z can be expressed in terms of the integrated contribution of the far field vorticity field

w(z) = z ′ q(z ′ )G(z, z ′ )dz ′ (3)
where G(z, z ′ ) is the Green function satisfying -k 2 G+ ∂ 2 G ∂z 2 = ikδ(z -z ′ ) with the appropriate boundary conditions. Here we apply the formulation to open flow whose perturbation vanishes at ±∞ so that G(z, z ′ ) = -i 2 e -k|z-z ′ | . Writing then

Dζ Dt = w = z ′ q(z ′ )G(z, z ′ )dz ′ (4)
formulates, together with (2), the mathematical description for the action in a distance interaction between the buoyancy b = -b z ζ, and the vorticity q.

The simplest nontrivial setup of vorticity-buoyancy interaction is of a single stably stratified density jump:

ρ = ρ 0 -|∆ρ|H(z) (5)
where H(z) is the Heaviside function and the density interface is at z = 0. The corresponding mean buoyancy gradient is then:

b z = - g ρ 0 dρ dz = g|∆ρ| ρ 0 δ(z) = ∆bδ(z) (6)
For finite vertical displacement ζ, (2) implies that the vorticity perturbation should be a delta function as well, q = q(x, t)δ(z). Equations ( 2) and ( 4) can be then written in the matrix form:

∂ ∂t q ζ = -ik U 0 ∆b 1 2k U 0 q ζ (7)
where U 0 is the mean velocity at z = 0. The solution of this simple eigen-problem is a superposition of the two Boussinesq internal gravity waves with the intrinsic phase speeds:

c ± = ± ∆b 2k ≡ ±c gr (8) q ζ = ζ + 0 2kc gr 1 e -ik(cgr+U 0 )t + ζ - 0 -2kc gr 1 e -ik(-cgr+U 0 )t (9)
where ζ ± 0 are the vertical displacements at z = 0 of the two modes. Hence, the vorticity and the vertical displacement of the eastward (westward) propagating wave, in the mean flow frame of reference, are in (anti) phase, ensuring coherent propagation of the buoyancy and the vorticity fields (Fig. 3). Next, we use these two waves as kernel gravity wave (KGW) building blocks, to examine the interaction between two density jumps (see Harnik et al., 2008, for more details on the KGW formulation).

3 Two density jump instability

Formalism

We consider the density profile of two density jumps located at z = ±h (Fig. 2a):

ρ = ρ 0 + |∆ρ|[1 -H(z -h) -H(z + h)] (10)
Normalizing length by h and time by 1/U z , equations ( 2) and (4) become:

( ∂ ∂t + ikz)q = -ikRiζ (11a) ( ∂ ∂t + ikz)ζ = - i 2 z ′ q(z ′ )e -k|z-z ′ | dz ′ (11b)
where all variables are non-dimensionalized and Ri = b z /(U z ) 2 is the Richardson number.

For the two density jump profile

Ri = ∆b h(U z ) 2 [δ(z -1) + δ(z + 1)] ≡ Ri[δ(z -1) + δ(z + 1)] ( 12
)
Denoting hereafter level z = 1 with the subscript 1 and z = -1 with 2, so that

q = q1 δ(z -1) + q2 δ(z + 1) (13)
we write equation set (11) explicitly for the two jumps:

( ∂ ∂t + ik)q 1 = -ik Riζ 1 (14a) ( ∂ ∂t + ik)ζ 1 = - i 2 (q 1 + q2 e -2k ) (14b) ( ∂ ∂t -ik)q 2 = -ik Riζ 2 (14c) ( ∂ ∂t -ik)ζ 2 = - i 2 (q 1 e -2k + q2 ) (14d)
and the intrinsic gravity wave speed in ( 8) is normalized to c gr = Ri/2k. By definition the two KGWs at each interface must preserve their eigen-structure of (9):

q± 1/2 = ± 2k Riζ ± 1/2 (15)
where as before the ± superscript indicates the sign of the KGW associated intrinsic propagation speed. Since the vorticity and the displacement at each interface are superpositions of the two KGWs (q 1/2 = q+

1/2 + q- 1/2 , ζ 1/2 = ζ + 1/2 + ζ - 1/2
), equation set ( 14) can be rearranged to be written as an equation set for the KGWs displacements:

∂ ∂t + ik(c gr + 1) ζ + 1 = -iσ(ζ + 2 -ζ - 2 ) (16a) ∂ ∂t + ik(-c gr + 1) ζ - 1 = -iσ(ζ + 2 -ζ - 2 ) (16b) ∂ ∂t + ik(c gr -1) ζ + 2 = -iσ(ζ + 1 -ζ - 1 ) (16c) ∂ ∂t + ik(-c gr -1) ζ - 2 = -iσ(ζ + 1 -ζ - 1 ) (16d)
Equation set ( 16) is straightforward to interpret. The mean normalized velocity at the

interfaces is U (±1) = ±1,

Results

The eigen-value solution of ( 16), for normal modes of the form of e ik(x-ct) , yields the dispersion relation

c = ± c 2 gr + 1 ± 2c gr 1 + σ k 2 (17)
Without interaction (σ = 0) the solution is reduced to the four uncoupled modes: c = ±(c gr ± 1). For σ = 0 instability is obtained when 2c gr 1 + σ k 2 > c 2 gr + 1. The unstable modes are stationary, c r = 0, experiencing a growth rate of

kc i = k -(c 2 gr + 1) + 2c gr 1 + σ k 2 .
Figure 4a shows the growth rate as a function of (k, Ri). As can be seen a global maximum for the growth rate is obtained when (k, Ri) ≈ (0.5, 1). The wavenumber range of instability shrinks as Ri increases and is shifted to higher wavenumbers. For Ri smaller than unity only a short wave cutoff exists whereas for Ri larger than unity both short and long-wave cutoffs exist. The unstable and stable modes come in conjugate pairs. The growth rate of these modes as well as the real phase speed of the four modes are shown for representative values of Ri in the two left columns of Fig. 5.

The amplitude ratios and the phase differences between the four KGWs for the unstable modes (and its continuation to the neutral regions at the two sides of the cutoffs)

are shown in the two right columns of Fig. 5 (normalized with respect to ζ - 1 ). In the unstable region the two counter-propagating waves (ζ - 1 , ζ + 2 ) have equal amplitudes. The pro-propagating waves (ζ + 1 , ζ - 2 ) are also of equal amplitude, smaller by a factor χ (Fig. 4b).

Each counter-propagating wave is aligned in anti-phase (both in vorticity and displacement) with the pro-propagating wave of the opposed level so that

(ζ + 1 , ζ - 2 ) = -χ(ζ + 2 , ζ - 1 ).
The two counter-propagating waves are in phase in their displacement (and anti-phased in their vorticity) at the long-wave cutoff and anti phased in their displacement (and in phase in their vorticity) at the shortwave cutoff. In between, the upper counter-propagating wave displacement (vorticity) is shifted eastward (westward) with respect to the lower counter-propagating one. A schematic of the KGW setup for the case where the counterpropagating waves are in quadrature is drawn in Fig. 6. Next we wish to interpret these results from the KGW interaction point of view.

The first question that might arise is why the pro-propagating kernels are present at all in the growth mechanism. For comparison, in the analogous barotropic Rayleigh and baroclinic Eady models, only counter-propagating Rossby waves exist and they are phase locked to yield unstable modes, similar to the counter-propagating waves (ζ - 1 , ζ + 2 ) of our setup. The reason lies in the fundamentally different action-at-a-distance mechanism between these problems. In the barotropic/baroclinic cases the induced velocity field generates vorticity at a distance by advecting the local mean (potential) vorticity gradient, so that the displacement and the vorticity anomalies are either in phase or anti-phased according to the sign of the mean vorticity gradient. In our problem, the induced velocity deforms the buoyancy in the far field, while vorticity is only subsequently induced in quadrature, by the new horizontal buoyancy gradients. This initial pure buoyancy anomaly necessarily excites both KGWs. Furthermore, in order to get growth the buoyancy (or displacement) fields cannot be fully in phase or in anti-phase with the vorticity field, hence they cannot be described solely by the counter-propagating kernel. To see this mathematically we write the displacement and vorticity in terms of their amplitudes and phases (ζ = Ze iθ , q = Qe iα ) and substitute a modal solution in the vorticity generation equation ( 2). Then by taking the ratio between the imaginary and the real parts of ( 2)

we obtain tan (θ -α) = c i c r -U (18)
In our setup tan (θ -α) = ∓c i for levels z = ±1, so that (everywhere outside the critical level) c i = 0 is obtained only when (θ -α) = (0, π). Since the pro-propagating waves must exist to obtain growth, but also must be phase locked, they must be hindered by the counter-propagating waves of the opposed level. The induced velocity by the latter on the former should be strong enough (hence the amplitude ratio χ) and in opposite signs (hence the anti-phased alignment) to overcome the pro-propagation tendency of the former (Fig. 6). Since for the unstable modes

(ζ + 1 , ζ - 2 ) = -χ(ζ + 2 , ζ - 1 )
, the four KGW equations ( 16) can be reduced (without approximation) into two equations for the two counter-propagating KGWs:

∂ ∂t + ik(-c gr + 1) ζ - 1 = -iσζ + 2 (19a) ∂ ∂t + ik(c gr -1) ζ + 2 = iσζ - 1 (19b)
where cgr = c gr -χ σ k , is the counter-propagating KGW intrinsic phase speed, reduced by the hindering effect of the anti-phased pro-propagating wave at the opposed boundary.

Writing ζ - 1 = Z - 1 e iǫ - 1 , ζ + 2 = Z + 2 e iǫ + 2
and substituting it in equation set ( 19), the real parts yield the instantaneous KGW growth rates:

Ż- 1 Z - 1 = Z + 2 Z - 1 σ sin (∆ǫ), Ż+ 2 Z + 2 = Z - 1 Z + 2 σ sin (∆ǫ) (20a,b)
while the imaginary parts yield the instantaneous KGW phase speeds:

c - 1 = - ǫ- 1 k = (1 -cgr ) + Z + 2 Z - 1 σ k cos ∆ǫ, c + 2 = - ǫ+ 2 k = -(1 -cgr ) - Z - 1 Z + 2 σ k cos ∆ǫ (21a,b) where ∆ǫ = ǫ + 2 -ǫ - 1 . For modal growth rate, kc i = Ż- 1 Z - 1 = Ż+ 2 Z + 2
, thus it is obvious that the KGWs have equal amplitudes so that kc i = σ sin (∆ǫ). In order for the waves to be phase locked, c r = c - 1 = c + 2 , the waves should be stationary with phase locking ∆ǫ, satisfying

cos (∆ǫ) = cgr-1 σ k
. Thus, the KGW dynamics is just like the two counter-propagating Rossby waves in barotropic and baroclinic instability, except that the intrinsic phase speed is proportional to the square root of the wavelength (rather than the wavelength itself for the Rossby kernels) and is somewhat hindered by an opposed anti-phased kernel.

For small wavenumbers the single counter-propagation phase speed is large (eq. 8) and the KGWs need to be anti-phased in order to hinder each other's counter propagation (Fig. 4c). Furthermore, the opposed pro-propagating kernel amplitude is large as well (Fig. 4b)

to amplify the hindering. Since the counter-propagation speed is also proportional to Ri, as it increases, such hindering can maintain phase locking only for larger wavenumbers, hence a long-wave cutoff is established for Ri > 1 (Fig. 4a). As the wavenumber slightly increases the additional hindering effect by the opposed anti-phased kernel allows the KGWs to be in a hindering-growing configuration (0 < ∆ǫ < π/2). Ignoring this effect (by taking χ = 0) prevents modal growth in this range (first column of Fig. 5). For larger wavenumbers the amount of hindering required for phase locking is reduced and the KGWs can be shifted towards a more amplifying configuration, which becomes optimal when ∆ǫ = π/2 (Fig. 6). This occurs for wavenumbers and Richardson numbers that satisfy cgr = 1, which is when each KGW (with its anti-phased opposed one), counterpropagates with a speed which balances exactly the mean flow speed. For even larger wavenumbers the counter-propagation speed is overwhelmed by the mean flow and the KGWs need to be in a helping-growing configuration (π/2 < ∆ǫ < π). In this range the hindering effect of the pro-propagating waves interferes with the phase locking and therefore χ decreases (Fig. 4b). Eventually the wavenumber becomes too large for effective phase locking (both because the counter-propagation speed and the interaction coefficient σ decrease with the wavenumber) and a shortwave cutoff is established. Increasing the

Richardson number increases the counter-propagation speed and hence shifts both the long-wave and shortwave cutoffs to higher wavenumbers but also shrinks the transition between them.

Summarizing this part of the analysis, we examine in detail Baines and Mitsudera's (1994) action-at-a-distance paradigm to stratified shear flow instability, in the simplest model we could construct. This model shares many similarities with the baroclinic Eady (1949) and barotropic Rayleigh (1880) models, but is fundamentally different in the sense that it does not satisfy the Rayleigh-Fjortoft conditions for instability, and that each interface supports two edge waves rather than one. Although the dynamics is more complex, the essence of the instability is the same -it is dominated by the phase locked constructive interaction between the two counter propagating edge waves.

These three models however, do not contain a critical level (they do contain a steering level in the middle, where U = c r , but the vorticity perturbation is zero there). Therefore, in the next session we proceed by adding a small ambient stratification to the two density jumps (Fig. 2b), to examine the effect of the CL on the instability from the vorticitybuoyancy action-at-a-distance perspective.

Connection with the Miles Howard Criteria

One motivation of this study, is to establish that the counter propagation of gravity waves explain the destabilisation of shear flows by stratification. Here destabilisation is in the sense that flows that do not satisfy the Rayleigh-Fjortoft conditions can be unstable when there are density vertical gradients. The example we treat in this section is to our knowledge the simplest configuration that explains this destabilization. The profiles we use are nevertheless discontinuous, whereas the M-H condition is related to continuously stratified flow. To circumvent this last difficulty, we can either adapt the M-H criteria to discontinuous profiles, or to find continuous density profiles that are unstable and that strongly resemble to our density jumps.

To adapt the M-H criteria to the discontinuous case, we can return to the M-H integral (Eq.4.3.3 in Baines 2004) and place into it our density profile. When doing so we the necessariy condition requires that

|φ z | 2 + k 2 |φ| 2 dz + Ri |φ| 2 |z -c| 2 z=-1 + |φ| 2 |z -c| 2 z=1 = 1 4 |φ| 2 |z -c| 2 dz , where φ = ψ √ (z-c)
.

EYAL: Here above I do not see exactly which conclusion you can get from this equation! But I put there what was a note after Eq. 12, in the previous version.

To find a continuous situation that resemble our density jumps, we will take for the density distribution

ρ = ρ 0 - |∆ρ| 2 tanh z -h ǫh + tanh z + h ǫh
where ǫ > 0 is a tunable parameter measuring the relative depth of the layers where the density varies compared to the distance between these layers. Hence, this continuous profile becomes equivalent to (10) in the limit ǫ → 0. In this profile the search for normal modes consist in solving the Taylor-Goldstein equation

(z -c) 2 ψ zz -k 2 ψ + Riψ = 0 where Ri = Ri 1 -tanh 2 ((z -1)/ǫ) 2ǫ + 1 -tanh 2 ((z + 1)/ǫ) 2ǫ .
To evaluate c we will next proceed numerically using a shooting method and use that for z >> 1 ψ ≈ e -kz . Starting for instance from z = 10, the Taylor Goldstein Eq. is integrated downward to z = -10 and using the Brankin et al. (1991) Runge-Kutta library subroutine. For a given k the values of c are then those for which the numerical integration gives ψ ≈ e kz within a very good approximation still at z = -10 (here the values z = ±10

for the start and the end points was found to be accurate enough, the results are not sensitive to an increase or to a moderate decrease of these bounds).

The results for the growth rates, of the unstable modes with c r = 0 and for various horizontal wavenumber are in Fig. 1a. Comparison with the exact results for the density jumps in the left column of Fig. 5, show excellent agreement. The Fig. 1b shows the vorticity distribution q for the most unstable mode (black line) together with the background gradient of density (green line). The comparison between these two curves shows that the vorticity anomalies stay confined to the layers where the density varies. The phase of q (red in Fig. 1b) shows that the the phase difference between the upper and lower is around 3π/2 (or -π/2) which means that the upper layer is almost in lead quadrature with the lower layer, as is the case for the corresponding mode for the density jumps in Fig. 4c (see also Fig. 8).

Suppression of instability by ambient stratification

We consider the density profile of Fig. 2b:

ρ = ρ 0 + |∆ρ|[1 -H(z -h) -H(z + h)] -z|ρ z | amb (22)
(|ρ z | amb = const) where the ambient Richardson number (outside the jumps) satisfies

0 < Ri amb = - g ρ 0 |ρ z | amb (U z ) 2 < 1 4 , so that Ri(z) = Ri amb + Ri[δ(z -1) + δ(z + 1)] (23)
We find the normal mode solutions of this setup both by numerical discretization of equations ( 2) and (4) (Appendix-A), and by converting them to the Taylor-Goldstein equation and solving the latter analytically (Appendix-B). The agreement between the two methods of solution is excellent (the two solutions for Ri = 2, Ri amb = 0.2 and k = 1.07 differ at most by 0.46%). As in the previous setup four discrete normal modes exist, for each zonal wavenumber, where the growing and decaying modes (if they exist) come in conjugate pairs.

Figure 7 shows the dispersion relation for the growing and decaying modes as Ri amb increases gradually from zero to 1/4. We choose to show the case where Ri = 2 since in the absence of ambient stratification the dispersion relation includes a distinct long wave cutoff and the role of the pro-propagating edge waves (hereafter EWs) is negligible, c.f. the upper row of Fig. 5 (the pro-propagating EWs were found to be unessential to the dynamics with ambient stratification). As Ri amb increases the growth rate decreases, and the spectrum emitting the instability shrinks and is shifted towards higher wavenumbers.

As expected from the M-H criterion the instability ceases when the Richardson number becomes larger than a quarter everywhere in the flow. In Fig. 8 we show the structure of the amplitudes and phases of the vorticity and the vertical displacements, for a specific wavenumber (k = 1.07), and for various values of Ri amb . The EWs hardly change with

Ri amb , the vorticity EWs remain tilted against the shear, by approximately a quarter of wavelength, and the displacement of the lower (upper) EW is almost in (anti) phase with its vorticity. Although the ambient stratification is constant everywhere in the domain the vorticity and the displacement are concentrated at the critical layer whose width decreases as the value of Ri amb approaches a quarter (e.g. Fig. 8).

This structure, which is drawn for Ri amb = 0.249 in Fig. 9, can be understood from the vorticity-buoyancy view when considering first the vertical displacement profile in the absence of ambient stratification. Although the vertical velocity resulting from the combined induction of the two EWs is minimized at the CL (due to the evanescence of the Green function), the vertical displacement is maximized there. This is because the fluid at the CL is phase locked with the two EWs and therefore its vertical displacement is amplified persistently (eq. ( 4) yields that ζ = w ik(U -cr)+kc i = w kc i at the CL). In the presence of ambient stratification the vertical displacement at the CL deforms the isopycnal surfaces and hence generates strong vorticity a quarter wavelength to its right in a standing wavelike fashion (q = i Ri amb c i ζ, eq. 2, Fig. 9, andc.f. Appendix-B in Harnik et al., 2008). Adding a small background stratification thus adds a third active region to the 2jump dynamics. It would simplify things, conceptually, to represent the CL region as a single standing-wave kernel with a structure at the exact center of the CL (q and ζ in quadrature both in space and time, c.f. appendix B of Harnik et al 2008). To verify that dynamically we can do that, we calculate the vertical velocity induced by the entire CL (the region where q and ζ are strongly phase-tilted with height), and compare it to the vertical velocity induced by a delta function vorticity at the CL center. While the delta function induces an untilted w field, the highly-structured CL region induces a w that is slightly phase tilted at the CL region, so that there is a small phase difference between the vertical velocity induced at the two jumps. However, this phase difference is tiny (0.012π for Ri amb = 0.225 and less than 0.0005π for Ri amb = 0.249), so that essentially the CL level region affects the flow like a single kernel at its center. The similar affect is due partly to the strong amplitude decay away from the CL center, and partly to the strong phase tilting which causes the top and bottom parts of the CL to largely cancel each other's influence at a distance.

Looking at a schematic illustration of a 3-kernel representation of the configuration of figure 9, drawn in figure 10, it is clear that if the induced vertical velocity field from the two kernels is shifted by a quarter wavelength, symmetry yields that the maximum induced displacement at the CL, should be one eighth wavelength in between the two edge waves. Then the vorticity anomaly at the CL, generated in quadrature in response to the induced displacement, induces a vertical velocity anomaly, which tends to oppose the existing displacement anomaly at all levels. Since this CL-induced vertical velocity is not exactly in anti-phase with the EW displacements, it also tends to shift the EW in the direction counter the shear (a hindering effect). Moreover, the CL-induced displacement anomaly also acts to decrease the existing CL vorticity anomaly, through the buoyancy restoring force.

As Ri amb approaches a quarter, the vorticity at the CL approaches infinity and its induced destructive velocity completely masks the constructive interference between the two EW vorticity delta functions 2 . The unstable mode which last survives is that for which the EWs are in the optimal constructive interference phase of π/2 (which is of wavenumber k = 1.07, shown in Figs. 8-10). As Ri amb crosses a quarter all modes become neutral. This schematic picture is slightly altered by the integrated effect of the non zero vorticity kernels outside the CL and the edges. These kernels are in phase with their adjacent edges (Fig. 8) and therefore act to help them to counter propagate against the shear. This effect overwhelms the opposed hindering by the CL and consequently the unstable regime is shifted toward the helping regime in the dispersion relation as Ri amb increases (Fig. 7). The other difference is that when Ri amb < 1/4 the CL is actually a layer with a finite width (Fig. 9). Nonetheless its overall effect remains destructive.

The destructive role of the CL seems surprising in light of the fact that this is the only place in the flow where energy is extracted from the mean flow (Lindzen and Tung, 1978).

For the barotropic/baroclinic case Harnik and Heifetz (2007) showed that the counter propagating Rossby wave kernels, adjacent to the CL, are governed by the action-at-adistance from remote dominating kernels. Next we show that this is also the case for the stably stratified shear setup. Lindzen and Tung (1978) and Lindzen and Barker (1985) looked at the energy flux carried by propagating waves (across the shear) at the vicinity of the CL and found that energy flux diverges at the CL when over-reflection occurs. To inspect the energy flux they applied the Eliassen-Palm (1960) first theorem (strictly valid for steady flow but they assumed a small growth rate):

pw = -ρ 0 (U -c)uw (24)
(where the bar denotes zonal averaging). For positive shear, energy flux divergence (pw > 0 above the CL, and pw < 0 below it) is obtained when the momentum flux uw < 0 at both its sides (and thus the energy growth, -uwU z , is positive there). This view of the CL and energy extraction from the mean flow arises naturally in the over-reflection view of things. Here we view this energy flux divergence at the CL using our kernel action-at-adistance approach to gain further insight into the relative contribution of the CL and EWs.

Toward this end we first develop an understanding of the the kernel action-at-a-distance contribution to the momentum flux uw.

By definition, the zonal and vertical velocities induced by a single kernel are in quadrature and therefore do not contribute to the momentum flux uw. Therefore, from the KGW perspective, the momentum flux is inherently non-local as it results from the interaction between different kernels. If we denote the contribution to the momentum flux at some level z, from two kernels located at z ′ and z ′′ , as [uw(z ′ , z ′′ )] z , then if z is sandwiched between the kernels, z ′′ > z > z ′ , figure 11 indicates that the momentum flux is negative (positive) if the kernels' vorticity is tilted westward (eastward) with height. Furthermore, the kernels' contributions are even in the sense that

[uw(z ′ , z ′′ )] z = [u(z ′ )w(z ′′ )] z + [u(z ′′ )w(z ′ )] z = 2[u(z ′ )w(z ′′ )] z . Moreover [uw(z ′ , z ′′ )] is in-
dependent of height everywhere between the two kernels (if z, for instance, is closer to the upper kernel from below the weaker contribution from the lower kernel is compensated by the stronger contribution from the upper one). In contrast, outside the kernels (z > z ′′ or

z < z ′ ) the contribution is vanished since [u(z ′ )w(z ′′ )] z = -[u(z ′′ )w(z ′ )] z .
Applying the cases of figure 11 to our problem, as presented for example in figures 8-10, it is clear that the interaction between the two EWs (tilted westward, against the shear) yields a negative contribution to uw in between them and hence a positive contribution to the energy growth. Thus, the mutual amplification between the two kernels can also be expressed by a positive contribution to the energy flux divergence at the CL. The kernel at the CL is tilted with the shear relative to both edge waves and therefore the interaction between the EWs and the CL contribute toward positive uw and hence toward energy flux convergence at the CL. Furthermore, since the vorticity at the critical layer (having a small finite width) is tilted with the shear it contributes toward a positive momentum flux.

In Appendix-C we derive the overall integrated expression for the momentum flux in terms of kernel interaction to be uw = -

z z ′ =-∞ ∞ z ′′ =z Q(z ′ )Q(z ′′ ) sin[α(z ′′ ) -α(z ′ )]e -k(z ′′ -z ′ ) dz ′ dz ′′ (25)
As a quantitative example we compute uw for the solution with k = 1.07, Ri amb = 0.24 (which is similar to that shown in Fig. 9), and calculate -uw at z = ±0.0308, which are the edges of the critical layer. Normalizing -uw to unity there we decompose the contribution into the upper-lower EW-EW interaction (which includes also the in phase ambient kernels in their vicinities); the overall upper and lower EWs interaction with the CL (EWs & CL); and the integrated kernels contribution within the critical layer (denoted as CL & CL). The budget becomes:

-uw = 1 = (EW & EW ) + (EW & CL) + (CL & CL) =

Summary and discussion

The dynamics of linear shear flow can be examined from a wave theory perspective or alternatively from a field theory one. These two perspectives are very different from each other.

The former views the dynamics in terms of the ability of waves to propagate/decay/become absorbed/become reflected or over-reflected, in different regions across the shear. The latter views the dynamics in terms of instantaneous action-at-a-distance interaction across the shear, between different remote vorticity waves propagating against the shear (but not across it). Each of the perspectives is valid, of course, and preferring one to the other is in some sense a matter of personal taste. et al. (2008), in an attempt to establish a coherent field theory perspective to stratified shear flow instability.

The need for such a perspective arises since the basic intuition associated with the Miles-Howard criterion -that the shear acts to destabilize the flow while the stratification acts to stabilize it -fails when considering a simple setup of Boussinesq plane Couette flow with two density jumps. Without the jumps the flow is sheared and unstratified, yet is neutral with respect to modal perturbations. Counter-intuitively, it is destabilized when stable stratification is added. This setup is analogous in some aspects to the baroclinic Eady model. There the over-reflection view (wave theory) seems a bit artificial, while the Rossby edge wave interaction view (field theory) seems more natural.

Opposed to the Eady (1949) model, where the Rossby EWs are generated by potential vorticity gradients, here the vorticity EWs are gravity waves generated by buoyancy gradients. In the Eady model each interface supports one counter-propagating Rossby wave, whereas here two gravity waves exist at each density jump, one counter propagating and one pro-propagating with the mean flow. Although this complicates the dynamics and although the wave generation mechanism is different, we showed in detail that the essence of the instability is the same as in the Eady model. The two counter propagating waves, one from each interface, dominate the dynamics as they phase lock each other in a growing configuration. The role of the pro-propagating waves was found to be minor (mainly in expanding the spectral range of instability in the low wavenumbers hindering regime).

These results are in line with the view of Baines and Mitsudera (1994).

This simple setup by itself emphasizes the subtleties associated with the M-H criterion.

A single density jump in a plane Boussinesq Couette flow satisfies the M-H necessary, but insufficient, condition for instability (since the Richardson number is a delta function there and zero elsewhere), nonetheless the flow is neutral. Only when another density jump is added, so that gravity waves can interact constructively in a distance, the flow becomes unstable. Furthermore, it seems strange that if we add an ambient stratification in between the density jumps, so that the ambient Richardson number exceeds the small value of a quarter (compared to infinity in the two jumps) the instability is shut down.

The explanation of Baines and Mitsudera (1994), that the waves are then absorbed by the critical level and hence are heavily damped, mixes arguments from the two perspectives.

Here we used the field-theory KGW formulation derived by Harnik et al (2008), which can be applied to continuous stratification as well, to show that the CL is the region which is persistently amplified by the edge waves, simply because fluid parcels there have no relative zonal motion with respect to the locked EWs (this is supported by initial value simulations that are not presented here due to lack of space). Hence, although the ambient stratification is small and uniform everywhere, at the CL the particle displacement anomalies are large, so that the resulting buoyancy generation of vorticity is also large.

The resulting continuous dynamics structure can be simplified to a three-kernel dynamics APPENDIX A

Numerical modal solution of the vorticity-displacement equations

Discretizing the normalized version of equations ( 2) and (4) yields:

( ∂ ∂t + ikz j )q j = -ikRi j ζ j (A1) ( ∂ ∂t + ikz j )ζ j = N n=1 q n G j,n (A2)
The domain is taken between (-3.5, 3.5) to assure no numerical boundary effect (the density jumps are located at z = ±1). N = 7 ∆z + 1 = 1751, is the number of points where sufficient accuracy is obtained for ∆z = 4 × 10 -3 . The Richardson number profile of ( 23) is discretized so that Ri(z = ±1) = Ri/∆z and Ri(z = ±1) = Ri amb , elsewhere. The discretized Green function is G j,n = -∆z • i 2 e -k|z j -zn| . Redefine G = i k G, (A1) and (A2) can be written then in the matrix form:

∂ ∂t             q 1 • • q N ζ 1 • • ζ N             2N ×1 = -ik •             z 1 0 • 0 Ri 1 0 • 0 0 • • 0 • • • • 0 • • 0 0 • 0 z N 0 • 0 Ri N G1,1 • • G1,N z 1 0 • 0 • • • 0 • • • • • • • 0 GN,1 • • GN,N 0 • 0 z N             A 2N ×2N             q 1 • • q N ζ 1 • • ζ N             2N ×1 (A3)
where the normal modes are the eigen-solution of (A3).

We hence decompose the setup of Fig. 2b into the three domains:

1. z > 1 ; 2. -1 < z < 1 ; 3. z < -1, and looking for the solution ψ j = √ y[A j H 1 (iy) + B j H 2 (iy)] where the subscript j=1,2,3 indicates the domain. We require continuity of the streamfunction:

z = 1 : ψ 1 = ψ 2 , & z = -1 : ψ 2 = ψ 3 (B6a,b)
Integration of the vorticity across z = 1 together with the definition (13) of the vorticity there yield:

∂ψ 1 ∂z - ∂ψ 2 ∂z z=1 = q1 (B7)
The T-G equation then suggests that

∂ψ 1 ∂z - ∂ψ 2 ∂z = - Ri (1 -c) 2 ψ z=1 (B8)
and similarly for z = -1 we obtain

∂ψ 2 ∂z - ∂ψ 3 ∂z = - Ri (1 -c) 2 ψ z=-1 (B9)
The outer radiation conditions of ψ 1 (z → ∞) = 0 and ψ 3 (z → -∞) = 0 set B 1 = A 3 = 0 Substitute (B5) in B(6,8,9) and write y 1 = k(1 -c) and y 2 = k(-1 -c) we get A 1 H 1 (iy 1 ) = A 2 H 1 (iy 1 ) + B 2 H 2 (iy 1 ) ; A 2 H 1 (iy 2 ) + B 2 H 2 (iy 2 ) = B 3 H 2 (iy 2 ) (B10a,b)

(A 1 -A 2 ) ∂ ∂z [ √ yH 1 (iy)] y 1 -B 2 ∂ ∂z [ √ yH 2 (iy)] y 1 -A 1 Ri (1 -c) 2 √ y 1 H 1 (iy 1 ) = 0 (B11) (B 2 -B 3 ) ∂ ∂z [ √ yH 2 (iy)] y 2 + A 2 ∂ ∂z [ √ yH 1 (iy)] y 2 -B 3 Ri (1 -c) 2 √ y 2 H 2 (iy 2 ) = 0 (B12)
Hence, in order to obtain a non trivial solution Det(M)=0, where This yields the phase speed c. The coefficients A 1 , A 2 , B 2 and B 3 are then obtained when solving the set of equations (B10,11). After deriving the solution for ψ we easily obtain the vertical displacement ζ = ψ z-c , and the vorticity q = -Ri (z-c) 2 ψ.

M =    H 1 (
then if z ′′ > z > z ′ , the integrand yields -1 8 Q(z ′ )Q(z ′′ ) sin[α(z ′′ ) -α(z ′ )]e -k(z ′′ -z ′ ) , which is independent of z and positive if the vorticity is tilted westward with height. Flipping between the dummy variables z ′ and z ′′ , yields the same contribution indicating that indeed [u(z ′ )w(z ′′ )] z = [u(z ′′ )w(z ′ )] z . If, however, z > z ′′ > z ′ or z < z ′ < z ′′ , then flipping between z ′ and z ′′ yields an exact cancellation and hence zero contribution. Taking these into account we obtain equation 25 as the explicit form of (C4). and for ǫ = 0.1. a) Growth rates kc i for the bulk Richardson numbers Ri = 0.5 (black), 1 (red), 1.5 (green), and 2 (red); b) vorticity distribution as a function of z of the unstable mode with k=0.43 when Ri =1: amplitude |q| (black), phase normalised so that it is nul in z = 1 (red); and density profile normalised by its value in z = 1(green dashed). .

  hence in the absence of interaction with the KGWs of the opposed interface (when the RHS are zeros), ζ + 1 and ζ - 2 tend to propagate in the direction of the mean flow whereas ζ - 1 and ζ + 2 propagate counter it. Each KGW, by construction, does not interact with the other KGW of the same interface, but with the two KGWs of the opposed one. The interaction coefficient σ ≡ e -2k 4 2k Ri, results from the vorticitydisplacement production relation (15), and the exponential evanescence of the induced vertical velocity as represented by the Green function, where the additional division by factor two results from the equi-partition of the vorticity at each interface between the two KGWs. The -i prefix on the RHS indicates that the vertical velocity is induced a quarter-wavelength to the right of the vorticity, recalling that for ζ + the vorticity and the displacement are in phase whereas for ζ -they are anti phased (hence the minus sign before ζ -, equation 15).

  For barotropic and baroclinic instabilities large bodies of literature exist for the two perspectives and recently Harnik and Heifetz (2007) related them, to obtain a deeper understanding of the instability mechanism. For stratified shear flow however, most of the research literature takes the wave theory point of view. This current study continues the work of Baines and Mitsudera (1994), and Harnik

  and Phase/π Ri=1, k=0.43, most unstable mode

Figure 1 :

 1 Figure1: Growth rate of the unstable in the continuous density profile given by (XX) and for ǫ = 0.1. a) Growth rates kc i for the bulk Richardson numbers Ri = 0.5 (black), 1 (red), 1.5 (green), and 2 (red); b) vorticity distribution as a function of z of the unstable mode with k=0.43 when Ri =1: amplitude |q| (black), phase normalised so that it is nul in z = 1 (red); and density profile normalised by its value in z = 1(green dashed). .

  EYAL: This is much better but still not quite clear. at which altitude this budget for uw apply. If it is for uwdz than it is wonderful, and should be really said. On top, say which altitudes z and z ′′ are taken to evaluate the separate contribution (EW & EW), ect... As this the corner stone of your discussion, it is mandatory to explain those thinks more clearly.

	78.08 -76.8 -0.28	(26)
	The (EW & EW) positive contribution is almost balanced by the (EW & CL) and indeed
	the mode is only slightly unstable.

  iy 1 ) -H 1 (iy 1 ) -H 2 (iy 1 ) 0 √ y 1 H 1 (iy 1 ) -∂ ∂z [ √ yH 1 (iy)]y 1 -∂ ∂z [ √ yH 2 (iy)]y 1 0 0 H 1 (iy 2 ) H 2 (iy 2 ) -H 2 (iy 2 )√ y 2 H 2 (iy 2 )

						
	∂ ∂z [ √ yH 1 (iy)]y 1 +	c Ri (1-c) 2				 
	0	∂ ∂z [ √ yH 1 (iy)]y 2	∂ ∂z [ √ yH 2 (iy)]y 2 -∂ ∂z [ √ yH 2 (iy)]y 2 +	c Ri (1+c) 2	(B13)

This setup is even simpler than the Holmboe model(1962), which is often used to address basic questions in stratified shear flows.

This argument is somewhat circular for normal modes since as the growth rate goes to zero, (2) indicates that the CL vorticity goes to infinity, therefore it is difficult to differ between cause and effect. For this reason we simulate the evolution towards the normal modes from various initial conditions, including beginning from the phase locked EWs with zero ambient vorticity. The simulations (not shown here due to lack of space) indicate indeed that the CL is the region which is persistently amplified by the edge waves, and as the CL vorticity grows it acts to decay the edge waves.

-two counter propagating EWs at the interfaces and one standing wave kernel at the CL.

From this perspective, the CL does not "absorb" the two interacting waves, rather it acts to induce a third kernel which decreases the two mutually amplifying kernels. This instability is allowed only when the CL effect is weak enough, that is when the Richardson number is smaller than a quarter. This destructive role of the CL may seem surprising, when recalling that the CL is the only place in the flow where energy flux is diverged, hence where energy is being yanked from the mean flow to the the perturbation. The kernel formulation provides a rigorous way to compute the action-at-a-distance contribution of the remote kernels to the Energy flux. It then turns out that it is the interaction between the two EWs that contributes toward positive energy flux divergence at the CL, while the interaction between the CL kernel and the EWs contributes toward energy flux convergence, therefore toward stabilization of the flow (a similar effect has been observed by Harnik and Heifetz in barotropic and baroclinic shear flows).

To conclude, we investigated one of the simplest setups for stably stratified shear flow instability and still the physical mechanism of the instability is far from being simple.

A possible reason for that is that the M-H theorem is necessarily spectral (one must assume a-priori a modal solution to obtain the condition), hence is not derivable from Hamiltonian considerations (Shepherd, 1990). This stands in contrast to barotropic and baroclinic shear flows whose instability conditions can be derived from globally conserved Hamiltonians (Heifetz et al., 2009).

For the transformation of variables:

the modified Bessel equation (e.g. Sneddon, 1956):

with ν = 1/4 -Ri, whose general solution is:

where I ν and I -ν are the modified Bessel functions and C 1 , C 2 are two constants. Equivalently we can obtain the Bessel function from T-G by writing y = iy:

whose solution is

where H 1 and H 2 are the Hankel functions and A, B are two constants. Since the two formulations are equivalent we arbitrary choose to use the Hankel functions for the matching solution.

APPENDIX C

Integrated kernel contribution to the momentu flux

We wish to express

in terms of the inversion of the vorticity profile q = Q(z)e i[kx+α(z)] . For clarity, we redefine the Green function to be the real, positive definite, and non-dimensional evanescent function; G(z, z ′ ) = iG(z, z ′ ) = 1 2 e -k|z-z ′ | , and use (3) to write

where the sine in the integrand simply indicates that the induced vertical velocity is located quarter of wavelength to the east of the inducing vorticity anomaly (Fig. 11).

Similarly we can write that

Since ∂ ∂z G(z, z ′′ ) > 0 for z < z ′′ , and ∂ ∂z G(z, z ′′ ) < 0 for z > z ′′ , the integrand indicates that the induced zonal velocity anomaly below (above) the inducing vorticity anomaly are in (anti) phase. Substitute (C2 and C3) in (C1) we obtain after zonal integration

It is clear from the integrand that a single kernel (where z ′ = z ′′ ) does not contribute to the momentum flux. Writing explicitly the Green function and its derivative G(z, z ′ ) = 1 2 e -k(z-z ′ ) for z > z ′ e k(z-z ′ ) for z < z ′ , ∂ ∂z G(z, z ′′ ) = k 2 -e -k(z-z ′′ ) for z > z ′′ e k(z-z ′′ ) for z < z ′′ (C5)


Ugo Boscain

Paolo Mason

Gianluca Panati

Mario Sigalotti

Controllability of spin-boson systems

come

1 Introduction and description of the models

The Rabi model alias the standard spin-boson model

In the Hilbert space H = L 2 (R, C) ⊗ C 2 , we consider the Schrödinger equation

i ∂ t ψ = H Rabi ψ, (1)
where

H Rabi = ω 2 (P 2 + X 2 ) ⊗ 1 + Ω 2 1 ⊗ σ 3 + g X ⊗ σ 1 , (2)
for some fixed positive constants ω and Ω, and g ∈ R, where X is the multiplication operator defined by (Xψ)(x) = xψ(x), and P = -i∂ x . Here and in the following we use σ 1 , σ 2 , σ 3 to denote the Pauli matrices

σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 , σ 3 = 1 0 0 -1 .
Using the natural isomorphism H ∼ = L 2 (R, C 2 ), one rewrites (1) as the system

i ∂ t ψ + = ω 2 (-∂ 2 x + x 2 )ψ + + Ω 2 ψ + + gxψ - i ∂ t ψ -= ω 2 (-∂ 2 x + x 2 )ψ -- Ω 2 ψ -+ gxψ + .
In the physics literature, it is customary to introduce the raising and lowering matrices

σ + = 1 2 (σ 1 + i σ 2 ) = 0 1 0 0 , σ -= 1 2 (σ 1 -i σ 2 ) = 0 0 1 0 . 1 Setting a = 1 √ 2 (X + i P ), a † = 1 √ 2 (X -i P ),
and omitting tensors, H Rabi appears in physics literature as

H Rabi = ω a † a + 1 2 + Ω 2 σ 3 + g √ 2 (a + a † )(σ -+ σ + ). (3)
Modeling of spin-boson interactions was initiated by Rabi in the 30's [START_REF] Rabi | On the process of space quantization[END_REF][START_REF] Rabi | Space quantization in a gyrating magnetic field[END_REF]. To our knowledge, the Hamiltonian (3) was first derived from a more fundamental model in the milestone paper [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF], where also the simpler Jaynes-Cummings Hamiltonian appears (see Section 1.2). A throughout discussion on various aspects of this model can be found in [START_REF] Allen | Optical Resonance and Two-level Atoms[END_REF]. Recent results on its integrability have been obtained in [START_REF] Braak | Integrability of the Rabi model[END_REF].

The Hamiltonian [START_REF] Allen | Optical Resonance and Two-level Atoms[END_REF] is ubiquitous in the literature, since it encodes one of the simplest possible couplings between a harmonic oscillator and a two-level system. The physical interpretation of the two factors appearing in the tensor product varies according to the context. We briefly outline the interpretation in cavity QED in the following digression.

Digression: the Rabi model in the context of cavity QED

While we refer to the specific literature for an exhaustive treatement of the problem [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF][START_REF] Cohen-Tannoudji | Photons & Atoms: Introduction to Quantum Electrodynamics[END_REF][START_REF] Raimond | Colloquium: Manipulating quantum entanglement with atoms and photons in a cavity[END_REF], we briefly sketch the derivation of the Hamiltonian (3) in the context of Cavity Quantum Electro Dynamics (CQED). In a typical CQED experiment, an atom moves across a Fabry-Perot cavity interacting with the quantized electromagnetic (EM) field of the cavity. The atom is initially prepared in a special state, to guarantee both a strong interaction with the cavity field and a decay time (to the atomic ground state) longer than the experiment timescale. The first goal is obtained by choosing a hydrogenoid atom (e.g., a Rubidium atom) prepared in a circular Rydberg state with high orbital angular momentum, and so with a large dipole moment. The second goal requires instead that transitions to the other atomic states are not resonant with the radiation from the environment, typically in the visible or infrared spectrum. Since the atomic energies are E n = -R Rb n 2 (where R Rb is the Rydberg constant, corrected for the Rubidium), one select n ∈ N in such a way that the transition frequency Ω at := 1 (E n -E n-1 ) is in the range of the microwaves (in the case of Rb, this is obtained by choosing n = 51). The cavity is designed in such a way that a normal mode of the EM field has frequency ω, which is almost resonant with Ω at , i.e., |Ω at -ω| ω, Ω at . Under these conditions, the only appreciable interactions are those among the two eigenstates |n and |n -1 , and the distinguished mode of the EM field with frequency ω. The linear space generated by |n and |n -1 corresponds to the factor C 2 in H, while the quantized EM mode yields a quantum harmonic oscillator, corresponding to the factor L 2 (R, C). This explains the structure of the first two terms in [START_REF] Allen | Optical Resonance and Two-level Atoms[END_REF].

A deeper understanding of (3) is obtained by deriving it, heuristically, from the Pauli-Fierz model, see [START_REF] Spohn | Dynamics of Charged Particles and their Radiation Field[END_REF]Chapter 13]. We consider a hydrogenoid atom, as, e.g., Rubidium, interacting with the quantized electromagnetic (EM) field in the cavity. The electron, with effective mass m * and charge -e, is described by a form factor ϕ ∈ C ∞ 0 (R 3 ), which is assumed to be a radial positive function, normalized so that R 3 ϕ(x)dx = 1. The atomic core is supposed at rest at the origin or, equivalently, moving across the cavity with constant velocity.

Within the Pauli-Fierz model, the quantized EM fields are described by operator-valued distributions acting on the Hilbert space H at ⊗ H f , where H at = L 2 (R 3 , dr) corresponds to the atom, while

H f = n∈N L 2 (R 3 , dk) ⊗ C 2 ⊗n sym
is the bosonic Fock space corresponding to the field [START_REF] Spohn | Dynamics of Charged Particles and their Radiation Field[END_REF]. In the Coulomb gauge, the vector potential (in free space) reads

A ϕ (r) = 2 λ=1 R 3 dk c 2 2ω(k) ϕ(k) e λ (k) e ik•r ⊗ a(k, λ) + e -ik•r ⊗ a(k, λ) † , (4)
where ω(k) = c|k|, ϕ denotes the Fourier transform of ϕ, the three vectors {k/|k|, e 1 (k), e 2 (k)} form an orthonormal basis at every point k ∈ R 3 , and a(k, λ) and a(k, λ) † satisfy the canonical commutations relations. Analogously, the transverse electric field reads

E ϕ,⊥ (r) = 2 λ=1 R 3 dk ω(k) 2 ϕ(k) e λ (k) i e ik•r ⊗ a(k, λ) -e -ik•r ⊗ a(k, λ) † . (5)
Actually, since the field is confined in the cavity, described by a compact set C ⊂ R 3 , the integral appearing in ( 4) and ( 5) can be replaced by an infinite sum k∈I , where I labels the solutions to the eigenvalue equation ∆A = λ k A satisfying ∇•A = 0 and appropriate boundary conditions on ∂C, see [START_REF] Spohn | Dynamics of Charged Particles and their Radiation Field[END_REF]Sec. 13.6]. The (spinless) Pauli-Fierz Hamiltonian is the self-adjoint operator

H = 1 2m * i ∇ r + e c A ϕ (r) 2 + 1 ⊗ H f -e 2 V ϕ,Coul (r) ⊗ 1, (6)
where

H f = 2 λ=1 R 3 dk ω(k) a(k, λ) † a(k, λ)
corresponds to the field energy, while

V ϕ,Coul (r) = R 3 dk | ϕ(k)| 2 |k| 2 e -ik•r = R 3 dx R 3 dy ϕ(x -r)ϕ(y) 4π|x -y|
is the Coulomb interaction smeared through ϕ.

To derive the Rabi Hamiltonian (3) from ( 6), several approximations are in order. (1) First, since the electron is bound to the atomic core, "one loses little by evaluating the vector potential at the origin instead of at r, the position of the electron " [START_REF] Spohn | Dynamics of Charged Particles and their Radiation Field[END_REF]Sec. 13.7]. Within such approximation, and after a suitable change of gauge, the Hamiltonian becomes

H dip = - 2 2m * ∆ r + V ϕ,Coul (r) ⊗ 1 + 1 ⊗ H f -e r • E ϕ,⊥ (0) + dip (r), (7)
( 1) Obviously, in the comparison, the order of factors in the tensor product should be reversed, so that L 2 (R, C) ⊗ C 2 is an effective reduced space for H f ⊗ H at . Since both conventions are well-established, we did not dare to reverse the order of factors in any of the two models, namely in (2) and [START_REF] Boscain | Adiabatic control of the Schrödinger equation via conical intersections of the eigenvalues[END_REF]. j=1 r j ⊗ E j (0). The expression for H dip involves the dipole operator D = -e r, thus explaining the name dipole approximation. Here a quadratic potential appears, namely

dip (r) = 1 2 2 3 R 3 dk e 2 | ϕ(k)| 2 r 2 ,
which is usually neglected in this context. As a second crucial approximation, only a distinguished mode of the cavity field is considered as interacting with the atom. If k * is its wavevector, the integrals (or sums) above reduce to the evaluation at k * . Accordingly, and assuming that ϕ(k * ) = 1, one has

A ϕ (0) = 2 λ=1 c 2 2ω(k * ) e λ (k * ) 1 ⊗ a(k * , λ) + a(k * , λ) † E ϕ,⊥ (0) = 2 λ=1 ω(k * ) 2 e λ (k * ) 1 ⊗ i a(k * , λ) -a(k * , λ) † . H f = 2 λ=1 ω(k * ) a(k * , λ) † a(k * , λ). (8)
Neglecting the sum over the polarization vectors e λ 's, we select an index λ * ∈ {1, 2} and we set ω(k * ) ≡ ω and a(k * , λ * ) ≡ a c , where the subscript c stands for cavity. Then one has the correspondences

A ϕ (0) a c + a † c = √ 2X c , and
E ϕ,⊥ (0) i(a c -a † c ) = - √ 2P c , (9)
where the last equality refers to the fact that the operators

1 √ 2 (a c + a † c ) and -i √ 2 (a c -a † c
) can be identified, up to unitary equivalence, with a canonical Schrödinger pair (X c , P c ) acting on L 2 (R, dq). Next, an inverse Fourier transform F -1 : L 2 (R, dq) → L 2 (R, dx) is performed, which intertwines the pair (X c , -P c ) with the pair (P, X), where P = -i∂ x and (Xψ)(x) = xψ(x). Correspondingly, F -1 a c F = -ia, so that H f is unitarily transformed into Hf = ω c a † a, which corresponds to the first addendum in (3), up to an irrelevant constant. Moreover,

A ϕ (0) -i(a -a † ) = √ 2P , and E ϕ,⊥ (0) a + a † = √ 2X, (10)
which shows that the interaction term (the third addendum in [START_REF] Albert | Genericity of simple eigenvalues for elliptic PDE's[END_REF] or in [START_REF] Allen | Optical Resonance and Two-level Atoms[END_REF]) is proportional to the electric field, which is assumed to be uniform in view of the dipole approximation. Third, one assumes that effectively only two atomic states |n -1 ≡ |g and |n ≡ |e are appreciably coupled to the radiation field. Then, the factor H at is replaced by a 2-dimensional Hilbert space

H 2L = Span C {|g , |e } C 2 ,
and the atomic Hamiltonian (the first term in [START_REF] Braak | Integrability of the Rabi model[END_REF]) by the 2-level Hamiltonian Ωat 2 σ 3 , where

Ω at = E n -E n-1 .
Finally, if the two selected atomic states are circular Rydberg states one has that the dipole operator is off-diagonal in the basis {|g , |e }, namely

g| D |g = 0, e| D |e = 0.
Moreover, there exist d > 0 and a versor e at ∈ C 3 such that d e at = e| D |g . We identify a state ψ at ∈ H 2L with the pair (ψ g , ψ e ) ∈ C 2 via ψ at = ψ g |g + ψ e |e , so that Dψ at = d (ψ g e * at |e + ψ e e at |g ) .

By using the latter fact, one obtains

D • E ϕ,⊥ (0) = ω(k * ) 2 d e λ * (k * ) • e * at σ + ⊗ i(a c -a † c ) + e at σ -⊗ i(a c -a † c ) , (11)
where, with a little abuse of notation, we wrote σ -= |g e| and σ + = |e g|. When adapting the previous formula to concrete physical situations, by specializing to the actual values of e at and e λ (k * ), one obtains some terms analogous to the ones appearing in (3), possibly with different weights. In particular, if e * at = e at , one obtains a coupling term proportional to σ 1 ⊗ P c , which, after the Fourier transform mentioned above, becomes σ 1 ⊗ X and so agrees with the third addendum in [START_REF] Albert | Genericity of simple eigenvalues for elliptic PDE's[END_REF] or in [START_REF] Allen | Optical Resonance and Two-level Atoms[END_REF].

This concludes the heuristic "derivation " of the Rabi Hamiltonian, starting from a fundamental model as the one given by ( 6). On the other hand, a rigorous derivation of the Hamiltonian (3) from a mathematical model of non-relativistic QED, as, e.g., the Pauli-Fierz model, in a suitable scaling limit, is to our knowledge absent from the literature.

The Jaynes-Cummings model

The Rabi model presented in the previous section has been considered difficult to treat in view of the widespread opinion that it is not explicitly solvable (cf. [START_REF] Braak | Integrability of the Rabi model[END_REF]).

For this reason Jaynes and Cummings [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF] proposed an approximated Hamiltonian, which is obtained from H Rabi by neglecting the so-called counter-rotating terms. More precisely, the interaction term in H Rabi , namely

g √ 2 (a † + a)σ 1 = g √ 2 (a † σ -+ aσ + + a † σ + + aσ -),
is replaced by a new interaction term where the last two terms in the right-hand side are neglected under the assumption that

|Ω -ω| Ω, g Ω, ω. (12)
According to this procedure, one obtains

H JC = ω a † a + 1 2 + Ω 2 σ 3 + g √ 2 (a † σ -+ aσ + ) (13) = ω 2 (P 2 + X 2 ) ⊗ 1 + 1 ⊗ Ω 2 σ 3 + g 2 (X ⊗ σ 1 -P ⊗ σ 2 ).
This model is considerably simpler than the Rabi one, since it admits the conserved quantity

C = a † a ⊗ 1 + 1 ⊗ 1 0 0 0 ,
which represents the total number of excitations. Indeed, the interaction term in H JC contains two parts: the term a † σ -creates an excitation of the bosonic mode and destroys one of the two-level system, while aσ + acts in the opposite way.

The heuristic justification of the approximation leading to ( 13) is based on separation of time scales. Indeed, when considering the dynamics in the interaction picture, the terms a † σ - and aσ + , which conserve the total number of excitations of the system, evolve periodically with frequency |Ω -ω|, while the remaining terms a † σ + and aσ -oscillate with frequency Ω + ω. If the almost-resonance condition (12) is satisfied, the latter oscillations average to zero on the longer time scale |Ω -ω| -1 . On the other hand, a rigorous mathematical justification for this approximation seems absent from the literature. However, as already noticed by Rouchon in [START_REF] Rouchon | Quantum systems and control[END_REF], it seems possible to obtain it by adapting the methods developed in [START_REF] Panati | Space-adiabatic perturbation theory in quantum dynamics[END_REF][START_REF] Panati | Space-Adiabatic Perturbation Theory[END_REF][START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF]. A similar task, by analytical methods, has been accomplished by Ervedoza and Puel on a related model [START_REF] Ervedoza | Approximate controllability for a system of Schrödinger equations modeling a single trapped ion[END_REF].

Notice that in applications to circuits QED the assumption g Ω, ω is in general not satisfied, since it is possible to achieve spin-filed interactions which are not much smaller than the field energy [START_REF] Schuster | Resolving photon-number states in a superconducting circuit[END_REF]. This is a motivation for the direct study of the Rabi model, which is pursued here.

Modeling an external control field

In a wide variety of experimental situations one can act on the system by an external field. The goal of the controller could be to lead the system from a given initial state to a prescribed final one. For spin-boson models, this amounts to study the control problem

i ∂ t ψ(t) = H 0 ψ(t) + H c (u(t))ψ(t), (14)
where H 0 represents either the Rabi or the Jaynes-Cummings Hamiltonian and H c is a selfadjoint operator describing the coupling between the system and the controlled external field.

The operator H c depends on the control function u which, in general, takes values in C m . In most cases the external field can act on the bosonic mode only, while the spin mode is not directly accessible. This leads to a control Hamiltonian of the form

H c (u(t)) = h c (u(t)) ⊗ 1. ( 15
)
where h c (u(t)) is a self-adjoint operator in L 2 (R). One of the simplest form for the operator

h c (u(t)) is the following h c (u(t)) = u(t)X. (16)
As argued in Section 1.1.1, in the context of cavity QED and within the dipole approximation, the operator X corresponds to the electric field at the origin, see [START_REF] Brouder | Adiabatic approximation, Gell-Mann and Low theorem, and degeneracies: A pedagogical example[END_REF]. Therefore, in this context, the control ( 16) corresponds to an external EM field which rescales the value of E ϕ,⊥ (0) by a time-dependent factor u(t). Notice that the linearity in X is a consequence of the dipole approximation.

Simmetrically, one might also consider a control consisting of an external field which modulates the value of A ϕ (0), as e.g., an externally-controlled magnetic field. In view of [START_REF] Brouder | Adiabatic approximation, Gell-Mann and Low theorem, and degeneracies: A pedagogical example[END_REF], this amounts to consider a control term h c (u(t)) = u(t)P , a problem which will be investigated elsewhere.

Other related controlled models

Controllability of finite-dimensional approximations of spin-boson systems have first been obtained in the physics literature, via constructive methods, in [START_REF] Law | Arbitrary control of a quantum electromagnetic field[END_REF] and [START_REF] Kneer | Preparation of arbitrary entangled quantum states of a trapped ion[END_REF]. The results have been extended to the full infinite-dimensional model in [START_REF] Yuan | Controllability of the coupled spin-1 2 harmonic oscillator system[END_REF] via analysis of non-resonances and in [START_REF] Adami | Controllability of the Schrödinger equation via intersection of eigenvalues[END_REF] by adiabatic arguments (see also [START_REF] Boscain | Adiabatic control of the Schrödinger equation via conical intersections of the eigenvalues[END_REF]).

A mathematical description of the controllability properties of a related model is given in [START_REF] Ervedoza | Approximate controllability for a system of Schrödinger equations modeling a single trapped ion[END_REF]. The main technical tool is an explicit estimate of the approximation error with respect to the Lamb-Dicke limit.

For the counterpart of the Eberly-Law model for more than one trapped ion an approximate controllability result is obtained in [START_REF] Bloch | Finite Controllability of Infinite-Dimensional Quantum Systems[END_REF] (see also [START_REF] Rangan | Control of trapped-ion quantum states with optical pulses[END_REF]), based on the analysis on a sequence of nested finite-dimensional systems. In the recent papers [START_REF] Keyl | Controlling Several Atoms in a Cavity[END_REF][START_REF] Paduro | Approximate controllability of the two trapped ions system[END_REF] the (approximate) controllability of the system is established by considering different families of controlled dynamics.

Approximate controllability of the controlled Rabi model

We consider here the approximate controllability problem for a system of the form [START_REF] Kato | Perturbation theory for linear operators[END_REF], where H 0 is the Rabi Hamiltonian and H c takes the form ( 15)-( 16).

The problem under consideration is then

       i∂ t ψ = H Rabi ψ + u(X ⊗ 1)ψ, u ∈ [0, δ], |Ω -ω| Ω, δ, Ω > 0. (17)
The main result of the paper is the following. The precise definition of approximate controllability will be given in the next section. It basically means that, for every choice of the initial and final state, there exists an admissible control law u, depending on the time, which steers the initial state arbitrarily close to the final one.

Theorem 1. System (17) is approximately controllable for almost every g ∈ R.

The proof of Theorem 1 goes by studying the applicability of a general approximate controllability result in dependence on the parameter g. In order to do so, we have to use perturbation theory in the parameter g up to order 4.

An approximate controllability result

We are going to recall a general controllability result for bilinear quantum systems in an abstract setting.

In a separable Hilbert space H, endowed with the Hermitian product •, • , we consider the following control system

d dt ψ = (A + u(t)B)ψ, u(t) ∈ U, (18)
where (A, B, U ) satisfies the following assumption. Assumption 1. U is a subset of R and (A, B) is a pair of (possibly unbounded) linear operators in H such that 1. A is skew-adjoint on its domain D(A);

2. there exists a Hilbert basis (φ k ) k∈N of H consisting of eigenvectors of A: for every k, Aφ k = iλ k φ k with λ k in R;

3. for every j in N, φ j is in the domain D(B) of B;

4. A + uB is essentially skew-adjoint for every u ∈ U ;

5. Bφ j , φ k = 0 for every j, k in N such that λ j = λ k and j = k.

If (A, B, U ) satisfies Assumption 1, then A + uB generates a unitary group t → e t(A+uB) . By concatenation, one can define the solution of ( 18) for every piecewise constant function u taking values in U , for every initial condition ψ 0 given at time t 0 . We denote this solution by t → Υ u t,t 0 ψ 0 . For every j, k ∈ N, we denote b jk = φ j , Bφ k . A pair (j, k) in N 2 is a non-resonant transition of (A, B) if b jk = 0 and, for every l, m, |λ j -

λ k | = |λ l -λ m | implies {j, k} = {l, m} or {l, m} ∩ {j, k} = ∅.
A subset S of N 2 is a chain of connectedness of (A, B) if for every j, k in N, there exists a finite sequence p 1 = j, p 2 , . . . , p r = k for which (p l , p l+1 ) ∈ S for every l and φ p l+1 , Bφ p l = 0 for every l = 1, . . . , r -1. A chain of connectedness S of (A, B) is non-resonant if every (j, k) in S is a non-resonant transition of (A, B). Definition 2. Let (A, B, U ) satisfy Assumption 1. We say that ( 18) is approximately controllable if for every ε > 0, for every ψ 0 , ψ 1 ∈ H, there exists a piecewise constant function

u ε : [0, T ε ] → U such that Υ u Tε,0 ψ 0 -ψ 1 < ε. Theorem 3 ([5]
). Assume that [0, δ] ⊂ U for some δ > 0 and let (A, B, U ) satisfy Assumption 1 and admit a non-resonant chain of connectedness. Then system (18) is approximately controllable.

Proof of Theorem 1

In this section, we consider the approximate controllability problem for (17) and we prove Theorem 1. The proof of the theorem is based on a suitable application of Theorem 3.

The strategy of the proof is the following. We first show in Subsection 4.1.1 that, for almost every g in R, some relevant pairs of eigenvalues of H Rabi satisfy the non-resonance condition, see [START_REF] Rabi | Space quantization in a gyrating magnetic field[END_REF]. This goal is achieved by exploiting the analyticity of the eigenvalues and by using perturbation theory. Then, in Subsection 4.1.2, we prove that these pairs of eigenvalues correspond to non-resonant transitions, according to the definition above. The degenerate case ω = Ω is treated separately in Section 4.2 Preliminarily, we introduce some additional notations. Denote by H Rabi,0 the Hamiltonian H Rabi where we set g = 0. The general form of H Rabi is then

H Rabi = H Rabi,0 + gV, ( 19
)
where V is the operator X ⊗ σ 1 . Let (ϕ j ) j∈N be the standard Hilbert basis of L 2 (R, C) given by real eigenfunctions of -∂ 2

x + x 2 , so that (-∂ 2 x + x 2 )ϕ j = (2j + 1)ϕ j and R xϕ j (x)ϕ j+1 (x)dx = (j + 1)/2 for j ≥ 0. Let ν 1 = (1, 0) T , ν -1 = (0, 1) T be the canonical orthonormal basis of C 2 . Based on (ϕ j ) j∈N , we obtain a Hilbert basis of factorized eigenstates Φ j,s = ϕ j ⊗ν s , j ∈ N, s ∈ {-1, 1}, of H Rabi,0 whose corresponding eigenvalues are

E j,s = ω j + 1 2 + s 2 Ω. ( 20
)
If Ω is not an integer multiple of ω then each eigenvalue E j,s is simple. (See Figure 2.) In the following, for ease of notations, we write in bold the elements of N × {-1, 1}, and for every j ∈ N × {-1, 1} we define n(j), s(j) in such a way that j = (n(j), s(j)).

For g ∈ R, denote by E g j , j ∈ N × {-1, 1}, the eigenvalues of H Rabi repeated according to their multiplicities, and by Φ g j , j ∈ N × {-1, 1}, an orthonormal basis of corresponding eigenstates. By a consequence of perturbation theory, detailed in the Appendix, we can assume that, up to the choice of a suitable labelling, each map

R g → (E g j , Φ g j ) ∈ R × L 2 (R, C) is analytic, with E 0 j = E j
, where E j is defined in [START_REF] Panati | Space-adiabatic perturbation theory in quantum dynamics[END_REF]. When E j is simple we can assume, in addition, that Φ 0 j = Φ j . Under the assumption |Ω -ω| Ω, the only case in which H Rabi,0 has multiple eigenvalues is when ω = Ω. If this is the case all eigenvalues different from the lowest one are double. Equation (32) allows to identify the splitting of the two-dimensional eigenspaces induced by the perturbation g → H Rabi,0 + gV , as detailed in Section 4.2.

In order to study the first and higher-order derivatives of g → E g j,s at g = 0, it is useful to introduce the quantities

V i,j = Φ i , (x ⊗ σ 1 )Φ j = δ n(i),n(j)-1 n(j) 2 + δ n(i),n(j)+1 n(j) + 1 2 (1 -δ s(i),s(j) ). (21)
4.1 Case ω = Ω Let us first prove that for almost every g ∈ R and every i, j, k, l ∈ N × {-1, 1}, with (i, j) = (k, l) and i = j, one has

E g i -E g j = E g k -E g l .
In order to do so, we observe that it is enough to show that for fixed i, j, k, l ∈ N × {-1, 1} as before, the set

S i,j,k,l = {g | E g i -E g j = E g k -E g l } (22)
is of full measure. By the analytic dependence on g of the eigenvalues of H Rabi , this is equivalent to say that g → E g i -E g j and g → E g k -E g l have different Taylor expansions at g = 0.

Let us consider the Taylor expansion

E g j = E j + ∞ m=1 g m E (m) j
.

The computation of the coefficients E (m) j carried on below is based on the Rayleigh-Schrödinger series (see, for instance, [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]Chapter XII]).

First of all we observe that E i -

E j = E k -E l is equivalent to n(i) -n(j) = n(k) -n(l) and s(i) -s(j) = s(k) -s(l), in the case in which Ω is not an integer multiple of ω/2. If Ω = (2m + 1)ω/2 for some non-negative integer m, then E i -E j = E k -E l implies n(i) + n(l) -n(j) -n(k) = 2m + 1 4 s(j) + s(k) -s(i) -s(l) ,
and thus, if the left-hand side is an integer number different from zero, it must be |s(j

) + s(k) -s(i) -s(l)| = 4, that is s(j) = s(k) = -s(i) = -s(l).
The term E

(1) j coincides with V j,j = Φ j , (X ⊗ σ 1 )Φ j , thus we deduce from ( 21) that

E

(1) j = 0 for every j.

Following [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF] we have that

E (2) j = - m =j (E m -E j ) -1 V j,m V m,j . Thus E (2) j = - m =j (E m -E j ) -1 (1 -δ s(j),s(m) ) 2 δ n(j),n(m)-1 n(j) + 1 2 + δ n(j),n(m)+1 n(j) 2 2 = -(E n(j)+1,-s(j) -E j ) -1 n(j) + 1 2 -(E n(j)-1,-s(j) -E j ) -1 n(j) 2 = -(ω -s(j)Ω) -1 n(j) + 1 2 + (ω + s(j)Ω) -1 n(j) 2 = ω + s(j)Ω(2n(j) + 1) 2(Ω 2 -ω 2 ) .
Notice that the computation above is correct also for n(j) = 0, even if in this case E n(j)-1,-s(j) is not a well defined eigenvalue of H Rabi,0 . Indeed, in this case the term

(E n(j)-1,-s(j) -E j ) -1 n(j) 2
counts as zero. Let us identify the values i, j, k, l such that E

(2)

i -E

(2) j = E

(2)

k -E

(2) l under the assumption that E i -E j = E k -E l . Recall that we also assume that (i, j) = (k, l) and i = j. From the above expression of E

(2) j we have

s(i) (2n(i) + 1) -s(j) (2n(j) + 1) = s(k) (2n(k) + 1) -s(l) (2n(l) + 1) . ( 23
)
If Ω = (2m + 1)ω/2 for some non-negative integer m and s(j) = s(k) = -s(i) = -s(l) then [START_REF] Raimond | Colloquium: Manipulating quantum entanglement with atoms and photons in a cavity[END_REF] gives n(i)+n(j)+n(k)+n(l)+2 = 0, which is impossible being the addends non-negative. The remaining case is when n(i) -n(j) = n(k) -n(l) and s(i) -s(j) = s(k) -s(l), in which case

s(i)n(i) -s(j)n(j) = s(k)n(k) -s(l)n(l). ( 24
)
Then, either s(i) = s(j), which implies s(k) = s(l) and then, by [START_REF] Rangan | Control of trapped-ion quantum states with optical pulses[END_REF],

s(i) = s(j) = s(k) = s(l),
or s(i) = -s(j), which implies s(k) = -s(l) and then, by [START_REF] Rangan | Control of trapped-ion quantum states with optical pulses[END_REF], n(i) + n(j) = n(k) + n(l). In the latter case it must be i = k and j = l, which is excluded by assumption. Therefore the nontrivial quadruples i,j,k,l satisfying both the equalities

E i -E j = E k -E l and E (2)
i -E

(2)

j = E (2)
k -E

(2) l are those for which

s(i) = s(j) = s(k) = s(l) and n(i) -n(j) = n(k) -n(l).
Let us now evaluate the terms E

(3) j

as in [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]. We have

E (3) j = m =j,p =j (E m -E j ) -1 (E p -E j ) -1 V j,m V m,p V p,j - m =j (E m -E j ) -2 V j,m V m,j V j,j .
Since V a,b = 0 only if s(a) = -s(b) it turns out that V j,m and V m,p are different from 0 only if s(j) = s(p) = -s(m), but then V p,j = 0. Thus, recalling that V j,j = 0, we have E

= 0 for every j ∈ N × {-1, 1}.

We are going to complete the proof that the set S i,j,k,l defined as in [START_REF] Rabi | Space quantization in a gyrating magnetic field[END_REF] has full measure by showing that if i, j, k, l ∈ N×{-1, 1} are such that (i,

j) = (k, l), i = j, n(i)-n(j) = n(k)-n(l) (which follows from E i -E j = E k -E l ) and s(i) = s(j) = s(k) = s(l) (which follows from E (2) i -E (2) j = E (2) k -E (2) l ), then E (4) i -E (4) j = E (4) k -E (4) l .
The general formula for E (4) j (see [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]) is

E (4) j = - m =j,p =j,q =j (E m -E j ) -1 (E p -E j ) -1 (E q -E j ) -1 V j,m V m,p V p,q V q,j (25)
+ m =j,p =j V j,j V j,m V m,p V p,j [(E m -E j ) -1 (E p -E j ) -2 + (E m -E j ) -2 (E p -E j ) -1 ] + m =j,p =j (E m -E j ) -2 (E p -E j ) -1 V j,m V m,j V j,p V p,j - m =j (E m -E j ) -3 V j,m V m,j V 2 j,j .
Since V j,j = 0, only the first and third term of the right-hand side must be evaluated.

Let us compute the first term in [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]. In order to avoid null terms we must assume s(j) = -s(m) = s(p) = -s(q) and thus n(j) = n(p). Therefore the only nonzero terms in the sum are given by n

(j) = n(m) + 1 = n(p) + 2 = n(q) + 1 (if n(j) > 1) and n(j) = n(m) -1 = n(p) -2 = n(q) -1. We have k =j,l =j,i =j (E k -E j ) -1 (E l -E j ) -1 (E i -E j ) -1 V j,k V k,l V l,i V i,j = (-ω -s(j)Ω) -1 (-2ω) -1 (-ω -s(j)Ω) -1 n(j) 2 n(j) -1 2 + (ω -s(j)Ω) -1 (2ω) -1 (ω -s(j)Ω) -1 n(j) + 1 2 n(j) + 2 2 .
Notice that the formula is correct also in the case where n(j) = 0 or n(j) = 1.

Let us now compute the third term in [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]. As before, to avoid null terms we assume s(j) = -s(m) = -s(p). The nonzero terms in the sum are given by n(m) = n(j) ± 1 and n(p) = n(j) ± 1 thus we have to sum four terms. We have

m =j,p =j (E m -E j ) -2 (E p -E j ) -1 V j,m V m,j V j,p V p,j = (-ω -s(j)Ω) -2 (-ω -s(j)Ω) -1 n(j) 2 2 + (-ω -s(j)Ω) -2 (ω -s(j)Ω) -1 n(j) 2
n(j) + 1 2 + (ω -s(j)Ω) -2 (-ω -s(j)Ω) -1 n(j) 2
n(j) + 1 2 + (ω -s(j)Ω) -2 (ω -s(j)Ω) -1 n(j) + 1 2 2 .
By summing up all the terms one sees that, for fixed s = s(j), the term E (4) j depends quadratically on n(j), i.e.

E (4) j = C 0 (s(j)) + C 1 (s(j))n(j) + C 2 (s(j))n(j) 2 ,
where the coefficient C 2 (s(j)) is given by

C 2 (s(j)) = s(j) Ω(ω 2 + 3Ω 2 ) 2(ω 2 -Ω 2 ) 3 = 0. So, if i, j, k, l are such that s(i) = s(j) = s(k) = s(l) = s and n(i) -n(j) = n(k) -n(l), we have E (4) i -E (4) j = E (4) k -E (4) l ⇐⇒ C 1 (s)(n(i) -n(j)) + C 2 (s)(n(i) 2 -n(j) 2 ) = C 1 (s)(n(k) -n(l)) + C 2 (s)(n(k) 2 -n(l) 2 ) ⇐⇒ C 2 (s)(n(i) 2 -n(j) 2 ) = C 2 (s)(n(k) 2 -n(l) 2 ) ⇐⇒ C 2 (s)(n(i) + n(j)) = C 2 (s)(n(k) + n(l)) ⇐⇒ n(i) = n(k) and n(j) = n(l).
This concludes the proof that for almost every g ∈ R and every i, j, k, l ∈ N × {-1, 1}, with (i, j) = (k, l) and i = j, one has The proof of Theorem 1 is then concluded, thanks to Theorem 3, if we show that the controlled Hamiltonian x ⊗ 1 couples, directly or indirectly, all the energy levels for almost all g ∈ R.

E g i -E g j = E g k -E g l . ( )/2 E j, 1 E j,1 (
More precisely, we show below that Φ g j , (x ⊗ 1)Φ g k = 0 for almost every g ∈ R for all j, k such that s(j) = s(k) and |n(j) -n(k)| = 1 or s(j) = -s(k) and n(j) = n(k). See Figure 2. As before, it is enough to show that the corresponding Taylor series in g is nonzero.

Set

Φ g j = Φ j + ∞ m=1 g m Φ (m)
j . We have

Φ j , (x ⊗ 1)Φ k = δ n(j),n(k)-1 n(j) + 1 2 + δ n(j),n(k)+1 n(j) 2 δ s(j),s(k) .
This is enough to say that Φ g j , (x ⊗ 1)Φ g k = 0 for almost every g for all j, k such that s(j) = s(k) and |n(j) -n(k)| = 1.

The term Φ

(1) j can be characterized through the relation

(H Rabi,0 + gx ⊗ σ 1 )(Φ j + gΦ (1) j + o(g)) = (E j + gE (1) j + o(g))(Φ j + gΦ (1) j + o(g)). (26)
Regrouping the first-order terms in [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF] we get

H Rabi,0 Φ (1) j + (x ⊗ σ 1 ) Φ j -E j Φ (1)
j -E

(1)

j Φ j = 0. ( 27
)
Denote by Π the orthogonal projection on the orthogonal complement to Φ j . Applying Π to [START_REF] Rouchon | Quantum systems and control[END_REF], we get (H Rabi,0 -E j 1)Φ

(1) j + Π(x ⊗ σ 1 )Φ j = 0. Notice that the orthogonal complement to Φ j is an invariant space for the operator H Rabi,0 -E j 1, which is invertible when restricted to it. We write (H Rabi,0 -E j 1) -1 to denote its inverse (whose values are in the orthogonal complement to Φ j ). Thus,

Φ (1) j = -(H Rabi,0 -E j 1) -1 Π(x ⊗ σ 1 )Φ j = l =j (E j -E l ) -1 Φ l , (x ⊗ σ 1 )Φ j Φ l .
The linear term in the Taylor expansion of Φ g j , (x ⊗ 1)Φ g k with respect to g is given by

Φ j , (x ⊗ 1)Φ (1) k + Φ (1) j , (x ⊗ 1)Φ k = l =k (E k -E l ) -1 Φ l , (x ⊗ σ 1 )Φ k Φ j , (x ⊗ 1)Φ l + l =j (E j -E l ) -1 Φ l , (x ⊗ σ 1 )Φ j Φ l , (x ⊗ 1)Φ k .
Taking into account only the nonzero terms in the first sum gives s(j) = s(l) = -s(k), |n(j) -n(l)| = 1 and |n(k) -n(l)| = 1, so for fixed j, k we only have (at most) two terms. We assume n(j) = n(k), thus we have

Φ j , (x ⊗ 1)Φ (1) k = (-ω + s(k)Ω) -1 n(k) + 1 2 + (ω + s(k)Ω) -1 n(k) 2 = ω + s(k)Ω(1 + 2n(k)) 2(Ω 2 -ω 2 ) .
Similarly, Φ

j , (x ⊗ 1)Φ k = ω -s(k)Ω(1 + 2n(k)) 2(Ω 2 -ω 2 ) . (1)
Hence, Φ

j , (x ⊗ 1)Φ k + Φ j , (x ⊗ 1)Φ

(1) k = ω Ω 2 -ω 2 ,
which is different from zero for every k. Thus the Taylor series of Φ g j , (x ⊗ 1)Φ g k is nonzero, which concludes the proof of the theorem in the case ω = Ω.

Case ω = Ω

Let us assume in this section that ω = Ω. Hence

E j,1 = E j+1,-1 = ω(j + 1)
for every j ∈ N and E 0,-1 = 0.

Recall that, as g varies in R, a complete set of eigenpairs for H Rabi can be given in the form (E g j , Φ g j ) j∈N×{-1,1} with each E g j and Φ g j depending analytically on g and E 0 j = E j . Notice that each Φ g j can be taken in L 2 (R, R) ⊗ C 2 and that we can take Φ 0 0,-1 = Φ 0,-1 = ϕ 0 ⊗ ν -1 . The eigenfunctions Φ 0 j,1 , Φ 0 j+1,-1 are characterized by the relation

Φ 0 j,1 , (X ⊗ σ 1 )Φ 0 j+1,-1 = 0.
(This follows by standard perturbation relations recalled in the Appendix, see equation (32).)

Since, moreover, Φ 0 j,1 and Φ 0 j+1,-1 are orthogonal and of norm one, one gets, up to signs,

{Φ 0 j,1 , Φ 0 j+1,-1 } = Φ j,1 + Φ j+1,-1 √ 2 , Φ j,1 -Φ j+1,-1 √ 2 .
Let us fix, by convention,

Φ 0 j,1 = Φ j,1 + Φ j+1,-1 √ 2 , Φ 0 j+1,-1 = Φ j,1 -Φ j+1,-1 √ 2 .
With this choice (32) yields

d dg E g j,1 g=0 = j + 1 2 , d dg E g j+1,-1 g=0 = - j + 1 2 .
Assume now that

E 0 j 1 -E 0 j 2 = E 0 j 3 -E 0 j 4 , (28)
d dg g=0 (E g j 1 -E g j 2 ) = d dg g=0 (E g j 3 -E g j 4 ), (29)
and let us prove that j 1 = j 3 and j 2 = j 4 . Write E 0 j l = ωk l for l = 1, . . . , 4 for some k 1 , . . . , k 4 ∈ N. Then [START_REF] Schuster | Resolving photon-number states in a superconducting circuit[END_REF] implies that k 1 = k 2 + p and k 3 = k 4 + p for some p ∈ Z.

Assume that k 1 , . . . , k 4 = 0 (the other cases being similar). Equality (29) then yields

s 1 k 2 + p -s 2 k 2 = s 3 k 4 + p -s 4 k 4 (30)
for some s 1 , . . . , s 4 ∈ {-1, 1} and we are left to prove that s 1 = s 3 , s 2 = s 4 and k 2 = k 4 .

Without loss of generality we can assume that k 2 + p = max{k 2 , k 2 + p, k 4 , k 4 + p} and that s 1 = 1. Then p and the left-hand side in [START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF] are nonnegative, which implies that s 3 = 1. We can then rewrite (30) as

k 2 + p -s 2 k 2 = k 4 + p -s 4 k 4 (31)
If s 2 = s 4 then, by monotonicity of k → √ k + p ± √ k we deduce that k 2 = k 4 and we are done. Finally, if s 2 = -s 4 then either s 2 = -1, and then k 4 > k 2 , or s 2 = 1 and then

√ k 2 + p = √ k 4 + p + √ k 4 + √ k 2 > √ p + √ k 2 ≥ √ k 2 + p, leading in both cases to a contradiction.
In order to apply Theorem 3 we are left to prove that the controlled Hamiltonian X ⊗ 1 couples, directly or indirectly, all the elements of the basis (Φ 0 j ) j∈N×{-1,1} . First notice that

Φ 0 j,1 , (X ⊗ 1)Φ 0 j+1,1 = 1 2 j + 1 2 + j + 2 2 = 0, j ∈ N,
meaning that X ⊗ 1 couples all basis elements of the type Φ 0 j,1 . Similarly, one has that Φ 0 j,-1 , (X ⊗ 1)Φ 0 j+1,-1 = 0 for every j ∈ N. Finally, one easily checks that Φ 0 j,1 , (X ⊗ 1)Φ 0 j,-1 = 0 for every j ∈ N, completing the proof.

Appendix: global real analyticity of eigenpairs

In this section we detail the proof of the analyticity of the eigenpairs (E g j , Φ g j ) as functions of g ∈ R. This result is somehow folklore and follows from well-known theorems in linear perturbation theory ( [START_REF] Kato | Perturbation theory for linear operators[END_REF][START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]). Since the latter are spread over the literature, we prefer to sketch here the main ideas of the argument, which is also useful in other contexts, as e.g. for the degenerate Gell-Man and Low theorem ( [START_REF] Brouder | Gell-Mann and Low formula for degenerate unperturbed states[END_REF][START_REF] Brouder | Many-body green function of degenerate systems[END_REF][START_REF] Brouder | Adiabatic approximation, Gell-Mann and Low theorem, and degeneracies: A pedagogical example[END_REF]).

Recall that, given two self-adjoint operators H 0 , W on a separable Hilbert space, W is said to be Kato-small with respect to H 0 if D(H 0 ) ⊂ D(W ) and for every a > 0 there exists b > 0 such that W ψ ≤ a H 0 ψ + b ψ for every ψ ∈ D(H 0 ). According to this definition, it is easy to check that V , defined as in [START_REF] Panati | Space-Adiabatic Perturbation Theory[END_REF], is Kato-small with respect to H Rabi,0 . Proposition 4. Let H 0 , W be self-adjoint operators on a separable Hilbert space H. Assume that W is Kato-small with respect to H 0 and that H 0 has compact resolvent. Then H g = H 0 + gW is self-adjoint with compact resolvent for every g ∈ R, hence it admits a complete orthonormal system of eigenfunctions {ψ j (g)} j∈N and the corresponding eigenvalues {E j (g)} j∈N have finite multiplicity and do not accumulate at any finite point. Moreover, up to a suitable choice of labelling, for each j ∈ N the function g → (E j (g), ψ j (g)) is analytic from R to R × H. Proof. Self-adjointness of H g follows from [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]Theorem X.12]. Let us first prove that H g has compact resolvent for every g ∈ R (we provide a direct proof; see also [START_REF] Kato | Perturbation theory for linear operators[END_REF]Theorem VII.2.4]). Since (2) (H g -i1) -1 = (H 0 -i1) -1 1 -g W (H g -i1) -1

and (H 0 -i1) -1 is compact, we are left to prove that W (H g -i1) -1 is bounded. Notice that for every a > 0 there exists b > 0 such that W ψ ≤ a (H g -i1)ψ + b ψ for every ψ ∈ D(H 0 ) = D(H g -i1). Setting ψ = (H g -i1) -1 φ for φ ∈ H, we get

W (H g -i1) -1 φ ≤ a φ + b (H g -i1) -1 φ ≤ (a + b) φ ,
where the last inequality follows from (H g -i1)η η ≥ | (H g -i1)η, η | ≥ η 2 , for η ∈ D(H g -i1).

If E j (g) is simple then, by Rellich's theorem [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]Theorem XII.8], up to a reordering of the spectrum of H g for g in a neighborhood of g, the eigenpair parameterization g → (E j (g ), ψ j (g )) is analytic near g. The previous result can be generalized to the case of eigenvalues of multiplicity m ≥ 2. Indeed, if E j 1 (g) = • • • = E jm (g) and j 1 , . . . , j m are distinct, we can assume, up to relabelling, that g → (E j 1 (g ), ψ j 1 (g )), . . . , g → (E jm (g ), ψ jm (g )) are analytic near g ([26, Theorem XII.13]). Moreover, ψ j l (g), W ψ j k (g) = δ lk d dg E j l (g), l, k = 1, . . . , m.

(See, for instance, [START_REF] Albert | Genericity of simple eigenvalues for elliptic PDE's[END_REF].) In order to describe the spectrum by a countable family of functions which are globally analytic on R, we should ensure that all such locally analytic functions E j (•) can be extended indefinitely.

Fix j and let I be the maximal interval containing 0 such that, up to relabelling, E j is analytic on I. Let us prove that, for every ∈ (0, +∞), [-, ] ⊂ I.

Set a = (2 ) -1 . Take b > 0 such that W ψ ≤ a H 0 ψ + b ψ for every ψ ∈ D(H 0 ). Hence, for ψ ∈ D(H 0 ) and g ∈ [-, ], It follows from (32) that, for every g ∈ I ∩ [-, ],

d dg E j (g) = | ψ j (g), W ψ j (g) | ≤ W ψ j (g) ≤ 2a H g ψ j (g) + 2b ψ j (g) . A standard application of Gronwall inequality yields that E j (•) is Lipschitz on I ∩ [-, ]. In particular the limits of E j (•) at the boundary of I ∩ [-, ] are well defined and they are eigenvalues of the corresponding perturbed operators (see e.g. [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]Theorem XII.7]). Therefore E j (•) can be extended indefinitely on [-, ] and thus also on the whole real line.

Figure 1 :

 1 Figure 1: Each vertex of the graph represents an eigenstate of A (when the spectrum is not simple, several nodes may be attached to the same eigenvalue). An edge links two vertices if and only if B connects the corresponding eigenstates. In this example, φ 1 , Bφ 2 and φ 1 , Bφ 3 are not zero, while φ 1 , Bφ 4 = φ 2 , Bφ 3 = φ 2 , Bφ 4 = 0.

4. 1 . 1

 11 Step I: Relevant eigenvalue pairs are non-resonant.

+ )/ 2 0Figure 2 :

 22 Figure 2: The dashed lines connect eigenvalues of H Rabi that are "coupled" by the controlled Hamiltonian when g = 0, while the dotted lines connect eigenvalues that are coupled by the controlled Hamiltonian for almost all g = 0.

  W ψ ≤ a (H 0 + gW )ψ + a|g| W ψ + b ψ ≤ a H g ψ + W ψ 2 + b ψand thereforeW ψ ≤ 2a H g ψ + 2b ψ , g ∈ [-, ].

  Hence,d dg |E j (g)| ≤ 2a|E j (g)| + 2b, for almost every g ∈ I ∩ [-, ].

[START_REF] Albert | Genericity of simple eigenvalues for elliptic PDE's[END_REF] Notice that (H g -i1) -1 maps H into D(H g ) ⊂ D(W ), so that W (H g -i1) -1 is well-defined on the whole Hilbert space H.

Acknowledgements

We are grateful to L. Pinna for a careful reading of the manuscript and some useful comments, and to the Institut Henri Poincaré for the kind hospitality in the framework of the trimester "Variational and Spectral Methods in Quantum Mechanics ", organized by M. J. Esteban and M. Lewin. This research has been supported by the European Research Council, ERC StG 2009 "GeCoMethods", contract number 239748, by the iCODE institute, research project of the Idex Paris-Saclay and by the EU FP7 project QUAINT, grant agreement no. 297861.