Workshop report: Photopatch testing—methods and indications

J. Fluid Mech. (2011), vol. 689, pp. 376–416. doi:10.1017/jfm.2011.421

c Cambridge University Press 2011

376

The onset of strongly localized thermal convection in rotating spherical shells Andrew P. Bassom1 , Andrew M. Soward2,3 † and Sergey V. Starchenko4 1 School of Mathematics and Statistics, The University of Western Australia,

Crawley 6009, Australia 2 School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne

NE1 7RU, UK 3 Mathematics Research Institute, University of Exeter, Exeter EX4 4QF, UK 4 IZMIRAN, Troitsk, Moscow Region, 142190, Russia

(Received 15 October 2010; revised 14 June 2011; accepted 22 September 2011; first published online 16 November 2011)

A Boussinesq fluid of kinematic velocity ν and thermal diffusivity κ is confined within a rapidly rotating shell with inner and outer sphere boundary radii ri∗ and ro∗ , respectively. The boundaries of the shell corotate at angular velocity Ωio∗ and a continuously varying stratification profile is applied which is unstable in 1/2 ri∗ < r∗ < rn∗ ≡ (1 + ε2 ) ri∗ and stable in rn∗ < r∗ < ro∗ . When ε  1, the unstable zone attached to the inner boundary is thin. As in previous small Ekman number E = ν/(2Ωio∗ ri∗ 2 ) studies, convection at the onset of instability takes on the familiar ‘cartridge belt’ structure, which is localized within a narrow layer adjacent to, but outside, the cylinder tangent to the inner sphere at its equator (Dormy et al. J. Fluid 2/9 Mech., 2004, vol. 501, pp. 43–70), with estimated radial width of order (Eε4 ) ri∗ . The azimuthally propagating convective columns, described by the cartridge belt, reside entirely within the unstable layer when E1/5  ε, and extend from the equatorial plane an axial distance εri∗ along the tangent cylinder as far as its intersection with the neutrally stable spherical surface r∗ = rn∗ . We investigate the eigensolutions of the ordinary differential equation governing the axial structure of the cartridge belt both numerically for moderate-to-small values of the stratification parameter ε and analytically when ε  1. At the lowest order of the expansion in powers of ε, the eigenmodes resemble those for classical plane layer convection, being either steady (exchange of stabilities) or, for small Prandtl number P ≡ ν/κ < 1, oscillatory (overstability) with a frequency ±Ω ∗ . At the next order, the axial variation of the basic state removes any plane layer degeneracies. First, the exchange of stabilities modes oscillate at a low frequency causing the short axial columns to propagate as a wave with a small angular velocity, termed slow modes. Second, the magnitudes of both the Rayleigh number and frequency of the two overstable modes, termed fast modes, split. When P < 1 the slow modes that exist at large azimuthal wavenumbers M make a continuous transition to the preferred fast modes at small M. At all values of P the critical Rayleigh number corresponds to a mode exhibiting prograde propagation, whether it be a fast or slow mode. This feature is shared by the uniform classical convective shell models, as well as Busse’s celebrated annulus model. None of them possess any stable stratification and typically are prone to easily excitable Rossby or inertial modes of convection at small P. By way of contrast these structures cannot † Email address for correspondence: [email protected]

Localized rotating convection

377

exist in our model for small ε due to the viscous damping in the outer thick stable region. Key words: buoyancy-driven instability, B´enard convection

1. Introduction The study of Boussinesq rotating convection in self-gravitating spheres and shells at a small Ekman number E has developed a long way since the pioneering studies of Roberts (1965, 1968) and Busse (1970a). That early work considered a uniform distribution of heat sources throughout the fluid and the key conclusion was that, provided that the inner solid core is sufficiently small, the onset of convection is localized at midlatitudes. Furthermore, this activity takes place in cells arranged on a cylinder whose generator is parallel with the rotation axis. This ‘cartridge belt’ therefore consists of axial rolls whose structure varies slowly between the equatorial plane and the intersection of the cylinder with the outer boundary. Busse (1970a) showed that the preferred mode of convection at onset has an axial velocity that is antisymmetric about the equatorial plane while the velocity normal to the rotation axis is symmetric with respect to the same plane. The rolls have a short E1/3 azimuthal length scale and propagate prograde with a Rossby wave character. Soward (1977) pointed out that in the E  1 limit the detail of the radial (that is normal to the rotation axis) structure is not straightforward and its correct form was only finally resolved much later by Jones, Soward & Mussa (2000). Their asymptotic solution describes rolls that extend radially in the direction normal to the cylinder on a relatively long E1/6 length scale; moreover these elongated structures also tilt progradely with increasing distance from the axis. The case of differential heating, when there are no heat sources in the fluid but instead a temperature difference is maintained between the inner and outer boundaries, was considered by Dormy et al. (2004). Then relatively large (at least compared with the internal heating structure) temperature gradients occur close to the inner boundary. As a result the cartridge belt now resides on the cylinder tangent to the inner boundary at its equator; this is commonly referred to as the tangent cylinder. The rolls now have the radial length scale E2/9 and tilt as before (see, for example, the equatorial cross-section portrayed in figure 4 of Dormy et al. 2004 for E = 10−7 ). Both Jones et al. (2000) and Dormy et al. (2004) include extensive discussions of high-resolution numerical results for convection at very small E. Unfortunately, similar computations for the onset of convection in thin shells have yet to be performed. Nevertheless, based on the description above, convection adjacent to the tangent cylinder is to be expected both in the case of internal and of differential heating; such motions have been found (see, for example, figure 6 of Takehiro & Hayashi 1995 for E = 10−3 , or figure 4 of Al-Shamali, Heimpel & Aurnou 2004 for E = 3 × 10−4 ). Certainly this general picture is now widely accepted (see figure 7 of Stanley et al. 2005, or figure 1a in Stanley & Glatzmaier 2010). Indeed the slow rotation case at large E was considered much earlier by Busse (1970b). Some cautionary remarks are appropriate here for it is well known, even in the classic configuration of uniform internal heating, that modes of convection with longer length scales attached to the outer boundary are possible at small Prandtl number P. This topic has received considerable attention recently (Zhang, Liao & Busse 2007;

378

A. P. Bassom, A. M. Soward and S. V. Starchenko

Net, Garcia & S´anchez 2008). Nevertheless, Zhang & Liao (2004) noted that the traditional cartridge-belt theory, relying on the short length scale approximations, applies for P  1 provided that E  P, a limit which is generally achieved by physically realistic models. 1.1. Applications to stellar convection zones and planetary convective cores The main reason for the continuing interest in rotating convection in both thick and thin shells is the possible application to stellar convection zones and planetary convective cores, for which nonlinear aspects are generally important. This includes the generation of differential rotation in the Sun as investigated by Busse (1970b) or in the giant planets (Aurnou & Olson 2001; Stanley & Glatzmaier 2010; Takehiro, Yamada & Hayashi 2011), while there are also many studies of planetary dynamo applications (Wicht & Tilgner 2010). The simple Boussinesq models of thermal convection are highly idealized. For example, in planetary applications it is thought that both thermal and compositional convection occur (Fearn & Loper 1981). Investigations of the conditions expected in the Earth’s core have been pursued by Braginsky & Roberts (1995) and Busse (2002) has found interesting results with two buoyancy ingredients on applying his annulus model (Busse 1970a). Indeed the Boussinesq approximation itself is inadequate to capture some features and the anaelastic approximation is often better as, for example, in applications to the giant planets (Stanley & Glatzmaier 2010). We remark in passing that Jones, Kuzanyan & Mitchell (2009) have undertaken the anaelastic extension of the Boussinesq study by Jones et al. (2000). Nevertheless, much of the basic dynamics can be investigated within the simple thermal Boussinesq model provided a non-uniform temperature gradient is imposed, which generally describes regions of both stable and unstable stratification (see, for example, Starchenko, Kotelnikova & ˇ Maslov 2006 or Simkanin, Hejda & Jankoviˇcov´a 2009). Such magnetohydrodynamic ˇ ˇ c´ık convection in a plane layer has been considered by Simkanin, Brestensk´y & Sevˇ (2003, 2006). Applications of partially stable planetary cores to the dynamos of Mercury and Saturn have been investigated by Christensen (2006) and by Christensen & Wicht (2008), while further planetary dynamo applications abound (see, for example, Stanley & Mohammadi 2008; Stanley & Glatzmaier 2010; Wicht & Tilgner 2010 and references therein). Planetary dynamo models with a thin unstable region lying above a thick stable region have been considered by Stanley & Bloxham (2004, 2006) as models for the ice giants Uranus and Neptune, as well as for the terrestrial planet Mercury (Stanley et al. 2005). Although in this paper we are also interested in thin unstable regions, our focus is with the opposite case where the unstable region lies below the stable region and this leads to the interesting phenomenon of penetrative convection. Takehiro & Lister (2001, 2002) noted its relevance to the Earth’s fluid core, where a thin stable region resides on top of the convective interior (Braginsky 1993); it is possible that penetrative convection also occurs in the Jovian atmosphere where the unstable and stable layers may be of similar thicknesses. Previously Zhang & Schubert (2000, 2002) identified a curious form of convection that can be set up when the Prandtl number P is small. They found that, although convection occurred in the inner unstable region, additional quite vigorous convection took place near the outer boundary well within the stable region. Zhang & Schubert (2000) coined the term ‘teleconvection’ for this effect. In this paper, we consider a model with continuously varying stratification dependent on the radius r∗ ; the inner and outer sphere boundaries are taken to

Localized rotating convection ∗

ri∗

379

ro∗ , ∗

be r = and respectively. The stratification is assumed to be unstable in a thin layer ri∗ < r < rn∗ adjacent to the inner boundary; above it, in the remaining thick layer rn∗ < r∗ < ro∗ , the fluid is strongly stably stratified. Relative to cylindrical polar coordinates s∗ , φ, z∗ , the tangent√cylinder s∗ = s∗t = ri∗ intersects the neutrally stable surface at z∗ = z∗n ≡ εri∗ (i.e. rn∗ = 1 + ε 2 ri∗ ), where ε is generally small. When E  1, and also E  P, convection at the onset of instability is manifest as travelling waves with real frequency Ω ∗ and possessing a large azimuthal wavenumber M, which are localized on the tangent cylinder exactly as described by Dormy et al. (2004). The only, but crucial, difference is that the convecting columns are now limited in axial extent to small |z∗ | = O(εri∗ ) due to the stable stratification of the exterior fluid |z∗ | > z∗n . The convection resembles that within a thin shell, save that the outer boundary r∗ = rn∗ can be thought of as being ‘soft’. We adopt this term to contrast with the situation for the usual thin shell model when both the inner and outer boundaries are ‘hard’, in the sense that the normal velocity to them vanishes. In our model the velocity simply decays to zero on leaving the unstably stratified region. 1.2. Outline of the paper Since convection is largely independent of conditions at the outer boundary r∗ = ro∗ , we adopt the inner boundary radius ri∗ as our unit of length in our non-dimensionalization of the problem in § 2. As this is in contrast to the more usual choice ro∗ made by Dormy et al. (2004), we compare the two formulations for the common case of differential heating in Appendix C (which is contained in the supplementary material available at journals.cambridge.org/flm for this paper). In § 2.3 we consider numerical solutions of the ordinary differential equations that describe the axial structure of the eigenmodes at onset. As anticipated above their axial extent |z∗ | = O(εri∗ ) decreases in concert with ε; the geometry of the convective region is identified in figure 1 labelled with the dimensionless variables introduced in § 2 below. Guided by those numerical results we construct asymptotic solutions of the ordinary differential equations in § 3, which build on the assumption ε  1. A complete understanding of the nature of the solutions requires taking our expansions in terms of the stratification parameter ε to two (zeroth and first) orders. The convective modes can then be of two types distinguished by the size of their frequency. One is a low-frequency (or slow phase angular velocity Ω ∗ /M) mode investigated in § 4; the other is a high-frequency (or fast) mode explored in § 5, which only exists for sufficiently small Prandtl number. The results of our zeroth-order theory detailed in §§ 3.1, 4.1 and 5.1 resemble those of the classical plane layer geometry described by Chandrasekhar (1961, Ch. 3), for which there exist steady (or exchange of stabilities) modes (linked to our slow modes) and oscillatory (or overstable) modes (related to our fast modes) with both positive and negative frequencies, say ±Ω ∗ , but with the same Rayleigh number. (In Appendix D within the supplementary material we confirm that the preferred mode with smallest Rayleigh number corresponds to our assumed symmetry; Busse 1970a.) As the azimuthal wavenumber M is decreased from large values the oscillatory modes ±Ω ∗ 6= 0 come into existence at some M = Mt . This is manifest by a pitchfork of the plot of Ω ∗ versus M, which we will refer to as a ‘frequency pitchfork’. For given M (< Mt ), the oscillatory mode of convection always possesses a Rayleigh number smaller than that for the steady mode. The overall conclusion is that the critical Rayleigh number corresponds to the oscillatory mode for sufficiently small P, less than Pc (≈ 0.65), and to the steady mode for P > Pc . The sharp distinction described between slow and fast modes is a feature of our small-ε theory.

380

A. P. Bassom, A. M. Soward and S. V. Starchenko 3

2 st

st

z N

1

0 –1

E 0

1

2

3

s

F IGURE 1. The shell geometry drawn for η = 0.35 used for the numerical results reported in § 2.3 for the particular case ε = 1. The inner and outer boundaries are indicated by the continuous circular arcs. The tangent cylinder s = st (long-dashed) touches the equator of the inner sphere at E (st , 0) and intersects the unstable–stable stratification interface r = rn at N (st , ε). The radial extent of the convective column is estimated by the cylinder s = st + δt , 2/9 where δt = (Eε4 ) , drawn dot-dashed for E = 10−4 .

The first-order analyses of §§ 3.2, 4.2 and 5.2 resolve the degeneracies of the zeroth-order results and the calculations of § 4.2 fix the frequency of the slow modes. A first hint that the frequency pitchfork at M = Mt from slow to fast modes is actually imperfect, arises in § 4.3. The detailed nature of the transition, from the slow modes for M ' Mt to the preferred fast modes for M / Mt , becomes apparent in § 5.1.3 (see also Appendix E in the supplementary material). The results of § 5.2 distinguish which of the zeroth-order frequencies ±Ω ∗ for the fast modes actually corresponds to the lower Rayleigh number. The smooth transition on varying M from large (slow modes) to small (fast modes) values is evidenced by the numerical results for small but finite ε described in § 2.3 and illustrated in figures 4. This smooth transition is akin to the results of Busse’s (1970a) annulus model, as we discuss in Appendix B.3. In contrast, on decreasing P our critical mode suffers a sudden jump in structure from slow to fast at P ≈ Pc , where the slow and fast mode minimum Rayleigh numbers coincide (see the asymptotic results portrayed in figure 4c). This discontinuous feature does not seem to arise in the more classical models of uniform internal heating and differential heating in spheres, shells and annuli (Busse 1970a; Jones et al. 2000; Dormy et al. 2004). With that proviso, our overall conclusion, that for all values of P the critical mode possesses a positive angular phase velocity describing prograde propagation (see also Appendix A), is consistent with the classical models. At small P, however, our comparison with the classical models in Appendix B uncovers a very remarkable difference which arises due to the distinction between ‘hard’ and ‘soft’ boundaries described at the end of § 1.1. This behaviour emerges from

Localized rotating convection

381

the asymptotic results of § 6, particularly the zeroth-order results of § 6.1. There it becomes apparent that the buoyancy force is present in the leading order force balance which does not occur in any of the pure uniform internal heating and differential heating models. In those, Rossby or inertial waves provide the leading force balance and the buoyancy force is only invoked at the next order so as to maintain the wave motion subject to small viscous dissipation. In our case, waves generated in the thin unstable layer propagate out into the stable layer above, where they are damped. That is why the buoyancy force in our model must be invoked at lowest order with the consequence that the Rayleigh number needed to sustain convection subject to our ‘soft’ outer boundary is an order of magnitude O(P−1 ) larger than a ‘hard’ outer boundary model would predict. Other aspects of our findings, such as the penetration depth of the convection into the overlying stable region, are finally summarized and discussed in § 7. However, we emphasize here (rather than later) that our curious small P results may shed new light on the teleconvection issue identified by Zhang & Schubert (2000, 2002). Since their unstable region is not asymptotically thin like ours, our results cannot be applied directly. Nevertheless, our results seem to suggest that convection must be driven somewhat harder than one might expect because of the effective upper ‘soft’ boundary. It may simply be that inertial waves driven quite hard in the unstable region focus near the outer boundary leading to the teleconvection. 2. Formulation We consider a shell of inner and outer radii ri∗ and ro∗ , respectively, which rotates at the constant angular velocity Ωio∗ of the bounding solids and is subject to the radial gravity field −gr∗ , where r∗ is the position vector. It is filled with Boussinesq fluid with coefficient of thermal expansion α, thermal conductivity κ and kinematic viscosity ν. The spherically symmetric basic density distribution is unstably stratified inside the inner shell ri∗ < r∗ < rn∗ and stably stratified in the remaining rn∗ < r∗ < ro∗ . We model that density by a temperature distribution T(r∗ ) whose gradient is ( < 0 (ri∗ < r∗ < rn∗ ), dT ∗ ∗ ∗ 0 (2.1a) T ≡ ∗ = −r Q (r ) > 0 (rn∗ < r∗ < ro∗ ). dr

The classical uniform heat source model concerns positive Q ∗ = constant and that of differential heating positive Q ∗ ∝ r∗ −3 . Since these two cases have been well studied in the small Ekman number limit (Jones et al. 2000; Dormy et al. 2004, respectively), we consider instead the choice

Q ∗ ∝ (r∗ −3 − rn∗ −3 ),

(2.1b)

consistent with (2.1a), for a basic state of mixed unstable–stable smoothly varying stratification. There may be different origins for the buoyancy force Rθ r in our dimensionless equation of motion (2.3a) below, such as thermal and compositional gradients, or it may be due to the compressibility of the fluid and associated thermodynamic processes that may need entropy interpretations. Whatever the source of buoyancy, its redistribution within the fluid is generally described by an advectivediffusion process as adopted by us in (2.3d) below. Furthermore, in the limit δn∗ ≡ rn∗ − ri∗ ↓ 0 addressed by our later asymptotic analysis, we only consider the first two terms of the Taylor expansion of Q ∗ , i.e. Q ∗ ≈ Qi∗ (rn∗ − r∗ )/δn∗ (see (2.15a,b) below), which is linear in r∗ and so essentially insensitive to the detailed nature of the applied stratification. (The subscripts i , n , o are used to denote values at r = ri , rn , ro

382

A. P. Bassom, A. M. Soward and S. V. Starchenko

respectively.) That said, Starchenko et al. (2006, equation (32)) proposed a buoyancy source that corresponds to that defined by (2.1a,b) on the basis of an almost adiabatic approximation. This proposal suggests the specific choice (2.1b), whereas our asymptotic analysis is based on the generic limit (2.15a). The linearized equations governing perturbations to the basic state are rendered dimensionless by adopting the inner sphere radius ri∗ for the unit of length and the viscous scale ri∗ 2 /ν for the unit of time, while for our unit of temperature we take −νri∗ Ti0 /κ, where Ti0 (< 0) is the temperature gradient at the inner boundary. Upon introducing the dimensionless parameters E=

ν , 2Ωio∗ ri∗ 2

R=−

gαTi0 ri∗ 5 , νκ

P=

ν , κ

η=

ri∗ , ro∗

(2.2a,b,c,d)

namely the Ekman, Rayleigh, Prandtl numbers and aspect ratio, respectively, our linearized equations of motion and continuity for the velocity v and temperature θ become ∂v + E−1 zˆ × v = −∇p + Rθr + ∇ 2 v, ∂t

∇ · v = 0,

(2.3a,b)

where zˆ is the unit vector parallel to the rotation axis. From these the vorticity equation ∂ ∂v ∇ × v − E−1 = −Rr × ∇θ + ∇ 2 ∇ × v ∂t ∂z

(2.3c)

follows, while the heat conduction equation for the temperature θ reduces to P

∂θ = Q v · r + ∇ 2 θ, ∂t

Q (r) ≡

Q ∗ r−3 − rn−3 = . Qi∗ 1 − rn−3

(2.3d,e)

The system (2.3a,b,d,e) corresponds to Starchenko et al. (2006, equations (26), (28), (32)), in which −∇θs is proportional to our Q (r)r). We solve the governing linear equations (2.3a–e) subject to isothermal θ = 0 and impermeable r · v = 0 boundary conditions on the inner and outer surfaces r = ri = 1 and r = ro = η−1 . The geometry of our flow domain is illustrated in figure 1. Since our interest is in convection close to the inner core boundary, the precise location of the outer boundary has little influence on the solution. It is for this reason that our non-dimensionalization is based on conditions at the inner boundary (see (2.2a,b)) and (2.3e)) rather than quantities at the outer boundary. In view of our slightly non-standard choice we relate our variables to those employed by Dormy et al. (2004) in Appendix C within the supplementary material. Furthermore, as the most unstable modes are located on the inner sphere tangent cylinder the height q ε = rn2 − 1, (2.4) of the column EN in figure 1, where Q > 0, is of considerable importance. 2.1. The small Ekman number asymptotic analysis In the small Ekman number limit

E  1,

(2.5)

Localized rotating convection

383

it is convenient to adopt cylindrical polar coordinates s, φ, z and use the axial toroidal–poloidal decomposition v = ∇ × Ψ zˆ + ∇ × ∇ × Φˆz

(2.6)

for the velocity. Separable solutions are sought proportional to

E˜ (φ, t) ≡ exp[i(Mφ − Ωt)].

(2.7)

The standard method for solving the remaining two-dimensional problem in s and z is to separate the variables again and seek solutions of the type S(s)Z(z). Such a procedure cannot satisfy the governing equations exactly for two reasons: first, the driving temperature gradient, as measured by Q , is dependent on both s and z and, second, the top and bottom boundaries z = ±zo (s), where q zo (s) = ro2 − s2 , ro = η−1 , (2.8a,b) are clearly s-dependent and R non-parallel. This difficulty is overcome by making a WKB-ansatz S(s) = exp(i K(s) ds), so that, for example, the Laplacian becomes ∇2 =

∂2 − A2 , ∂z2

where A2 = (M/s)2 +K 2 .

(2.9a,b)

This leads to an ordinary differential equation for Z(z) whose solution determines a local eigenvalue problem for the radial wavenumber K. This procedure can only be justified when the radial s-length scale is short (K large); in other words it is necessary that the s variations of Q and zo are negligible on that short radial scale. On application of the WKB-ansatz two types of solution are possible. One form, typified by the internal heating case of Jones et al. (2000), is a mode localized in the neighbourhood of some cylinder s = sM outside the tangent cylinder s = st = 1 of the inner sphere. The other structure, characterized by the differential heating problem studied by Dormy et al. (2004), is a mode localized adjacent (but exterior) to the tangent cylinder. The complete asymptotic solutions of these problems are rather technical and the only issue that we need to take on board is that where the localization occurs the actual radial scale lengthens and is manifest by the approximation K ≈ 0. The detailed radial structure requires a higher-order theory which for the case sM 6= st needs to be addressed to resolve the eigenvalue problem (see Jones et al. 2000). Since our problem is dominated by differential heating, we reasonably anticipate that the tangent cylinder modes considered by Dormy et al. (2004), for which the local radial structure is determined by Airy functions with complex argument, pertain to our case. Then, since the higher-order theory is not needed to resolve the eigenvalue problem, we do not pursue the radial structure in this paper, other than to make the estimates (2.18a,b) that lead to restriction E1/5  ε (see (2.19c) below) on the validity of our strong stratification theory formulated in § 2.2. The leading-order approximation in powers of E of the vorticity equation (2.3c) is E−1 ∂v/∂z ≈ 0. To resolve the ensuing geostrophic degeneracy we seek short length scale quasigeostrophic solutions localized radially in the vicinity of the tangent cylinder s = st = 1, for which Ψ and Φ vary rapidly in the azimuthal φ-direction, somewhat less rapidly in the radial s-direction but slowly in the axial z-direction such

384

A. P. Bassom, A. M. Soward and S. V. Starchenko

that the respective length scales δφ , δt , δz satisfy ( ) zo (1) δφ  δt  δz = min implying ε

1 ∂ ∂ ∂   , s ∂φ ∂s ∂z

(2.10a,b)

where we have noted that the column EN in figure 1, drawn for ε < zo (1), has height zn (1) = ε (see (2.4)). These assumptions allow us to invoke the quasigeostrophic approximation v ≈ ∇ × Ψ zˆ + W zˆ

supplemented by

with W = zˆ · ∇ × ∇ × Φˆz ≈ A2 Φ,

∇ × v ≈ Υ zˆ + ∇ × W zˆ with

Υ = zˆ · ∇ × ∇ × Ψ zˆ ≈ A2 Ψ,

(2.11a,b) (2.11c,d)

in which z-derivatives are neglected. This provides us with sufficient accuracy for our asymptotic analysis and is consistent with the approximations made by Roberts (1968), Busse (1970a) and subsequent researchers. When we make the additional assumption K  M/s (see (2.10b)), our approximations (2.11a–d) lead to the expressions r · v = sU + zW

with sU ≈ iMΨ

(2.12a,b)

for the radial velocity, where U is the velocity in the radial s-direction, and r · ∇ × v ≈ zA2 Ψ + iMW

with

A2 ≈ M 2 /s2

(2.12c,d)

(see (2.9b)) for the radial vorticity. A conventional Ekman layer of width O(E1/2 ) (Greenspan 1968) forms at the intersection of the tangent cylinder with the outer boundary. The Ekman suction velocity produced by this layer depends on whether the boundary is rigid or stress free. In either case its magnitude is negligible at leading order in powers of E and can safely be ignored. Accordingly, since r · v ≈ iMΨ + zW (see (2.12a,b)), it follows from (2.14c) below that the remaining isothermal and impermeable boundary conditions reduce to the single requirement r ·v ≈ 0

(2.13)

at r = ro on the mainstream flow exterior to the Ekman layers. In our E  1 quasigeostrophic limit, axial z-components of the vorticity equation (2.3c) and its curl, together with the heat conduction equation (2.3d) yield, correct to lowest order in E, the well-known equations dW = −iMERθ, dz dΨ E(A2 − iΩ)W − = zERθ, dz (A2 − iPΩ)θ = Q (iMΨ + zW),

E(A2 − iΩ)(A2 Ψ ) −

(2.14a) (2.14b) (2.14c)

respectively, where the only z-derivatives retained stem directly from the large Coriolis acceleration. We note that (2.14b) is simply the axial component of the equation of motion under the quasigeostrophic approximation E−1 zˆ × v ≈ −∇p: Ep ≈ −Ψ.

(2.14d)

We solve the coupled equations (2.14a–c) in the Northern Hemisphere 0 < z < zo (s) subject to the boundary condition r · v = 0 on z = zo (s) (see (2.13)) and the

385

Localized rotating convection

symmetry condition W = 0 on the equatorial plane z = 0. (In Appendix D within the supplementary material we verify that this is the appropriate symmetry choice since modes with other characteristics are more stable.) As has already been explained, we seek solutions localized outside but adjacent to the tangent cylinder so that sM = st = 1. On the tangent cylinder the fluid is only unstably stratified (Q > 0) within the column 0 < z < δz (see (2.10a)); elsewhere, when ε < zo (1), the stratification is stable (Q < 0) for zn (1) = ε 6 z < zo (1). 2.2. The scaling for strong stratification; ε  1 For the strongly stratified case, ε  1, the axial extent 0 < z < ε of the unstable layer is small. Indeed, when s = 1 + O(ε2 ) and z = O(ε), (2.3e) reduces to

Q ≈ (rn − r)/δn ≈ 1 − [2(s − 1) + z2 ]/ε2 ,

δn = rn − 1,

(2.15a,b)

which suggests the introduction of the stretched axial coordinate z ≡ εζ.

(2.16)

Like both Jones et al. (2000) and Dormy et al. (2004), we scale our other variables with respect to E but, unlike them, introduce a further scaling with respect to ε: R ≡ (Eε)−4/3 R , θ ≡ θ,

(M, A) ≡ (Eε)−1/3 (m, a),

(U, W) ≡ (Eε)−2/3 (u, w),

Ω ≡ (Eε)−2/3 ω, (2.17a,b,c)

Ψ ≡ (Eε)−1/3 ψ

(2.17d,e,f )

(see (C 7a–f ) within the supplementary material for the relation to the variables of Dormy et al. 2004). The scaling for M anticipates a φ length scale δφ = (Eε)1/3 , which is short compared with the z length scale δz = ε as required by (2.10a), when E1/2  ε. However, this by itself is not sufficient to ensure the validity of our analysis and a tighter restriction needs to be fulfilled, which emerges from the higher-order amplitude equation for the radial s structure. On regarding ω at fixed R as a complex function of s and m, arguments similar to those of Dormy et al. (2004, p. 52) lead to the estimates ω,s (s − 1) ∼ (Eε)2/3 ω,mm

∂2 , ∂s2

ω,s Q,s 1 ∼ ∼ = O(ε −2 ) ω,mm Q δn

(2.18a,b)

(see (2.15)) at s = st = 1 (ω,s ≡ ∂ω/∂s etc.). The former balance (2.18a) is the direct origin of the Airy function structure mentioned earlier. Together (2.18a,b) determine 2/9 the s length scale s − 1 ∼ δt = (Eε 4 ) of the convective columnar rolls, equivalently 1/2

(δt /ε2 )

= δφ /δt = (E/ε 5 )

1/9

.

(2.19a)

Our asymptotic analysis is only valid if δt is short compared with the radial width δn = O(ε 2 ) of the unstable region at z = 0 (see (2.15)), a condition that must be supplemented to (2.10a) giving δφ  δt  δn ( δz ).

(2.19b)

E1/5  ε  1.

(2.19c)

From (2.19a) it follows that both requirements δφ  δt and δt  δn are met when As we are only interested in moderately small ε, the condition (2.19c) is readily met in physically interesting situations, for which E is very small indeed. So, throughout our analysis, we envisage a limit in which E ↓ 0 independent of the size of ε and, for that matter, the Prandtl number P.

386

A. P. Bassom, A. M. Soward and S. V. Starchenko

From (2.14c) the temperature perturbation may be expressed as θ = Q (su + zw)/(a2 − iPω)

in terms of the scaled radial velocity

with a2 ≈ m2 /s2 ,

(2.20a,b)

(Eε)2/3 r · v = su + εζ w with su ≈ imψ.

(2.20c,d)

The substitution of (2.20a,b) into the equations of motion (2.14a,b) determines dψ = −iεmC ζ ψ − B w, dζ

dw = iεmC ζ w − F ψ, dζ

(2.21a,b)

where

C=

RQ , a2 − iPω

B = ε2 C ζ 2 − (a2 − iω),

F = m2 C − a2 (a2 − iω).

(2.21c,d) (2.21e)

For modes with the symmetry of interest w=

dψ =0 dζ

at ζ = 0,

(2.22a)

while according to (2.20c,d) the condition su + εζ w = 0 on the outer boundary r = ro is met when q imψ + εζ w = 0 at ζ = ε−1 zo (s) = ε−1 ro2 − s2 . (2.22b) 2.3. Numerical results for the tangent cylinder neutral modes Since the neutrally stable modes (Im{ω} = 0) associated with the smallest Rayleigh numbers are positioned in the vicinity of the tangent cylinder s = 1, we approximate Q (see (2.3e), but cf. (2.15)) and a2 (see (2.20b)) by

Q (r) = [(1 + z2 )

−3/2

−3/2

− (1 + ε 2 )

]/[1 − (1 + ε 2 )

−3/2

],

a2 = m2 . (2.23a,b)

Accordingly, we solved numerically the system (2.21a,b) subject to the boundary conditions (2.22a,b) using standard spectral methods. Solutions were written as series of suitably transformed Chebyshev polynomials and, given m, Newton iteration was applied to the values of R (m) and ω(m) so that all of the boundary conditions were fulfilled. The radius ratio (2.2d) was taken to be η = 0.35 for consistency with the choices made by Jones et al. (2000) and by Dormy et al. (2004). The real Rayleigh number R (m) and frequency ω(m) for the tangent cylinder s = 1 p −2 (zo = η − 1 ≈ 2.676) neutral modes were determined for various values of the wavenumber m and stratification parameter ε and results for various Prandtl numbers P are illustrated in figures 2 and 4. The onset of instability is determined by the minimum value of R over m. That is given by the stationary value Rc = R (mc ), at which dR /dm = 0; the corresponding frequency is ωc = ω(mc ). For various values of ε we plot the eigenfunctions uc = imψ(z) and wc = w(z) versus z/zo of the critical modes at m = mc , normalized by ψ(0) = −i/m in figures 3 and 5. The values of zn /zo , at which the switch (Qn = 0) between stable and unstable stratification occurs, are listed in table 1 (also table 2). The corresponding tangent cylinder geometry for the case ε = 1 is drawn to scale in figure 1.

387

Localized rotating convection (a) 4

(b) 0.2

3

0.1

2

0 1 – 0.1 0

0.5

1.0

(c) 4

(d )

0

0.5

0

0.5

1.0

5 4

3

3 2

2 1

1

0 0

0.5

1.0

m

1.0

m

F IGURE 2. Numerical solutions of (2.21a,b) subject to the boundary conditions (2.22a,b). The left-hand plots show the forms of R (m) (real) for the choices ε = 1 (line with • marks), ε = 10−0.5 (line), 10−1 (dash), 10−1.5 (dot-dash) and 10−2 (dots). Also included is the asymptotic result Rs0 for ε ↓ 0 given by (4.3c) and denoted by the line with × marks. The minimum (mc , Rc ) on each curve identifies the critical mode. The right-hand plots show ε−1 mω using the same choices of ε with the asymptotic form mωs1 (real) defined by (4.6a). The critical frequency ωc is identified by the value of ε−1 mω at m = mc . The Prandtl number is P = 5 in (a, b) and P = 1 in (c, d).

We begin by considering the case P = 5, because at that largish value of P the Rayleigh number R plotted versus m in figure 2(a) exhibits a relatively mild dependence on the value of ε chosen. Indeed for m greater than roughly unity, the plots for the various values of ε 6 0.316 · · · are virtually indistinguishable. As m is decreased below unity, the sensitivity of the results on the value of ε increases and substantial variations become visible; although in the limit ε ↓ 0, the dependence of R (m) evidently assumes a well-defined form. In the same small-ε limit, the frequency ω approaches zero and so, in order to visualize this quantity, we plot ε −1 mω versus m in figure 2(b). Once more a well-defined frequency, scaled with ε −1 , appears to emerge as ε ↓ 0. In figure 3 we plot the critical mode eigenfunctions uc and wc at m = mc , for various values of ε with the eigenvalues Rc and ε−1 mc ωc listed in table 1. Since, for all m, the imaginary part of u and the real part of w both tend to zero as ε ↓ 0, we plot in figure 3(c,d), for the particular critical cases, the scaled values Im{uc }/ε and Re{wc }/ε, respectively. All plots in figure 3(a–d) take on a self-similar form relative to the stretched variable z/ε as ε ↓ 0, which is evidenced by a comparison of the ε = 10−1 and 10−1.5 dashed and dot-dashed curves, respectively. The corresponding plots of R versus m from computations for P = 1, are illustrated in figure 2(c,d). They show similar features to the P = 5 case, although the sensitivity to ε at small m becomes more evident.

388

A. P. Bassom, A. M. Soward and S. V. Starchenko

(a)

(b) 1.0

Im{wc}

Re{u c}

1.0

0.5

0.5

0

0 0

(c)

0.5

0

1.0

(d)

0

1.0

0 –0.5

Re{wc} ε

Im{uc} ε

–0.2

0.5

–0.4 –0.6

–1.0 –1.5

–0.8 –2.0

–1.0

–2.5 0

0.5

1.0

0

0.5

z z0

1.0

z z0

F IGURE 3. The P = 5 critical mode eigenfunctions uc , wc with the eigenvalue data listed in table 1, subject to the normalization uc (0) = 1, plotted versus z/zo , for the choices ε = 1 (line with • marks), ε = 10−0.25 (line with + marks), ε = 10−0.5 (line), ε = 10−1 (dash), ε = 10−1.5 (dot-dash): (a) Re{uc }, (b) Im{wc }, (c) Im{uc }/ε, (d) Re{wc }/ε.

ε zn /zo mc Rc ε −1 mc ωc

1

10−0.25

10−0.5

10−1

10−1.5

0

0.3736 0.44 0.6344 0.020

0.2101 0.50 1.1425 0.050

0.1182 0.66 1.8559 0.093

0.0374 0.808 2.2606 0.076

0.0118 0.822 2.2984 0.072

0 0.8238 2.3025 0.07143

TABLE 1. The P = 5 critical mode data for the eigenfunctions plotted in figure 3. The right-hand ε = 0 column gives the asymptotic ε ↓ 0 values of msc0 , Rsc0 and msc0 ωsc1 defined by (4.5d,c) and (4.8) respectively.

When P < 1 a different solution structure emerges from the results in figure 4. That a new structure exists is readily concluded from the plots of mω for P = 0.8, 0.65, 0.2 in figure 4(b,d,f ), respectively. As the limit ε ↓ 0 is approached, we see that, for some mt (P) (see (4.14c)), mω tends to a well-defined non-zero limit when m < mt , but to zero when m > mt , as before in the P > 1 cases illustrated in figure 2(b,d). This numerical evidence suggests there are at least two possible modes distinguished by the behaviour of the frequency as ε ↓ 0. Henceforth, we call the first family of lowfrequency waves with ω = O(ε) illustrated in figure 2 for all m (P > 1) and figure 4 for m > mt (P) (P < 1) as slow modes. The second family of relatively high-frequency waves with ω = O(1) appearing in figure 4 for m < mt (P) will be referred to as fast modes.

389

Localized rotating convection (a) 4

(b)

0.3 0.2

3

0.1 2 0 1

0

–0.1 0.5

–0.2

1.0

0

0.5

1.0

0

0.5

1.0

(d) 0.4

(c) 4 3

0.2

2 0 1 –0.2 0

0.5

1.0

(e) 4

(f ) 0.4

3

0.2

2

0 1 –0.2 0

1.0

0.5

m

0

0.5

1.0

m

F IGURE 4. As in figure 2 but for the Prandtl number P = 0.8 (mt ≈ 0.6233) in (a,b), P = 0.65 (mt ≈ 0.6997) in (c,d) and P = 0.2 (mt ≈ 0.8736) in (e,f ). Note, however, that the right-hand plots are of mω rather than the scaled ε−1 mω. Moreover, whereas the slow mode asymptotic values Rs0 and mω = 0 for m > mt continue to be denoted by the lines with × marks, the fast mode asymptotic values R0 and mω0 6= 0 for m < mt , given by the results of § 5.1.1, are denoted by the lines with ◦ marks.

These terms ‘slow’ and ‘fast’ refer to the relative magnitudes of the angular phase speeds of the waves and we discuss them later in §§ 4 and 5, respectively. The analyses of §§ 4.1 and 5.1 show that, in the limit ε ↓ 0, the fast mode branch originates from a frequency pitchfork of the slow mode branch. Upon closer inspection, the asymptotic analysis of § 5.1.3 reveals that at small but finite ε the frequency pitchfork is imperfect, which for the values P = 0.8 and 0.65 used in figure 4(a–d) is illustrated in figure 8 later in this article. Those plots indicate the existence of a second solution branch not apparent from the graphs in figure 4, which at given m only capture the mode with the lowest Rayleigh number. In fact, the second solution branch consists of

390

A. P. Bassom, A. M. Soward and S. V. Starchenko (a)

(b) 0.4

1.0

0.3

Im{wc}

Re{u c}

0.8 0.6 0.4 0.2

0.2 0.1 0

0 0

0.5

–0.1

1.0

(c)

(d) 0.1

0

0.5

1.0

0

Im{uc}

Re{wc}

–0.2

0

–0.4 –0.6 –0.8 –1.0

0

0.5

1.0

0

0.5

z z0

1.0

z z0

F IGURE 5. The P = 0.2 critical mode eigenfunctions with the eigenvalue data listed in table 2, as in figure 3: (a) Re{uc }, (b) Im{wc }, (c) Im{uc }, (d) Re{wc }.

ε zn /zo mc Rc mc ωc

1

10−0.25

10−0.5

10−1

10−1.5

0

0.3736 0.24 0.0779 0.111

0.2101 0.26 0.2330 0.129

0.1182 0.34 0.4061 0.24

0.0374 0.404 0.6202 0.434

0.0118 0.430 0.6994 0.472

0 0.4395 0.7379 0.499

TABLE 2. The P = 0.2 critical mode data for the eigenfunctions plotted in figure 5. The right-hand ε = 0 column gives the asymptotic ε ↓ 0 values of mfc0 , Rfc0 and mfc0 ωfc0 obtained from the method outlined in § 5.1.2.

the slow mode branch for m < mt (not plotted in figure 4) together with a fast mode branch almost identical to that plotted but with a frequency of opposite sign. In this way the frequency plots of figure 4(b,d,f ) would exhibit the pitchfork form visible in figure 8(b). An examination of figures 2(a,c) and 4(a,c,e) suggests that as ε ↓ 0 the slow mode limit for R is independent of the value of P. In contrast, the fast mode limit changes and the value m = mt (P), defined by (4.14c), at which it bifurcates from the slow mode increases with P. When P = 0.8, figure 4(a) shows that the fast mode branch bifurcates from the left side of the slow mode branch and that the most dangerous mode (that is, the m that corresponds to the smallest R ) is actually of slow type. As the value of P is decreased, so the fast mode develops a local minimum whose value, determined in the analysis of § 5.1.2 below, coincides with the minimum of the slow mode at P = Pc ≈ 0.65; the case illustrated in figure 4(c). On decreasing P yet

Localized rotating convection

391

further, the minimum of the fast mode continues to decrease and the identity of the preferred instability switches from slow to fast. Indeed for sufficiently small P the bifurcation point mt moves to the right side of the slow mode branch, with the fast mode minimum lying well below the slow mode minimum as figure 4(e) for the case P = 0.2 indicates; the corresponding critical mode eigenfunctions for various values of ε are illustrated in figure 5 with the eigenvalues Rc and mc ωc listed in table 2, and as ε ↓ 0 these are clearly of fast mode type. To sum up then, the numerical results for ε ↓ 0 illustrated in the pairs of figures 2 and 3 together with figures 4 and 5 reveal many properties of the solutions, the most significant of which are identified by the frequency plots. Those clearly indicate distinct slow and fast mode structures, whose relative importance depends on the size of P. The resolution of the small-ε characteristics by asymptotic methods is the subject of the remainder of the paper. 3. Asymptotic theory for ε  1 In order to understand the nature of the two classes of waves identified in § 2.3, we now cast the governing equations (2.21a,b) with Q (r) given by (2.23a) applicable on the tangent cylinder s = 1 in a form more suited to the limit ε ↓ 0, for which

Q (r) = 1 − ζ 2 + O(ε2 ),

(3.1)

where ζ = z/ε. Asymptotic solutions are sought of the form ˆ ) exp(−βζ 2 /2), ψ(ζ ) = ψ(ζ

w(ζ ) = ζ w(ζ ˆ ) exp(−βζ 2 /2),

(3.2a,b)

ˆ ), w(ζ where ψ(ζ ˆ ) have polynomial-like behaviours and β is restricted by the condition Re{β} > 0

(3.2c)

to ensure boundedness as ζ ↑ ∞ but is otherwise to be chosen at our convenience. The substitution of (3.2a,b) into (2.21a,b) shows that ψˆ and w ˆ satisfy

ζ

1 dψˆ = (β − iεmC )ψˆ − B w, ˆ ζ dζ

(3.3a)

dw ˆ ˆ = [−1 + (β + iεmC )ζ 2 ]w ˆ − F ψ, dζ

(3.3b)

for which substitution of (3.1) and a2 = m2 into (2.21c–e) determines

C=

R (1 − ζ 2 ) + O(ε 2 ), B = −m2 + iω + O(ε 2 ), m2 − iPω F = m2 C − m2 (m2 − iω) + O(ε2 ).

(3.3c,d) (3.3e)

The key characteristics of the solutions of the system (3.3a,b) describing the neutral modes Im{ω} = 0 are the real values of the Rayleigh number R , the frequency ω and the complex√value of β. Here Re{β} is important as it measures the e-folding ζ distance Λ = 1/ (Re{β/2}) (see (7.3b) below) of the convective columns. In addition to that axial decay, the azimuthal tilt of the cell boundaries is determined by the phase surfaces Mφ = Im{β/2}ζ 2 + constant (see (7.2)), because the leading order approximations to ψˆ and w ˆ are just complex constants. We consider asymptotic expansions based on

R = R0 + εR1 + O(ε2 ),

ω = ω0 + εω1 + O(ε2 )

(3.4a,b)

392

A. P. Bassom, A. M. Soward and S. V. Starchenko

with β = O(1). Correct to O(ε), we write (3.3c–e) as

R0 (1 − ζ 2 ) + O(ε), B = −m2 + iω0 + iεω1 + O(ε 2 ), (3.5a,b) m2 − iPω0 m2 R0 (1 − ζ 2 ) − m2 (m2 − iω0 ) F= 2 m − iPω0    m2 R1 Pm2 R0 2 2 2 +ε (1 − ζ ) + iω1 (1 − ζ ) + m m2 − iPω0 (m2 − iPω0 )2

C=

+ O(ε2 )

(3.5c)

and seek solutions of (3.3a,b) in the form ψˆ = ψˆ 0 + iεψˆ 1 βζ 2 + O(ε2 ),

w ˆ =w ˆ 0 + iε(w ˆ 10 + w ˆ 12 βζ 2 ) + O(ε 2 ).

(3.6a,b)

We consider the general nature of the zeroth- and first-order expansions in the following sections, before applying them explicitly to the slow and fast modes in §§ 4 and 5. 3.1. Zeroth order

On substitution of ψˆ ≈ ψˆ 0 and w ˆ ≈w ˆ 0 into (3.3a,b) the lowest-order terms determine   m2 R0 2 2 − m (m − iω ) ψˆ 0 , (3.7a,b) −β ψˆ 0 = (m2 − iω0 )wˆ 0 , −wˆ 0 = 0 m2 − iPω0 −β w ˆ0 = which together yield

m2 R0 ψˆ 0 , (m2 − iPω0 )

(3.7c)

m2 − iω0 = β 2. (3.8) m2 − iPω0 We refer to the first and second equalities as (3.8)1 and (3.8)2 respectively. Together the real and imaginary parts of these two complex equations provide four real equations for R0 , ω0 , Re{β}(> 0) and Im{β} at given m and P. 2

β + m2 (m2 − iω0 ) = m2 R0

3.2. First order At O(ε), we take advantage of the two alternative expressions for R0 given by (3.8). We substitute into (3.5a,c) the value given by (3.8)2 , when R0 is multiplied by ζ 2 , and the value given by (3.8)1 , when multiplied by ζ 0 (= 1). This determines the representations

1 − βζ 2 + β −1 m2 (m2 − iω0 ), m2 − iω0   1 − βζ 2 R1 1 − βζ 2 −1 −1 2 2 β F ≈ 2 +ε + β m (m − iω0 ) m − iω0 R0 m2 − iω0   iεω1 P(1 − βζ 2 ) −1 2 2 + 2 + β m ((1 + P)m − 2iPω0 ) . m − iPω0 m2 − iω0 β −1 m2 C ≈

(3.9a)

(3.9b)

We substitute (3.6a,b), (3.9a,b) and (3.5b) into (3.3a,b) and equate coefficients of the various powers of ζ . With the help of ψˆ 0 = −β −1 (m2 − iω0 )wˆ 0 (from (3.7a)), the O(ε)

393

Localized rotating convection 0

2

0

2

terms proportional to ζ and ζ in (3.3a), and ζ and ζ in (3.3b) give, respectively,   mψˆ 1 mw mω1 1 ˆ 10 −1 2 2 2 = 2 − + β m (m − iω ) (= a), (3.10a) + 0 w ˆ0 m − iω0 m2 − iω0 ψˆ 0 mψˆ 1 mw ˆ 12 1 − =− 2 (= b), (3.10b) ˆ w ˆ m − iω0 ψ0 0 mω1 mw ˆ 10 = 2 [P + β −1 m2 (m2 − iω0 )((1 + P)m2 − 2iPω0 )] w ˆ0 m − iPω0 imR1 2 − [1 + β −1 m2 (m2 − iω0 ) ] (= c), (3.10c) R0 ˆ 10 − 3w ˆ 12 ) Pmω1 imR1 mψˆ 1 m(w + = 2 − ˆ w ˆ0 m − iPω0 R0 ψ0   1 −1 2 2 − + β m (m − iω0 ) (= d). (3.10d) m2 − iω0

(Note that (3.10b) is also recovered from the O(ε) terms proportional to ζ 4 in (3.3b).) The set of four equations (3.10a–d) only have a solution when their right-hand sides a, b, c, d satisfy a − 3b − 2c + d = 0 implying   mω1 (1 − P)m2 −1 2 2 2 − 2β m ((1 + P)m − 2iPω0 )(m − iω0 ) m2 − iPω0 m2 − iω0 imR1 2 [1 + 2β −1 m2 (m2 − iω0 ) ] + R0 1 = − 2 + 2β −1 m2 (m2 − iω0 ). (3.11) m − iω0 The real and imaginary parts of this equation determine the values of R1 and ω1 .

4. Slow modes The numerical evidence outlined in § 2.3 suggests that there are slow modes with ω0 = 0. Under this simplifying assumption, the zeroth-order result (3.8) reduces to (4.3a) below, which implies that β is real and so Im{β} = 0. Consequently, at first order, the right-hand side of (3.11) is real. Then (3.11) can only be solved when R1 = 0, leading to the formula (4.6a) below. Thus, the more precise reduced forms of the expansions (3.4a,b) are now

R = Rs ≡ Rs0 + ε2 Rs2 + O(ε4 ),

ω = ωs ≡ εωs1 + O(ε 3 ),

(4.1a,b)

ωsc = εωsc1 + O(ε 3 ),

(4.2a,b)

where we identify slow mode values by the subscript s and have anticipated that factors of ε 2 generate the higher-order terms in the expansions. We refer to the mode that minimizes the Rayleigh number as the critical mode and the critical values will be distinguished by the additional subscript c . Then, for example, the critical values Rs = Rsc , ωs = ωsc and ms = msc have expansions

Rsc = Rsc0 + ε2 Rsc2 + O(ε4 ),

msc = msc0 + O(ε 2 ).

(4.2c)

394

A. P. Bassom, A. M. Soward and S. V. Starchenko

4.1. Zeroth order With ωs0 = 0 the pair of complex equations (3.8) simplifies to the pair of real equations

βs + m6 = m2 Rs0 = βs2

subject to βs > 0 (see (3.2c)) with solutions p βs = 12 (1 + 1 + 4m6 ) (> 1), p Rs0 = 21 (1 + 2m6 + 1 + 4m6 )/m2 .

(4.3a) (4.3b) (4.3c)

The relationship (4.3c) for Rs0 (m) is plotted as the asymptotic lines in figure 2(a,c), which show excellent agreement with the numerical solutions in the limit ε ↓ 0. Of particular significance is the verification that this leading-order solution is independent of the Prandtl number P. Differentiation of (4.3a) with respect to m yields the two results m dβs 6(βs − 1) = (> 0), βs dm 2βs − 1

m dRs0 2(4βs − 5) = , Rs0 dm 2βs − 1

(4.4a,b)

from which it is straightforward to minimize Rs0 over m. The vanishing of the righthand side of (4.4b) and use of βs−1 m6 = βs − 1 from (4.3a) immediately determines the lowest-order critical values βsc = 5/4,

which on further use of (4.3a) yields

Rsc0 = 5m4sc0 ,

βsc−1 m6sc0 = 1/4,

(4.5a,b)

m2sc0 = (5/16)1/3 .

(4.5c,d)

4.2. First order With ωs0 = 0, Rs1 = 0 and βs real, the relation (3.11) for ωs1 reduces to

mωs1 =

1 βs − βt‡ 1 − 2βs−1 m6 = − (P − 1) + 2(P + 1)βs−1 m6 1 + P βs − βt

on use of βs−1 m6 = βs − 1, where βt (P) =

3+P 2(1 + P)

3 and βt‡ ≡ βt (0) = ; 2

(4.6a)

(4.6b,c)

hereafter the superscript ‡ is used to denote P = 0 values. Of course, since βs (m) is given by (4.3b), the relation (4.6a) yields mωs1 as a function of m, and it is this function that is plotted as the asymptotic limit in figure 2(b,d) as ε ↓ 0. The key values of m are mt (when it exists) at which βs = βt implying that mωs1 is infinite, and m‡t = (3/4)1/6 = 0.95318 . . . at which βs = βt‡ = 3/2 implying that mωs1 vanishes. Significantly the formula (4.6a) exhibits the limiting behaviours ( 0 as P ↑ ∞, mωs1 → (4.7) ‡ −1 (m 6= mt ) as P ↓ 0. At critical, (4.6a) together with the result βsc−1 m6sc0 = 1/4 (see (4.5b)) determines msc0 ωsc1 =

1 ≶0 3P − 1

for P ≶ Psct ≡ 1/3

(4.8)

395

Localized rotating convection

consistent with the limiting forms (4.7). Then use of (3.10a–d) with (4.5a,b) determines m3sc0 ψˆ sc1 1 5P − 2 , =− ˆ 2 3P − 1 ψsc0

where

m3sc0 w ˆ sc10 1 5P + 1 = , w ˆ sc0 4 3P − 1

(4.9a)

m3sc0 w ˆ sc12 1 P = , w ˆ sc0 2 3P − 1

wˆ sc0 = −4m4sc0 ψˆ sc0 .

(4.9b,c)

(4.9d)

The small-ε critical eigenfunction plots (ε 6 10 ) of Re{uc } (uc = imc ψc ) and Im{wc } at m = mc in figure 3(a,b) conform to the latter result (4.9d), while the corresponding plots of Im{uc }/ε and Re{wc }/ε agree with the former results (4.9a–c). Likewise at mc the eigenvalues Rc and ε−1 mc ωc listed in table 1 approach the asymptotic values msc0 , Rsc0 and msc0 ωsc1 predicted by (4.5c,d) and (4.8). The formula (4.6a) reveals some interesting features of the frequency that depend on the following properties. First, since dβs /dm > 0 (see (4.4a)), m and βs (> 1 see (4.3b)) increase together in concert. Second, since dβt /dP = −1/ (1 + P)2 < 0 (from (4.6b)), βt decreases monotonically with increasing P. The key values −1

βt (0) ≡ βt‡ = 3/2,

show that βt lies in the range So, when

βt‡

βt (1) = 1

> βt > 1/2.

and

βt (∞) = 1/2

P > 1,

(4.10a,b,c)

(4.11)

(4.10b) implies that βt 6 βt (1) = 1, with the consequence that βs > 1 > βt . It follows from (4.6a) that mωs1 ≷ 0 when βs ≶ βt‡ for m ≶ m‡t with limiting values ( 1/(P − 1) when βs ↓ 1 as m ↓ 0, mωs1 → (4.12) −1/(P + 1) when βs ↑ ∞ as m ↑ ∞. These various behaviours of mωs1 , as well as the critical value msc0 ωsc1 = 1/(3P − 1) (see (4.8)) at m = msc0 = (5/16)1/6 = 0.82377 . . . (< m‡t ) (see (4.5d)), agree well with the direct numerical solutions for the case P = 5 portrayed in figure 2(b). The agreement is less impressive in figure 2(d) when P = 1 but this can be explained by noting that mωs1 ↑ ∞ as m ↓ 0 (see (4.12)). 4.3. The role of the parameter βt > 1, when P < 1 The above conclusions change somewhat when

P < 1,

(4.13)

for then βt > 1 with the added complication that mωs1 also changes sign at m = mt , where βs = βt . Use of the formula (4.3a) determines Rt ≡ Rs0 (mt ) = βt2 /m2t and m6t = βt (βt − 1). Since βt = (3 + P)/[2(1 + P)] (see (4.6b)), it is convenient to write y2t ≡ 2(βt − 1) = (1 − P)/(1 + P),

(4.14a)

as in (5.24b) below, and so obtain the compact expressions 1/3

Rt ≡ (2βt5 /y2t ) ,

m2t ≡ (y2t βt /2)

1/3

.

(4.14b,c)

396

A. P. Bassom, A. M. Soward and S. V. Starchenko

We note that like βt both mt and yt are decreasing functions of P. Thus, mt lies in the range 0 < mt 6 m‡t (0 < y2t 6 1), where 1 < βt 6 βt‡ . In fact m2t starts with m2t ↓ 0 as P ↑ 1, reaches the critical value m2t = m2sc0 = (5/16)1/3 , where βt = βsc = 5/4 1/3 (see (4.5a,d)), when P = Psct ≡ 1/3 (y2t = 1/2) and finishes with m2t ↑ m‡2 t = (3/4) as P ↓ 0. In view of the inequalities just described we note that we still have mωs1 ≶ 0 when βs ≶ βt‡ for m ≶ m‡t , provided that βs > βt (i.e. m > mt ). However, the sign reverses across m = mt so that mωs1 < 0 when βs < βt (i.e. m < mt ). Nevertheless, this latter property is of no real interest because the preferred mode of convection at given m is of fast rather than slow type whenever m < mt (see figure 6a below). Indeed figure 4(b,d,f ) for P = 0.8, 0.65, 0.2, respectively, illustrate the slow mode behaviour described for m ' mt with the vanishing of ω at m ≈ m‡t . Evidently the frequency diverges as m ↓ mt signalling an imperfect frequency pitchfork to relatively large frequency modes, which are the fast modes noted in the caption of figure 4 and analysed in § 5. In order to capture the smooth transition from the slow to the fast modes in the immediate neighbourhood of m = mt , we require the Taylor expansions of βs − βt and Rs0 − Rt . To that end, we note that the derivatives (4.4a,b) take the values     m dβs 3 m dRs0 = (1 − P), = 1 − 3P, (4.15a,b) βs dm t 2 Rs0 dm t which combined with (4.6a) determine   m − mt 1 4P mωs1 + ≈ , mt 1+P 3(1 − P2 )(3 + P) Rs0 − Rt m − mt ≈ (1 − 3P) Rt mt

(4.16a) (4.16b)

for |m − mt |  1. Despite the fact that the term 1/(1 + P) on the left-hand side (4.16a) is negligible as m → mt , we have retained it to emphasize that its neglect is only valid for |m − mt |  P, when P  1. 5. Fast modes: 1 > P = O(1) The second family of solutions identified in figure 4 was seen to possess an O(1) frequency as ε ↓ 0. To describe these the entire development of § 3 is needed. We denote the critical values of these fast modes by the subscripts fc . Then, for example, R = Rfc , ω = ωfc and m = mfc have the expansions

Rfc = Rfc0 + εRfc1 + O(ε2 ),

ωfc = ωfc0 + εωfc1 + O(ε2 ),

mfc = mfc0 + εmfc1 + O(ε 2 ).

(5.1a,b) (5.1c)

5.1. Zeroth order In this section we analyse the formula (3.8) for real R0 , ω0 and complex β at given m. We begin by casting (3.8)1 and (3.8)2 in the suggestive forms  2 2 β (m2 − iω0 ) R0 m2 − iω0 β R0 m2 − iω0 + = , = , (5.2a,b) m6 m4 m4 m2 − iPω0 m6 m10 m2 − iPω0

397

Localized rotating convection

where the complex number β is restricted by Re{β} > 0. Their reduction is aided by writing m2 − iω0 ≡ m2 (Rν /P)e−iθν ,

m2 − iPω0 ≡ m2 Rκ e−iθκ .

(5.3a,b)

Consideration of the product (m2 ± iω0 )(m2 − iPω0 ) leads to the useful identities

together with

1 + P i(θν −θκ ) 2P Rκ e , = e−i(θν +θκ ) + 1 − P Rν 1−P  
Pω2 1 R2ν + PR2κ = P 1 + 40 = Rν Rκ cos(θν − θκ ), 1+P m R2κ − R2ν = 1 − P2 ,

R2ν − P2 R2κ = (1 − P2 )

P2 ω02 . m4

The substitution of (5.3a,b) into (5.2b) yields s β Rν D e−i(θν −θκ )/2 , where R0 = m10 D 2 , = 6 m PRκ

(5.4a) (5.4b)

(5.4c,d)

(5.5a,b)

for some positive D (since Re{β} > 0). In turn, substitution of this value of β into (5.2a) determines (m2 − iPω0 )

1/2

m4 D + (m2 − iω0 )(m2 − iPω0 ) (m2 − iω0 )1/2   m4 Rν Rκ P3/2 D i(θν −θκ )/2 2P Rκ 1 + P i(θν −θκ ) = , e + − e P 1 − P Rν 1−P Rν3/2 Rκ1/2

R0 =

(5.6)

where, to derive (5.6)2 , we have made use of (5.4a). The vanishing of the imaginary part of the right-hand side of (5.6)2 then leads to P3/2 D

1+P cos [(θν − θκ )/2] 1−P

(5.7)

implying (π/2 >) |θν | > |θκ |.

(5.8a,b)

Rν3/2 Rκ1/2

=2

provided that θν 6= θκ . Since D > 0, it follows that the fast modes under consideration only exist when P<1

This is in accordance with the earlier numerical evidence, which suggested that only the slow mode solution is possible when P > 1. 5.1.1. A parametric representation in terms of x = Pω0 /m2 The result (5.7) for D provides the key to further reductions. For on substitution of this value into the two formulae (5.5b) and (5.6), we obtain a coupled pair of equations for R0 and m that ultimately leads to their parametric representation in terms of x = Pω0 /m2 (see (5.13c) below). The first substitution into (5.5b) determines

(1 − P)2 P3 R0 = (1 + P)R3ν Rκ [1 + cos(θν − θκ )] = R2ν R2κ ΣΣP , 1 + P 2m10 where, to obtain (5.9a)2 in which Σ ≡ (Rν + Rκ )/Rκ

and ΣP ≡ (Rν + PRκ )/Rκ ,

(5.9a)

(5.9b,c)

398

A. P. Bassom, A. M. Soward and S. V. Starchenko

we have made additional use of (5.4b). On making the second substitution into the right-hand side of (5.6)2 , the remaining real equality gives PR0 = R2κ (PΣ + ΣP ). m4 Together (5.9a,d) determine the more useful relations (1 − P)

P R2κ (PΣ + ΣP )2 m R0 = , 2(1 + P) R2ν ΣΣP 2

m6 =

(5.9d)

P2 (1 − P) PΣ + ΣP . 2(1 + P) R2ν ΣΣP

Finally substitution of (5.7) into (5.5a) gives   1 − P P2 β = R2ν [1 + e−i(θν −θκ ) ]. 1 + P m6

(5.10a,b)

(5.11)

The substitution of m as given by (5.10b) into the expression (5.10a) for m2 R0 leads to the parametric representations 1/3 5/3  1/3 1/6  G P G P and m = (5.12a,b) R0 = 2/3 1+P F 1+P F 1/6

for R0 and m in terms of the functions F(P; x) ≡ where from (5.3a,b)

2R2ν R2κ ΣΣP (1 − P2 )2

Rν (P; x) ≡

p

,

P2 + x 2

are functions of the single parameter

G(P; x) ≡

R2κ (PΣ + ΣP ) , 1 − P2

and Rκ (x) ≡

p

1 + x2

x ≡ Pω0 /m2 .

Equivalently, ω0 ≡ xm2 /P or

mω0 = x (G/F)1/2 /(1 + P),

(5.12c,d)

(5.13a,b) (5.13c) (5.13d,e)

so that R0 , m and ω0 are determined as functions of the Prandtl number P and the parameter x by (5.12a,b) and (5.13d), respectively. Significantly, however, whereas R0 and m only depend on the magnitude of x through x2 , ω0 can take positive or negative values depending on the sign of x. In identifying the critical Rayleigh number this sign distinction becomes a crucial issue and so we introduce the notation ω0± ≡ ±|ω0 | = (sgn x)|ω0 |.

(5.13f )

Finally, we note that further use of (5.10b) and (5.11) determines, in the notation of (5.13f ),    1 1 P  β± = β = + 1 + e−i(θν −θκ ) for x ≷ 0. (5.14) 2 Σ ΣP From (5.13a,b) we note that P 6 Rν /Rκ < 1 and so, on the range 0 6 P < 1, use of the definitions (5.9b,c) shows that 1/2 < Σ −1 + PΣP−1 < 3/2. With this simple bound and recalling that |θν − θκ | < π/2, we deduce from (5.14) that 1/4 < Re{β± } < 3/2

and

0 < ∓Im{β± } < 3/4,

where Re{β+ } = Re{β− } and Im{β+ } = −Im{β− }.

(5.15a,b)

399

Localized rotating convection

The limiting behaviours as x → ±∞ may be determined upon noting that Rν ≈ Rκ ≈ |x|, Σ ≈ 2, ΣP ≈ 1 + P and (θν − θκ ) ≈ 0 (both θν and θκ ≈ ±π/2). Then (5.10a,b) give   1 + 3P 2 m3 x [(1 − P)(1 + 3P)]1/2 m2 R0 ≈ P , ≈± , (5.16a,b) 2(1 + P) P 2(1 + P)

while (5.13d), (5.14) yield

m3 x 1 + 3P , β≈ , (5.17a,b) P 2(1 + P) all of which pertain to the limit m ↓ 0 by (5.16b). As the magnitude of x is reduced, the two solutions linked to x ≷ 0 merge as x → 0. At x = 0 they provide a frequency pitchfork of the slow mode solution anticipated in § 4.3. There the values Rν = P, Rκ = 1, Σ = 1 + P, ΣP = 2P and θν = θκ = 0 lead as above to the limiting values   3+P 2 [(1 − P)(3 + P)]1/2 2 2 m R0 = mt Rt ≡ , (5.18a,b) , m3 = m3t ≡ 2(1 + P) 2(1 + P) 3+P mω0 = 0, β = βt ≡ (5.19a,b) 2(1 + P) mω0 ≈

(see (4.6b) and (4.14a–c)). For small x(6= 0), expansions taken to O(x2 ) yield

R0 − Rt x2 ≈ −Υ 2 , Rt P

m − mt x2 ≈ −Γ 2 , mt P

Υ (P) =

Γ (P) =

5 + 3P − 10P2 − 2P3 , 3(3 + P)

(1 + P)(13 − 4P − P2 ) 24(3 + P)

(> 0).

(5.20a,b)

(5.20c,d)

Since x2 /P2 ≈ ω02 /m4t , (5.20a,c) imply that m − mt ω02 ≈− , 4 mt Γ mt

R0 − Rt m − mt ≈ −Υ Rt Γ mt

(5.21a,b)

for 0 6 mt − m  1. Used in conjunction with of the asymptotic relation (4.16b) for Rs0 − Rt , we can write (5.21b) in the alternative form

R0 − Rs0 ω2 ≈ −Υs 40 , Rt mt

(5.22a)

where Υs (P) = Υ − (1 − 3P)Γ = (1 − P)(3 + 8P + P2 )/8

(> 0).

(5.22b)

In the next section we analyse the solutions elsewhere more carefully and show that, because m is a monotonically decreasing function of |x| (x−1 dm/dx < 0, see (5.26b,c)), the fast mode solutions are restricted to the range 0 < m 6 mt . As m ↑ mt , (5.21a) shows that the fast mode merges (ω0 → 0, see (5.21a)) with the slow mode from below (R0 ↑ Rs0 , see (5.22a,b)). Plots of R0 versus m obtained from our parametric representations are shown in figure 6(a), as well as the slow mode solution Rs0 . The various bifurcations of the fast mode solution R0 from it for different values of P (< 1) are clearly visible, including that from the steady mode

400

A. P. Bassom, A. M. Soward and S. V. Starchenko

(a) 3

(b) 0.8 0.6

2 0.4 1 0.2

0

0.2

0.4

m

0.6

0.8

1.0

0

0.2

0.4

m

0.6

0.8

1.0

F IGURE 6. The marginal solutions: (a) R0 versus m. The slow solution Rs0 (see (4.3c)) is illustrated by the continuous curve (upper right). The fast solutions (see (5.12a,b)), which bifurcate from the steady solutions, are illustrated from left to right by the cases P = Pfct ≈ 0.8 (see (5.28b)) (short-dashed), P = Pc ≈ 0.65 (dotted), P = Psct ≈ 0.33 (chain), P = 0.1 (dots-dashed), P = 0.01 (long-dashed), P = 0.001 (dotted) and P = 0 (the R0‡ -curve defined by (6.4a,b); continuous). The curve Rfc0 versus mfc0 , through the minima (determined by the zero |y| = |yfc0 | for the vanishing of the right-hand side of (5.26a)) as P varies over the ‡ range 0 < P < Pfct , is shown continuous. (b) mω0 (= mω0+ only) (see (5.13c)), mω0‡ = mω0+ (see (6.3a) and (6.4b)) with the line style of (a). (a) 1.6

(b) 0

1.2

Im{ }

Re{ }

1.4

1.0 0.8

–0.2

Im{

}

–0.4 Re{

0.6

Re{ }

}

0.4 0

0.2

0.4

0.6

m

0.8

1.0

–0.6 0

0.2

0.4

0.6

Im{ } 1.0 0.8

m

F IGURE 7. The marginal solutions: β (= β+ only) (see (5.14)), β ‡ (= β+‡ ) (see (6.4c)). (a) Re{β}, (b) Im{β} versus m with the line styles of figure 6.

critical values m = msc0 , R0 = Rsc0 , which occurs when P = Psct ≡ 1/3 (see § 4.3). The corresponding plots of mω0+ for the same values of P used in figure 6(a) are shown on figure 6(b); each curve exhibits the m = 0 value (5.17a) with m3 x/P given by (5.16b). Likewise we plot the real and imaginary parts of β+ on figures 7(a,b) finding them to lie inside the bounds identified in (5.15a,b). Their end point values (5.17b) and (5.19b) are both attained. 5.1.2. An alternative parametric representation Although the introduction of the parameter x enables a direct parametric solution of the dispersion relation (3.8), or its equivalent (5.2a,b), it is cumbersome to work with because the key parameters Rν and Rκ defined by (5.13a,b) involve square roots. Nevertheless the relation (5.4c) between them allows us to take this parametric form of solution further and express Rν and Rκ in terms of a new parameter in the following

Localized rotating convection

401

way. We set Rν =

p

1 − P2

1 − y2 , 2|y|

Rκ =

p 1 + y2 1 − P2 , 2|y|

(5.23a,b)

which we substitute into (5.9b,c) to obtain Σ=

2 1 + y2

ΣP =

(1 + P) − (1 − P)y2 . 1 + y2

(5.23c,d)

Furthermore substitution of the formulae (5.23a,b) for Rν and Rκ into (5.4d) determines, via (5.13c), p r 2 p (y2t − y2 )(y−2 1−P t −y ) 2 x= 1−P , where yt = , (5.24a,b) 2y 1+P

sgn y = sgn x and 0 < |y| 6 yt 6 1, since 0 6 P < 1. Here yt is the value of |y| at x = 0, namely the end of the fast mode branch where it merges with the slow mode ω0 = 0. The alternative limit |ω0 | ↑ ∞ corresponds to |y| ↓ 0. The values of F and G (see (5.12c,d)), with Rν , Rκ , Σ and ΣP given by (5.23a–d), are 2

2 F (1 − y2 ) (y−2 t −y ) = , 1−P 4y4

2 G (1 + y2 )(2y−2 t −1−y ) = . (5.25a,b) 1−P 4y2

In turn they determine R0 and m given by (5.12a,b) as functions of y, whose logarithmic y-derivatives give 3R00 1 5 4 2 5 =− 2 + + + − −2 , 2 2yR0 y 1 + y2 1 − y2 y−2 − y − 1 − y2 2y t t 1 3m0 1 1 2 1 − −2 (> 0), = 2+ + + −2 2 2 2 ym y 1+y 1−y yt − y 2yt − 1 − y2

(5.26a) (5.26b)

where here and henceforth the prime is used to denote the y-derivative. Furthermore the product of x2 y and the logarithmic y-derivative of x (defined by (5.24a)) gives xyx0 = −Rν Rκ

(< 0).

(5.26c)

By (5.26b) m and y increase in concert (and by (5.26c) |x| decreases) monotonically from m = 0 as |y| ↓ 0 (|x| ↑ ∞) to m = mt (> 0) at |y| = yt (|x| = 0). The minimization of R0 over m determines the lowest-order fast mode critical values R0 = Rfc0 and m = mfc0 . To that end, we write |y| = |yfc0 | for the value of |y| at which the right-hand side of (5.26a) vanishes. We refrain from adding the subscript ± to our lowest-order critical values, as only one sign actually corresponds to the minimum. In fact, we establish in § 5.2 that the minimum actually corresponds to both ωfc0 > 0 and yfc0 > 0; the sign chosen in all our figures. With this proviso in mind, (5.26a,b) together imply ( m dR0 R00 /R0 < 0 when 0 < m < mfc0 , = 0 (5.27) R0 dm m /m > 0 when mfc0 < m < mt .

402

A. P. Bassom, A. M. Soward and S. V. Starchenko (a) 10

(b) 3 2

5 1 P 0

0 P –1

–5 –2

P –10 –2

–1

0

1

–3

–2

0

–1

1

m

m

F IGURE 8. Plots versus m ˜ = ε−2/3 (m − mt ) of the imperfect frequency pitchfork, as determined by (5.29a,b), for the cases P = Pfct ≈ 0.8 (short-dashed), P = Pc ≈ 0.65 (dotted), P = 0.1 (dots-dashed) (same line style as in figure 6a,b): (a) R˜ = ε−2/3 (R − Rt ) and (b) mω˜ ≈ ε−1/3 mt ω.

Of course, the second range mfc0 < m < mt only exists when |yfc0 | < yt . From (5.26a.b) we deduce that   m dR 0 2(−2 + 17y2t + 8y4t − 3y6t ) = R0 dm t 2 + 7y2t + 4y4t ( Υ > 0 when 0 < P < Pfct , = (5.28a) Γ < 0 when Pfct < P < 1, where we have noted that Γ > 0 (see (5.20d)), and Pfct = 0.79857303 · · ·

(5.28b)

is the solution of 3(3 + P)Υ (P) ≡ 5 + 3P − 10P2 − 2P3 = 0 (see (5.20b)) on 0 < P < 1. Of course, the result (5.28a) follows directly from (5.20a,c). We denote the value of P at which the fast mode minimum Rayleigh number Rfc0 equals the slow mode minimum Rsc0 by Pc ≡ 0.65176592 · · ·. The plot of R0 versus m for P = Pc as well as other values of P is illustrated on figure 6(a). The fast mode is preferred when Rfc0 < Rsc0 and this occurs on 0 < P < Pc . The fast mode still has a minimum Rfc0 on Pc < P < Pfct , but, since Rsc0 < Rfc0 , the slow mode is favoured. Note that fast mode minimum ceases to exist for Pfct 6 P after the slow and fast mode minima have blended at P = Pfct . Evidently the transitional value Pc lies in the range 0.333 · · · ≡ Psct < Pc < Pfct . We also plot in figure 6(a,b), respectively, the curves Rfc0 and mfc0 ωfc0 versus mfc0 as P varies between 0 and Pfct , where mfc0 ωfc0 = 0. The asymptotic behaviours of Rfc0 and mfc0 ωfc0 for P  1 are given by (6.12a) and (6.13a) below. Critical data for the particular case P = 0.2, not illustrated on figure 6, is listed in the right-hand column of table 2. It provides the ε ↓ 0 limiting values of mc , Rc and mc ωc for the finite ε eigenfunctions plots in figure 5(a–d). 5.1.3. The imperfect frequency pitchfork We now consider the nature of the frequency pitchfork at m = mt suggested by the numerical results portrayed in figure 4. For small m − mt the natural blend of the slow

Localized rotating convection and fast mode formula (4.16a) and (5.21a) for the frequencies εωs1 and ω0 is   Γ 2 4εP ω + (m − mt ) ω ≈ , 3 3(1 − P2 )(3 + P) mt

403

(5.29a)

while the fast mode formula (5.22a) for the Rayleigh number continues to hold and on use of (4.16b) may be expressed in the alternative form m − mt R − Rt ω2 + Υs 4 ≈ (1 − 3P) . Rt mt mt

(5.29b)

In Appendix E within the supplementary material, we develop asymptotic expansions valid when m − mt = O(ε 2/3 ), ω = O(ε 1/3 ) and R − Rt = O(ε 2/3 ). Under that scaling the terms in (5.29a) and (5.29b) are O(ε) and O(ε 2/3 ), respectively. The relation (5.29a) is rendered plausible on writing m ˜ = ε −2/3 (m − mt ), ω˜ = ε −2/3 ω and by considering its large |m| ˜ behaviour. Then √ we recover the slow mode ω˜ → 0 as m ˜ → ±∞ and the two fast modes ω˜ → ± (m3t m/Γ ˜ ) as m ˜ ↓ −∞. In the same spirit, we write R˜ = ε−2/3 (R − Rt ) and the slow and fast Rayleigh number limits are recovered from the large |m| ˜ asymptotics of (5.29b). Nevertheless, although the former assertion (5.29a) is plausible, the latter (5.29b) is not obvious because the fast mode theory makes no distinction between Γ ω02 /m4t and (mt − m)/mt . In the supplementary material, we argue why the natural expansion of the Rayleigh number R leading to (5.29b) is in terms of ω2 about the slow mode value Rso as in (5.22a). The nature of the frequency pitchfork is determined by (5.29a) which shows that the slow mode with ω = O(ε) has ω > 0 for m > mt . On decreasing m it makes a smooth transition to a fast mode with ω = O(1), still with ω > 0, for m < mt . In contrast, the slow mode has ω < 0 for m < mt . On increasing m it doubles back on itself in the vicinity of m = mt becoming a fast mode, still with ω < 0. The Rayleigh number R determined by (5.29b) is slave to the behaviour of ω. We plot the tilde variables R˜ , ω˜ versus m, ˜ as determined by (5.29a,b), on figure 8(a,b). The continuous ˜ transition of R in figure 8(a) identified by the ω˜ > 0 branches in figure 8(b) for the cases P ≈ 0.8, P ≈ 0.65 and P = 0.1 compare well with the numerical results for R at P = 0.8, P = 0.65, P = 0.2 portrayed in figure 4(a,c,e) respectively. We should emphasize that, for P  1, (4.16a) shows that the explicit form (5.29a) is only valid when |m − mt |  P and this condition requires ε  P1/2 ( 1). The form of (5.29a) corrected by (4.16a), when ε = O(P1/2 ), is     4P m − mt 13 2 √ ω + (m − mt ) ω ≈ ε − , (5.30) 9 mt 36 3 valid for m − mt = O(P), ω = O(P1/2 ). 5.2. First order A feature illustrated by figure 8(a,b) is that, on decreasing m after the frequency pitchfork at mt , the frequencies of the two fast modes are almost equal in magnitude but possess opposite signs. More significantly, their corresponding Rayleigh numbers are also almost the same but do differ slightly. So in this section we study the first-order corrections R1 and ω1 , which solve (3.11) or, more compactly on use of (3.8),   2β − 3 P(2β − 1) imR1 2β − 3 + 2 + (2β − 1) = 2 . (5.31) −mω1 2 m − iω0 m − iPω0 R0 m − iω0

404

A. P. Bassom, A. M. Soward and S. V. Starchenko 0.5

0

×

–0.5

–1.0

–1.5 0

0.2

m

0.4

0.6

F IGURE 9. Plots of 102 (R0 − R0‡ ) versus m. The P = 0 case is the straight line m-axis (continuous), while other small P (6= 0) cases are P = Pe ≈ 0.016 (chain), P = 0.01 (longdashed) and P = 0.001 (dotted). The curve Rfc0 − R0‡ (> 0) (as P varies) is shown continuous (upper left).

The objective then is to identify which of the two fast modes actually corresponds to the lower Rayleigh number at various values of m and, in particular, at the critical wavenumber mfc . Recall from (5.3a,b) that the sign change ω0 7→ −ω0 leads to θν 7→ −θν , θκ 7→ −θκ , Rν 7→ Rν , Rκ 7→ Rκ and β 7→ β ∗ (see (5.11)), where the superscript ∗ denotes complex conjugate. So, if we denote the values of R1 and ω1 corresponding to ω0± (see (5.13f )) by R1± and ω1± , we see upon taking the complex conjugate of (5.31) that

R1+ = −R1− ,

ω1+ = ω1− .

(5.32)

The sign of ω0 , which gives negative R1 , identifies the mode with the smaller Rayleigh number. The transition from one form of the preferred mode to the other occurs when

R1+ = −R1− = 0.

(5.33)

Interestingly, the neutral (m, R0 ) curves at given P, are nested in the usual way for moderate P, but begin to overlap each other for small P. The overlapping only occurs for very small P(< 0.0163499 . . .) and is not clearly visible in figure 6(a). For that reason, values of R0 − R0‡ , where R0‡ ≡ R0 |P=0 (recall that the superscript ‡ denotes P = 0 values; see below (4.6b,c)) are plotted with considerable magnification on figure 9. From that figure, we see that an envelope of the (m, R0 )-curves exists between A‡el : (m, R0‡ ) = (0, 0) and A‡er : (m, R0‡ ) = (0.57854686 . . . , 0.67221071 . . .) with |ω0‡ | = 1.29635133 . . ., i.e. the two points at which the curves P = 0 and P ↓ 0 intersect. As P is increased from zero, the left-hand point Ael : (mel , Rel ), at which the (m, R )-curve touches the envelope, moves to the right along the envelope. Likewise, the right-hand point Aer : (mer , Rer ) moves to the left. The points coalesce at Ae : (me , Re ), when P = Pe ≈ 0.0163499, me ≈ 0.424203, Re ≈ 0.371354 with |ω0 | ≈ 1.608375.

405

Localized rotating convection

In Appendix A we show that the condition for the points Ael and Aer to lie on the envelope is the same as (5.33), namely R1 = 0. In fact, the normal situation is that R1+ < 0 < R1− . The only exception, for which R1− < 0 < R1+ , occurs when P < 0.0163499 . . . and then only between the points Ael and Aer . In that case mfc0 lies in the interval 0 < mfc0 < mel (see figure 9) outside the interval mel < m < mer between Ael and Aer and so R1+ |fc < 0 < R1− |fc . Consequently for all P the minimizing Rayleigh number Rfc always corresponds to values of mfc , for which the normal situation applies. That critical mode is characterized by a positive frequency ωfc0 > 0:

Rfc = Rfc0 + εR1+ |fc ,

i.e.

Rfc1 = R1+ |fc .

(5.34)

In the low-Prandtl-number limit of the following § 6, this result is confirmed analytically by the results (6.19) and (6.20). 6. Fast modes in the low-Prandtl-number limit P  1 The nature of the fast modes in the low-Prandtl-number limit,

P  1,

(6.1)

ω0 /m2 = O(1),

(6.2)

warrants special attention. To begin, when

we write ‡ ω0± ≈ ω0± ≡ ±m2

p Y 2 − 1,

that is,

x = ±P

p

Y 2 − 1,

(6.3a,b)

see (5.13c,f ), where, of course, Y > 1. Then from (5.3a,b) under the assumptions (6.1) and (6.2), we obtain p (Rν /P)eiθν = 1 ± i Y 2 − 1, Rκ eiθκ = 1 + O(P). (6.3c,d)

These determine Rν = PY and Rκ ≈ 1, with the consequence from (5.9b,c) that Σ ≈ 1 and ΣP ≈ P(Y + 1). In turn (5.12a–d) yields  1/2  1/6 Y +2 R0 R0‡ (Y + 2)3 ≈ ≡ where m ≈ ; (6.4a,b) m m 2Y 2 (Y + 1) 2Y 2 (Y + 1) a solution that breaks away from the steady solution branch at Y = Yt = 1. As Y increases so both m and R0‡ decrease approaching zero as Y ↑ ∞, as illustrated on ‡ figure 6(a). The corresponding behaviour of mω0+ defined by (6.3a,b) is shown on figure 6(b). Finally, the substitution of (6.3c,d) into (5.14) gives " # r Y +2 Y −1 ‡ 1∓i (6.4c) β± ≈ β± ≡ 2Y Y +1

with the real and imaginary parts of β+ plotted on figure 7(a,b). However, the solution (6.3a,b), (6.4a–c) only remains valid when 1 6 Y  P−1 , for which 0 6 |ω0 |/m2  P−1 . Then use of (6.4b) enforces the limitation m  P1/3 . 6.1. Zeroth theory for Pω0 /m2 = O(1)

Since R0 ≈ R0‡ decreases with m throughout the range m  P1/3 we need to consider m = O(P1/3 ) in order to find the minimizing value of the Rayleigh number, which

406

A. P. Bassom, A. M. Soward and S. V. Starchenko

corresponds to ω0 /m2 = O(P−1 ).

(6.5)

For such values the formulation in terms of y rather than x is now particularly helpful with (5.13c) and (5.24a,b) reducing to Pω0 1 − y2 ≡ x ≈ m2 2y

and yt ≈ 1.

(6.6a,b)

This has the useful consequence that (5.3a,b) reduce to Rν e

iθν

iθκ

≈ ix,

Rκ e

(1 − iy)2 = 1 + ix ≈ i 2y

(6.7a,b)

giving (but see also (5.23a,b)) Rν ≈

1 − y2 , 2|y|

eiθν ≈ i(sgn y),

Rκ ≈

1 + y2 , 2|y|

eiθκ ≈ i(sgn y)

(6.7c,d)

1 − iy , 1 + iy

(6.7e,f )

where, from (6.6a), sgn y = sgn x = sgn ω0 . In addition (5.25a,b) with yt ≈ 1 (see (6.6b)) determine 3

F≈

(1 − y2 ) , 4y4

G≈

1 − y4 . 4y2

(6.8a,b)

The expressions (5.12a,b) for the zeroth-order value of the Rayleigh number and wavenumber simplify with (6.8a,b) to 5/3

R0 (1 + y2 ) ≈ P1/3 4y2/3 (1 − y2 )1/3

1/6

and

which, on elimination of P, yields

m |y|1/3 (1 + y2 ) ≈ P1/3 (1 − y2 )1/3

R0 (1 + y2 ) ≈ m 4|y|

,

(6.9a,b)

3/2

.

(6.9c) 2

On multiplication of m6 /P2 , determined by (6.9b), first with (Pω0 /m2 ) as defined by (6.6a), and second with P2 β/m6 , given by (5.11) (but also using (6.7c,e,f )), we obtain the compact expressions (1 + y2 ) mω0 ≈ (sgn y) 2

1/2

and β ≈

1 − iy , 2

(6.10a,b)

respectively. On setting R00 = 0 in (5.26a) with yt = 1, the critical values yfc0 are readily seen to satisfy 3y4fc0 − 6y2fc0 + 1 = 0. The required solution on the permitted range |yfc0 | 6 1 immediately follows as p √ |yfc0 |2 ≈ 1 − 2/3 = 1/(3 + 6). (6.11)

407

Localized rotating convection The substitution of (6.11) into (6.9a–c) determines √ √ 1/3 √ 5/3 1/6 Rfc0 ( 6 − 1) (3 + 6) mfc0 31/12 ( 6 − 1) ≈ , ≈ , √ P1/3 24/3 32/3 P1/3 21/12 (3 + 6)1/6 √ √ 1/2 3/2 Rfc0 ( 6 − 1) (3 + 6) ≈ , mfc0 25/4 33/4 while the values of (6.10a,b) at y = ±|yfc0 | are √ 1/2 ( 6 − 1) 1 mω0± ≈ ± , β± ≈ 3/4 1/4 2 3 2

1∓

(3 +

i √

(6.12a,b) (6.12c)

! 1/2

6)

.

(6.13a,b)

The substitution of (6.12b) into (6.13a) fixes the frequency as √ √ 1/6 1/3 ( 6 − 1) (3 + 6) 1/3 P ω0± ||y|=|yfc0 | ≈ ± . (6.14) 22/3 31/3 Significantly, the critical Rayleigh number Rfc0 given by (6.12a) is an order of magnitude O(P1/3 ) smaller than the critical value Rsc0 given by (4.5c) for the almost steady modes. So for small P the fast modes are certainly preferred to the slow ones. 6.2. First-order theory The result (6.14) identifies two frequencies ω0± = ±|ωfc0 | but only one actually corresponds to the critical mode. To determine which sign fixes the minimizing Rayleigh number, we must consider the next order correction ε R1 given by (5.31). Specifically, aided by (6.6a) and (6.7a,b), the formulae (5.3a,b) reduce to

m2 − iω0 = −iP−1 m2

1 − y2 , 2y

m2 − iPω0 = −im2

(1 + iy)2 , 2y

(6.15a,b)

while (6.10b) determines 2β = 1 − iy. The substitution of these results into (5.31) yields   iy(1 − y2 ) im3 R1 mω1 2 + iy + + (1 − y2 ) = 2 + iy (6.16) 2 2P R (1 + iy) 0 and the vanishing of the real and imaginary parts gives 2

(1 + y2 ) mω1 = − , 1 + 3y2

m3 R1 2y(1 − 2y2 ) =− . P R0 1 + 3y2

(6.17a,b)

Upon further use of the expressions (6.10a) for mω0 and (6.9b) for m/P1/3 , (6.17a,b) may be written as 3/2

ω1± 2 (1 + y2 ) =− |ω0 | 1 + 3y2

,

R1± 2(1 − y2 )(1 − 2y2 ) =∓ R0 (1 + y2 )1/2 (1 + 3y2 )

(6.18a,b)

respectively, where the notation R1± and ω1± identifies the values of R1 and ω1 corresponding to sgn y = ±1 as in (5.32). √ Evidently R1 = 0 at |y| = |yel | = 1/ 2, which with (6.6a) determines |x| = |xel | = 2−3/2 . The corresponding value of m = mel is then fixed by (6.9b): mel ≈ 31/6 P1/3 ,

(6.19)

408

A. P. Bassom, A. M. Soward and S. V. Starchenko

which, upon comparison with mfc0 given by (6.12b), clearly satisfies mfc0 < mel . Furthermore (6.9b) and (6.18b) together show that R1+ ≶ 0 ≶ R1− when m ≶ mel . Consistently at critical, where |y| = |yfc0 | (see (6.11)), (6.18b) yields  3/4 2 R1± 1 =∓ . (6.20) √ 1/2 R0 |y|=|yfc0 | 3 ( 6 − 1) From this result we can conclude that Rfc1 = R1+ |fc < 0 (see (5.34)) linked to ω0 > 0, and so the minimizing Rayleigh number corresponds to ωfc0 > 0 exactly as predicted by the finite P results of § 5.2. A comprehensive first-order theory for the critical values (5.1a–c) requires us to consider the corresponding expansion of yfc . This takes the form yfc = yfc0 + εyfc1 + · · · ,

(6.21a)

0 = ε−1 R 0 (yfc ) = R000 (yfc0 )yfc1 + R10 (yfc0 ).

(6.21b)

where, correct to lowest order, yfc1 solves

The solution yfc1 = −R10 (yfc0 )/R000 (yfc0 ) fixes the first-order critical value corrections

Rfc1 = R1 (yfc0 ),

ωfc1 = ω00 (yfc0 )yfc1 + ω1 (yfc0 ),

mfc1 = m0 (yfc0 )yfc1 ,

(6.22a,b) (6.22c)

which appear in the expansions (5.1a–c). Since evaluation of these quantities is tedious and not particularly illuminating, we do not pursue the matter any further here. 7. Discussion The analysis of rotating Boussinesq convection undertaken in this paper invokes the two limits of small E and small ε taken sequentially. Thorough investigations of the onset of instability in the former rapid rotation limit (E  1), when the unstable stratification is caused by a temperature gradient T 0 (r∗ ) ≡ −r∗ Q ∗ (r∗ ), have been undertaken by both fully numerical and semi-analytical methods for the cases Q ∗ ∝ 1 of a uniform distribution of heat sources (Jones et al. 2000) and Q ∗ ∝ (ri∗ /r∗ )3 of differential heating maintained by a temperature drop between the inner and outer boundary (Dormy et al. 2004). In the latter case the temperature gradient is relatively steep near the inner core boundary (r∗ = ri∗ ) with the consequence that the ‘cartridge belt’ of convection rolls resides adjacent to (but outside) the inner core tangent cylinder. There an analytic reduction of the problem is possible to ordinary differential equations with the axial coordinate z∗ as the independent variable, which may be solved numerically. In our corresponding small-E semi-analytic study, we have blended the models and reversed the sign of the heat sources so that they become ‘heat sinks’. This artificial device is adopted to provide a simple model of a smoothly varying stratification and should not be taken as a literal interpretation of the underlying physics. The relative sizes of differential heating and the heat sinks has been adjusted to yield an unstable layer adjacent to the inner core out to a radius rn∗ . In the small-ε thin layer limit δn∗ = rn∗ − ri∗ = O(ε2 ri∗ ), our asymptotics depends only on the leading-order approximation Q ∗ ≈ Qi∗ (rn∗ − r∗ )/δn∗ (see (2.15a)), which is insensitive to the precise form of the stratification. The onset of convection continues to be confined to the tangent cylinder but its axial extent away from the equatorial plane is restricted to be O(εri∗ ) (see figure 1). Our fully analytic theory takes

Localized rotating convection

409

advantage of that short length and we have argued that the correct ordering of the various length scales is (2.10a) for the leading-order asymptotic theory to be valid. At finite-ε, our numerical solution in § 2.3 of the ordinary differential equations (2.21a,b) subject to the boundary conditions (2.22a,b) reveals trends towards the small-ε regime, for which the later analytic asymptotic studies apply. In the case of the critical modes at M = Mc (say), the plots of Re{uc } in figures 3(a) and 5(a), show that for ε = 1, the classical ‘cartridge belt’ picture of Taylor columns independent of z∗ is roughly achieved. (Note that, although Im{uc } in figures 3(c) and 5(c) is obviously z∗ -dependent, its magnitude relative to Re{uc } is small confirming that the entire uc is almost z∗ -independent.) However, once ε is as small as 10−0.5 the solution is evidently z∗ -dependent decaying to zero before the outer boundary is reached. With further decrease in ε the eigenfunction continues to change independent of the precise form of the outer boundary condition. It is only at smallish ε below roughly 10−1 that an asymptotic structure emerges and only then do the R (m) and ω(m) plots in figures 2 and 4 suggest any distinction between slow and fast modes. Indeed the asymptotic development of § 5 reveals that, in the limit ε ↓ 0, fast modes can only occur for a restricted range of azimuthal wavenumbers M(< Mt , say, dependent on the Prandtl number P). Since Mt → 0 as P ↑ 1, the fast modes only come into play for P < 1 and this sharp bound is a facet of the small-ε limit. At finite, but small, ε the fast and slow modes blend one into the other in the vicinity of M = Mt . Such a smooth transition is a feature common with the classical sphere and shell models. Owing to the z∗ -variation of Q on the tangent cylinder, our solutions at small ε are manifestly z∗ -dependent. The consequences of this z∗ -dependence are however fundamentally different at the zeroth and first order of the asymptotics. At zeroth order, motion is driven by the axial torque −iMRθ (see (2.14a)), while the temperature perturbations θ are achieved by the radial-s∗ convection Q sU (see (2.14c), (2.12b)) of heat. This, at any rate, corresponds to the plane layer analogues, particularly Chandrasekhar (1961), described in Appendix B.1. Nevertheless, at first order, the generation of the axial pressure gradient dp/dz = −E−1 dΨ/dz (see (2.14d)) driven by the axial buoyancy force zRθ (see (2.14b)) and the axial convection Q zW (see (2.14c))

(7.1a) (7.1b)

of heat, both ignored completely throughout Appendix B, together contribute to subtle corrections of the lowest-order solution. The two effects (7.1a,b), crucial to our firstorder development, may be traced to the factor iεmC appearing in both (2.21a) and (2.21b). Indeed, both our numerical results for ε as large as 10−0.5 and our asymptotic results for ε  1 indicate that the flow amplitude at the outer boundary is negligible with the consequence that our solutions are insensitive to the nature of the outer boundary conditions. The revelation, that the onset of instability in our model is consistently manifest by prograde travelling waves despite the lack of dependence on the outer boundary topography, is significant. For it is often stated that the prograde movement of the cartridge belt convective cells in rotating shells relies on topography and so is best described by the term ‘thermal Rossby wave’, following Busse’s (Busse 1970a) pioneering annulus model studies (see Appendix B.2). We should add that our results do not depend directly on the cylindrical geometry either. It is only when the radial-s∗ structure is considered, as hinted at by (2.18a,b), that cylindrical effects (see particularly the s∗ -dependence of A ≈ M/s, (2.12d)) influence the solution. That consideration is outside the scope of our present study, but its influence on limiting the radial-s∗ extent and on the azimuthal tilting of the convecting columns in the cartridge

410

A. P. Bassom, A. M. Soward and S. V. Starchenko

belt is illustrated by the equatorial cross-section portrayed in figure 4 of Dormy et al. (2004). An intriguing feature of our small-ε results is that, at lowest order in our expansions in powers of ε, the axial structure is controlled by the complex exponential

E (φ, z, t) = exp[i(Mφ − Im{β/2} (z/ε)2 −Ωt) − Re{β/2} (z/ε)2 ]

(7.2)

(see (2.7), (2.16) and (3.2a,b)), where the dimensionless variables, z, t and Ω determine the axial distance z∗ = ri∗ z, time t∗ = (ri∗ 2 /ν)t and frequency Ω ∗ = (ν/ri∗ 2 )Ω, respectively. Whereas the radial-s∗ velocity (U) is proportional to E (see (3.2a)) the axial velocity (W) is proportional to zE (see (3.2b)) and so its modulus |W| achieves its maximum at p √ √ ε(Λ/ 2)ri∗ with Λ/ 2 = 1/ Re{β}, (7.3a,b)

where Λ is the ratio of the e-folding height of U to the axial height z∗n = εri∗ of the unstable region. In the case of the critical slow modes Ω = Im{β} = 0, this e-folding ratio for the convecting columns, determined by the result (4.5a), is p Λsc = 8/5 = 1.26 . . . (7.4a)

with the implication that there is moderate penetration into the stable region above. For the critical fast modes, the penetration height ratio Λfc at critical increases and reaches its maximum Λfc ↑ 2

as P ↓ 0

(7.4b)

(see figure 7a and (6.13b)). Another feature of the fast modes revealed by (7.2) is that the axes of the propagating columns cease to be aligned in the z∗ -direction but bend to reveal an azimuthal component. With M > 0, figures 6(b) and 7(b) demonstrate that Ω and Im{β} possess opposite signs. Furthermore all critical modes are characterized by Ωfc > 0 for which Im{βfc } < 0. Since Ωfc /Mfc > 0 the critical modes have prograde propagation, while the column boundaries Mfc φ − Im{βfc /2} (z/ε)2 = constant bend in the retrograde (−φ)-direction with increasing |z| away from the equatorial plane. At the e-folding height the phase angle shift achieved is |Im{β}/Re{β}|. Again figure 7(a,b) and (6.13b) indicate that the maximum phase angle shift reaches its maximum √ −1/2 = 0.428 · · · as P ↓ 0. (7.5a,b) |Im{βfc }/Re{βfc }| = (3 + 6)

This corresponds to a relatively small phase angle of roughly 25◦ . So the axial bending predicted by the asymptotic theory is not particularly pronounced and largely hidden by the decaying convection amplitude with increasing |z|, even when P  1. This bending is also likely to occur in the classical shell models of Jones et al. (2000) and Dormy et al. (2004) but as they did not comment upon the matter, it is presumably not particularly noticeable in those models either. Perhaps the most intriguing feature of our results is the critical Rayleigh number power law behaviour Rfc = O(P1/3 ) (see (6.12a)), as P ↓ 0. As we argue in Appendix B.3, if we were to replace the ‘soft’ interface r∗ = rn∗ between the stable and unstable regions by a rigid or ‘hard’ boundary, the critical Rayleigh number would be O(P) smaller: Rfc = O(P4/3 ). This is a potentially disturbing result, which prompts us to question whether we really have identified the mode that minimizes the Rayleigh number. Reassuringly our numerical integrations of the z∗ -dependent second-order

Localized rotating convection

411

ordinary differential equation tend support our findings and the physical explanation of our result is relatively simple. The traditional picture of convection between rigid boundaries is that of gently excited inertial waves. They are selected by their ability to advect heat efficiently and this is sufficient to produce a weak buoyancy force that maintains the motion against the small viscous dissipation. In our system the situation is more complicated, for any inertial wave existing in the unstable region propagates energy out of it into the effectively infinite stably stratified region above, where the energy is viscously dissipated. There is essentially no wave reflection at the ‘soft’ boundary interface between stable and unstable layers. In consequence the buoyancy force needs to be involved in the zeroth-order force balance in the thin unstable region (note the R0 in (3.7c)) in order to maintain the convection. This important distinction between the ‘hard’ and ‘soft’ boundary configuration accounts for the O(P) disparity mentioned. Acknowledgements We are particularly grateful to Dr E. Dormy for his continued interest and constructive input. We thank Professor F. H. Busse, Professor C. A. Jones and Professor K. Zhang for useful discussions and anonymous referees for helpful comments. S. V. S. is grateful for the support of RFBR grant No. 09-05-00979.

Supplementary material is available online at journals.cambridge.org/flm. Appendix A. The envelope In this Appendix we demonstrate that the condition for the points Ael : (mel , Rel ) and Aer : (mer , Rer ) on the stability boundary for given P to lie on the envelope of those curves generated by varying P is the same as (5.33), namely R1 = 0. The values of the Rayleigh number R0 and the frequency ω0 for the fast mode at given wavenumber m are determined by their real value solution of (3.8). As m varies so a curve is generated in the (m, R0 )-plane. Within a restricted parameter range, under variation of the Prandtl number P, an envelope is generated and to find it we consider the relation between the differentials dR0 , dω0 and dP at fixed m. Differentiation of (3.8) and its further use yields

idω0 , m2 − iω0 dR 0 idω0 i(Pdω0 + ω0 dP) 2β −1 dβ = − 2 + . R0 m − iω0 m2 − iPω0

(2β − 1)dβ = −2im2 (m2 − iω0 ) dω0 = −2β(β − 1)

(A 1a) (A 1b)

Elimination of dβ determines

dR 0 dω0 Pdω0 − (2β − 3) 2 − (2β − 1) 2 R0 m − iω0 m − iPω0 ω0 dP . = (2β − 1) 2 m − iPω0

i(2β − 1)

(A 1c)

The envelope is determined by requiring dR0 = 0 for some non-zero real dω0 and dP. This determines the condition   2β − 3 m2 − iPω0 Im = 0. (A 2) 2β − 1 m2 − iω0

412

A. P. Bassom, A. M. Soward and S. V. Starchenko

Interestingly, the first-order corrections ε R1 and εω1 fixed by (5.31) satisfy i(2β − 1)

mR1 mω1 Pmω1 1 − (2β − 3) 2 − (2β − 1) 2 = (2β − 3) 2 . R0 m − iω0 m − iPω0 m − iω0 (A 3)

Evidently, the condition for R1 = 0 for some real ω1 is again (A 2). This remarkable result confirms that R1 = 0 is achieved on the envelope of the (m, R0 ) curves under variation of P. Appendix B. A plane layer model Here we consider the relationship between our work and some of the classical models for convection. First consider a plane layer model of our system, in which the s-direction is vertical and the temperature gradient is constant. This configuration continues to be described by (2.14) but with the terms zREθ and zW, important to our first-order theory (see (7.1a,b)), omitted from (2.14b,c) as in our zeroth-order development. Accordingly (2.21a,b) become   dw RQ m2 dψ = (a2 − iω)w, = a2 (a2 − iω) − 2 ψ (B 1a,b) dζ dζ a − iPω

together with the simplification

Q = 1.

(B 1c)

Since the coefficients in the second-order system (B 1a,b) are all constants, the solution with the symmetry dψ/dz = w = 0 at the equator z = 0 has the form ψ = ψˆ cos(nζ ),

w=w ˆ sin(nζ )

(B 2a,b)

for some complex constant n. We make the long radial (now vertical) length scale approximation a ≈ m. Then the equations (B 1a,b) are satisfied by (B 2a,b) when

m2 − iω (B 3a,b) , where β = n2 . m2 − iPω The formula (B 3a) is simply our earlier fundamental zeroth-order relation (3.8)1 . Whereas the non-uniform stratification studied in our paper leads to (3.8)2 , which may be regarded as a supplementary formula for β, here we consider two distinct boundary conditions at ζ = z/ε = H (say), that provide alternative formulae for the value of n in (B 3b). The discussion of these cases in §§ B.1 and B.2 helps to place our solutions of §§ 4 and 5 in context, as we expound in § 7. 2

β + m2 (m2 − iω) = m2 R

B.1. A rectangular annular duct The simplest lateral closure of our plane layer model occurs when a plane boundary is located at the latitude z = εH, at which w = 0, so generating a rectangular annular duct. This boundary condition is satisfied when

sin(nH) = 0 2

⇒

nH = Nπ

(N = 0, 1, 2, · · ·),

(B 4a,b)

where, significantly, β = n is real. For the steady modes ω = 0, the relation (B 3a,b) determines

Rs ≈ (n/m)2 +m4 ,

(B 5)

413

Localized rotating convection from which we obtain the critical values

Rsc = 3m4sc ,

1/3

m2sc = (n2 /2)

,

(B 6a,b)

analogous to (4.5c,d). For fixed n 6= 0, the results (B 5) and (B 6a,b) compare qualitatively with those for the zeroth-order slow modes theory of § 4.1. More subtle effects, such as the small non-zero frequency predicted by the first-order shell theory of § 4.2, are not captured by our duct model which neglect the crucial terms (7.1a,b). It should not be overlooked, however, that the minimizing Rayleigh number for the rectangular duct is determined by the n = 0 mode corresponding to geostrophic motion. Then, according to (B 5), the n = 0 critical values are reached under the limit R ↓ 0 as m → 0. Despite that cautionary comment, physically realistic models with z-dependent stratification, manifest by non-constant Q in (B 1b), do not admit a geostrophic n = 0 mode solution. Although the shell fast modes of § 5 clearly correspond to complex β, it is worth exploring the consequences of the duct model result (B 3a) with real β = n2 for the overstable modes ω 6= 0. The imaginary and real parts of (B 3a) determine " #  P n 2 1 − P  n 2 4 R = 2(1 + P) + m , ω2 = − m4 . (B 7a,b) 1+Pm 1+P m Obviously, like the zeroth-order fast mode theory of § 5.1, the modes (B 7) only exist for P < 1 as in (5.8a). Furthermore, since the stability boundary (B 7a) only exists for   1 − P 2 1/3 m2 < m2t ≡ n (P < 1), (B 7c) 1+P

its features are qualitatively similar to those of the shell fast modes portrayed in figure 6(a) (cf. also ω in figure 6(b)). Furthermore the critical overstable duct mode predicted by (B 7a) is characterized by

Rfc = 6(1 + P)m4fc , m2fc =

which only exists when P <

√ 2/3.



ωfc2 = [(2/P2 ) − 3]m4fc ,

P2 n2 2 (1 + P)2

1/3

,

(B 8a,b) (B 8c)

B.2. The Busse annulus In view of the geostrophic degeneracy noted below (B 6), Busse (1970a) proposed his annulus model obtained by tilting the z = εH boundary by a small angle tan−1 (ηH). To realize the features of his model within our formulation, it is necessary that H itself is small and that η = O(1) is a constant. Then the boundary condition (2.22b) on this tilted surface (z/H) − ε = η(1 − s) becomes imηHψ + w = 0, which on substitution of (B 2a,b) yields

imηH ψˆ cos(nH) + w ˆ sin(nH) = 0.

(B 9a)

(n/H) tan(nH) = iηm(m2 − iω).

(B 9b)

with −nψˆ = (m − iω)w ˆ from (B 1a). Hence (B 9a) becomes 2

When H  1, (B 9b) determines the spectrum n ≈ Nπ/H (N = 0, 1, 2, . . .) as before in (B 4b). Nevertheless, Busse’s point is that the tilted boundary removes the

414

A. P. Bassom, A. M. Soward and S. V. Starchenko

degeneracy of the geostrophic mode, so that the lowest mode nH ≈ 0 is determined more precisely by (β =)n2 ≈ iηm(m2 − iω).

The marginal modes have

R≈



ηP (1 + P)m

2

+ m4 ,

ω≈

η (1 + P)m

(B 10)

(B 11a,b)

and propagate with the character of Rossby waves. They determine the critical values  1/3 η P2 η 2 4 2 Rc = 3mc , ωc = , mc = . (B 12a,b,c) (1 + P)mc 2 (1 + P)2 B.3. Appraisal The rectangular duct model of Appendix B.1 captures features of the zeroth-order theories of the slow (§ 4.1) and fast (§ 5.1) shell modes. It should not be overlooked that the first-order shell theories of §§ 4.2 and 5.2 actually show that the critical slow and fast modes lie on the same stability boundary. That is a characteristic of the Busse annulus result (B 11a), for which the P  1 critical values have the slow mode features mc = O(1), ωc = O(P−1 ) (see (4.5d) and (4.8)) and the P  1 critical values have fast mode features mc = O(P1/3 ), ωc = O(P−1/3 ) (see (6.12b) and (6.14)). On the one hand, the symmetries of the rectangular duct model of Appendix B.1 prevent any sign selection of the frequency. On the other hand, it is encouraging to find that positive value of the phase velocity (ω/m)c at critical, evinced by Busse’s topographic waves of Appendix B.2 for all P, also happens in our shell model due instead to the z-dependencies (7.1a,b). It is also of interest to note that when P  1 the leading-order force balance in both the duct and Busse annulus models corresponds to that for inertial waves, which is manifest by the approximation at critical of (B 3b) to n2 = βc ≈ (mc ωc )2 , a real positive constant. Despite the fact that the scalings of mc from both the rectangular duct and Busse annulus models at small and large P agree with those determined by our shell model, that is not the case for the critical Rayleigh number R . Although the results (B 6a,b) and (B 12a,c) determine Rc = O(1) for P  1 consistent with (4.5c), similar consistency is not achieved at small P. In fact, for P  1 both (B 8a,c) and (B 12a,c) predict Rc = O(P4/3 ) in striking contrast to the shell model estimate Rc = O(P1/3 ) determined by (6.12a). This surprising discrepancy of the small factor P draws attention to the fact that the outer stable layer in the shell model does not behave like a rigid boundary. Rather it acts to absorb and damp any inertial waves generated in the thin unstable layer beneath it. Indeed, since βfc0 defined by (6.13b) is fully complex, it cannot be balanced by the real (mfc0 ωfc0 )2 alone and so is not simply an inertial wave. By necessity therefore the buoyancy force is involved in the leading-order force balance. REFERENCES

A L -S HAMALI, F. M., H EIMPEL, M. H. & AURNOU, J. M. 2004 Varying the spherical shell geometry in rotating thermal convection. Geophys. Astrophys. Fluid Dyn. 98, 153–169. AURNOU, J. M. & O LSON, P. L. 2001 Strong zonal winds from thermal convection in a rotating spherical shell. Geophys. Res. Lett. 28 (13), 2557–2559.

Localized rotating convection

415

B RAGINSKY, S. I. 1993 MAC-oscillations of the hidden ocean of the core. J. Geomagn. Geoelectr. 45, 1517–1538. B RAGINSKY, S. I. & ROBERTS, P. H. 1995 Equations governing convection in Earth’s core and the geodynamo. Geophys. Astrophys. Fluid Dyn. 79, 1–97. B USSE, F. H. 1970a Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441–460. B USSE, F. H. 1970b Differential rotation in stellar convection zones. Astrophys. J. 159, 629–639. B USSE, F. H. 1975 A model of the geodynamo. Geophys. J. R. Astr. Soc. 42, 437–459. B USSE, F. H. 2002 Is low Rayleigh number convection possible in the Earth’s core? Geophys. Res. Lett. 29 (7), 1105. doi:10.1029/2001GL014597 3pp. C HANDRASEKHAR, S. 1961 Hydrodynamic and Hydromagnetic Stability. Oxford, Clarendon Press. C HRISTENSEN, U. R. 2006 A deep dynamo generating Mercury’s magnetic field. Nature 444, 1056–1058. C HRISTENSEN, U. R. & W ICHT, J. 2008 Models of magnetic field generation in partly stable planetary cores: applications to Mercury and Saturn. Icarus 196, 16–34. D ORMY, E., S OWARD, A. M., J ONES, C. A., JAULT, D. & C ARDIN, P. 2004 The onset of thermal convection in rotating spherical shells. J. Fluid Mech. 501, 43–70. F EARN, D. R. & L OPER, D. E. 1981 Compositional convection and stratification of the Earth’s core. Nature 289, 393–394. G REENSPAN, H. P. 1968 The Theory of Rotating Fluids. Cambridge University Press. J ONES, C. A., K UZANYAN, K. M. & M ITCHELL, R. H. 2009 Linear theory of compressible convection in rapidly rotating spherical shells, using the anelastic approximation. J. Fluid Mech. 634, 291–319. J ONES, C. A., S OWARD, A. M. & M USSA, A. I. 2000 The onset of thermal convection in a rapidly rotating sphere. J. Fluid Mech. 405, 157–179. ´ N ET, M., G ARCIA, F. & S ANCHEZ , J. 2008 On the onset of low-Prandtl-number convection in rotating spherical shells: non-slip boundary conditions. J. Fluid Mech. 601, 317–337. ROBERTS, P. H. 1965 On the thermal instability of a highly rotating fluid sphere. Astrophys. J. 141 (1), 240–250. ROBERTS, P. H. 1968 On the thermal instability of a rotating-fluid sphere containing heat sources. Phil. Trans. R. Soc. Lond. A 263, 93–117. Sˇ IMKANIN, J., B RESTENSK Y´ , J. & Sˇ EV Cˇ ´I K, S. 2003 Problem of the rotating magnetoconvection in variously stratified fluid layer revisited. Stud. Geophys. Geod. 47, 827–845. Sˇ IMKANIN, J., B RESTENSK Y´ , J. & Sˇ EV Cˇ ´I K, S. 2006 On hydromagnetic instabilities and the mean electromotive force in a non-uniformly stratified Earth’s core affected by viscosity. Stud. Geophys. Geod. 50, 645–661. ˇ ´ , D. 2009 Convection in rotating non-uniformly Sˇ IMKANIN, J., H EJDA, J. & JANKOVI COV A stratified spherical fluid shells: a systematic parameter study. Contrib. Geophys. Geod. 39 (3), 207–220. S OWARD, A. M. 1977 On the finite amplitude thermal instability of a rapidly rotating fluid sphere. Geophys. Astrophys. Fluid Dyn. 9, 19–74. S TANLEY, S. & B LOXHAM, J. 2004 Convective-region geometry as the cause of Uranus’ and Neptune’s unusual magnetic fields. Nature 428, 151–153. S TANLEY, S. & B LOXHAM, J. 2006 Numerical dynamo models of Uranus’ and Neptune’s magnetic fields. Icarus 184, 556–572. S TANLEY, S. & G LATZMAIER, G. A. 2010 Dynamo models for planets other than Earth. Space Sci. Rev. 152, 617–649. S TANLEY, S., B LOXHAM, J., H UTCHISON, W. & Z UBER, M. 2005 Thin shell dynamo models consistent with Mercury’s weak observed magnetic field. Earth Planet. Sci. Lett. 234, 27–38. S TANLEY, S. & M OHAMMADI, A. 2008 Effects of an outer thin stably stratified layer on planetary dynamos. Phys. Earth Planet. Inter. 168, 179–190. S TARCHENKO, S. V., KOTELNIKOVA, M. S. & M ASLOV, I. V. 2006 Marginal stability of almost adiabatic planetary convection. Geophys. Astrophys. Fluid Dyn. 100, 397–427. TAKEHIRO, S. & H AYASHI, Y.-Y. 1995 Boussinesq convection in rotating spherical shells – a study on the equatorial superrotation. In The Earth’s Central Part: Its Structure and Dynamics (ed. T. Yukutake). pp. 123–156. Terra Scientific Publishing Co.

416

A. P. Bassom, A. M. Soward and S. V. Starchenko

TAKEHIRO, S. & L ISTER, J. R. 2001 Penetration of columnar convection into an outer stably stratified layer in rapidly rotating spherical fluid shells. Earth Planet. Sci. Lett. 187, 357–366. TAKEHIRO, S. & L ISTER, J. R. 2002 Surface zonal flows induced by thermal convection trapped below a stably stratified layer in a rapidly rotating spherical shell. Geophys. Res. Lett. 29 (16), 1803. doi:10.1029/2002GL015450 4pp. TAKEHIRO, S., YAMADA, M. & H AYASHI, Y.-Y. 2011 Retrograde equatorial surface flows generated by thermal convection confined under a stably stratified layer in a rapidly rotating spherical shell. Geophys. Astrophys. Fluid Dyn. 105, 61–81. W ICHT, J. & T ILGNER, A. 2010 Theory and modelling of planetary dynamos. Space Sci. Rev. 152, 501–542. Z HANG, K. & L IAO, X. 2004 A new asymptotic method for the analysis of convection in a rapidly rotating sphere. J. Fluid Mech. 518, 319–346. Z HANG, K., L IAO, X. & B USSE, F. H. 2007 Asymptotic solutions of convection in rapidly rotating non-slip spheres. J. Fluid Mech. 578, 371–380. Z HANG, K. & S CHUBERT, G. 2000 Teleconvection: remotely driven thermal convection in rotating stratified spherical layers. Science 290 (5498), 1944–1947. Z HANG, K. & S CHUBERT, G. 2002 From penetrative convection to teleconvection. Astrophys. J. 572, 461–476.