Unified Linear Stability Analysis for Thermal Convections in Spherical Shells Under Different Boundary Conditions and Heating Modes

Layered interior structure of most planets may be well approximated by spherical shells. An important subject for analyzing planet thermal structure and evolution is to determine the transitional condition from heat conduction dominated mechanism to thermal convection dominated mechanism inside a spherical shell, which can be accomplished through linear stability analysis. In this paper, we propose an analytical solution in terms of Frobenius series to the linear stability analysis that systematically unifies various boundary conditions, heating modes, and shell geometries characterized by inner‐out radius ratio. The analysis is eventually reduced to solving the critical Rayleigh number together with its corresponding spherical harmonic degree from a vanishing determinant of sixth order, and formulas for calculating the elements of the determinant are explicitly provided. The obtained high‐accuracy results clearly show that thermal instability depends strongly on the constraining strength posed by the boundary conditions, heating modes, and shell geometry. The critical Rayleigh number generally decreases with increasing internal heating fraction and inner‐out radius ratio, as well as with switches from fixed to free‐slip or from isothermal to constant heat flux in boundary conditions, whereas remarkable exceptional cases do exist. A thorough analysis is provided to fully explain the general trends and exceptional cases, which uncovers the complicated interplay between shell geometry, heating modes, and boundary conditions, and thus substantially improves our understanding of the underlying physics for thermal convections in spherical shells. Flow streamlines and perturbation temperature distributions for some typical cases of axisymmetric thermal convections are presented as well.


Introduction
As layered interior structure of most planets may be well approximated by spherical shells, determining the transitional condition from heat conduction dominated mechanism to thermal convection dominated mechanism inside a spherical shell forms an important subject for analyzing planet thermal structure and evolution. Therefore, linear stability analysis for thermal convections in spherical shells attracts a lot of research interests and has profound influence on many geodynamic topics, including mantle convection pattern and large-scale mantle structure (e.g., Bercovici, 2015;Davies, 1977;Oxburgh & Turcotte, 1978;Zhong & Liu, 2016).
Historically, Sir Benjamin Thompson (Count von Rumford) is often credited as the first researcher to experimentally observe convections (Brown, 1967), and Henri Bénard's first observations of cellular convection inspired Lord Rayleigh to develop linear stability analysis theory for thermal convection, which was focused on the first phase or the transient onset (Bercovici, 2015). Up till the early 1960s, the onset of convective instability from infinitesimal perturbations was well studied across several fields, as summarized in the classical treatise by Chandrasekhar (1961). Recent studies have shown the increasing power of numerical simulation in investigating thermal convections, but it remains a challenge to understand the underlying physics of the modeling results (Ogawa, 2008;Zhong & Liu, 2016). A classical linear stability analysis is still important for inferring dominant physical mechanism forming convective structure (Turcotte & Schubert, 2014;Zhong & Liu, 2016). Furthermore, results of linear stability analysis usually provide appropriate starting values for iterative nonlinear convection calculations (e.g., Hsui et al., 1972;Palm, 1975;Roberts, 1967;Tveitereid & Palm, 1976;Zebib, 1993).
Previous studies have developed various methods for linear stability analysis of thermal convection in spherical shells. It is common to expand the longitude and latitude components of the involved functions in terms of spherical harmonics, whereas the radial component is represented by various approaches. For example, Feffreys and Bland (1950) used Taylor series expansion, Chandrasekhar (1952) and Backus (1955) used the Bessel function expansion, Avila et al. (2013) adopted the Chebyshev polynomial expansion, while Zebib et al. (1980Zebib et al. ( , 1983 approximated the radial temperature variations by sinusoidal Fourier series and the radial poloidal function by numerical integration. With the improvement of computing power in recent years, propagative matrix method (e.g., Buffett et al., 1994;Liu & Zhong, 2013) and finite difference method (e.g., Kameyama et al., 2013) were also used in linear stability analysis.
In spite of some interesting cases being successfully investigated based on existing methods, a systematic linear stability analysis for thermal convections in spherical shells has never been conducted, and a unified approach that is applicable to various boundary conditions, heating modes, and shell geometry has not been established. Traditional expansion methods, such as those based on Taylor series, Bessel functions, Chebyshev polynomials, and Fourier series, are appropriate only for solving problems subject to particular type boundary conditions combined with specific heating modes, while the numerical methods cannot ensure sufficient computational accuracy, especially for thin spherical shells.
In this paper, we propose a new approach based on Frobenius series (a special case of Laurent series) expansion and the classical Frobenius method (Asmar, 2004). By analyzing the principal terms for the characteristic equation, we prove in Appendix A that there exists a logarithmic branch point (Courant & Hilbert, 1989;Ponnusamy & Silverman, 2006) in solutions to thermal convections in spherical shells. The difficulty raised by the logarithmic branch point cannot be effectively and efficiently dealt with by the existing expansiontype or numerical methods, whereas Frobenius series expansion with an additional logarithmic term is the natural choice for solving such problems. As a result, our new method is applicable for solving problems associated with various shell geometry subject to various boundary conditions and heating modes. Similar to linear stability analyses for convections in porous medium (Beck, 1972;Hong et al., 2008;Kim et al., 2008;Kim & Kim, 2002;Nield, 1968Nield, , 1975Nield & Bejan, 2017;Noghrehabadi et al., 2013), we systematically investigate the effects of boundary condition by considering each of the outer and inner boundaries of a spherical shell as either an isothermal or constant heat flux one as well as either a fixed or free-slip one, resulting in 16 combinations in total. In addition, we consider general combined heating that allows continuous variation from pure basal heating to pure internal heating, and adopt an arbitrary precision computation package (Bailey, 2017;Bailey et al., 2002) to ensure that all the involved computations are sufficient to provide critical Rayleigh numbers accurate to at least six decimal places for all the shell geometries continuously varying from a sphere to a thin shell with an inner-outer radius ratio of 0.95, facilitating a thorough investigation of the shell curvature effect (Zhang & Busse, 1997).
In addition to the novelty and universality of our proposed new method, this study contributes significantly in three aspects. First, our results clearly reveal the general trends and some remarkable exceptional cases for critical Rayleigh number depending upon shell geometries, heating modes, and boundary conditions. A thorough analysis is provided to fully explain the general trends and exceptional cases, which uncovers the complicated interplay between shell geometry, heating modes, and boundary conditions, and thus substantially improves our understanding of the underlying physics within thermal convections in spherical 10.1029/2019EA000672 Earth and Space Science shells. Second, the explicit formulas given in Appendix C and FORTRAN codes given in the supporting information serve as convenient tools and provide access with no technical barrier for researchers in need. Finally, the presented results of high accuracy offer powerful means for benchmarking numerical simulations of thermal convections, which are nowadays widely employed but extremely difficult to quantify computational errors without comparison to analytical solutions.

Methods
Prior to the onset of thermal convection, heat conduction is the only heat transfer mode inside the modeled spherical shell. The conductive temperature in the shell is considered as spherically symmetric and governed by the heat conduction equation where ρ 0 is the fluid density, k is the thermal conductivity, and H is the volumetric heat production rate, all considered as constants. The conductive temperature T c and its radial gradient β c (r)can be obtained by successively integrating equation (1), that is, where c and c' are the integration constants.
It is readily obtained from equation (3) that the temperature rise from the outer boundary (r = R 1 ) to the inner boundary (r = R 2 ) of the spherical shell can be expressed as Thermal convection occurs when this temperature difference is sufficiently large.
In order to establish a unified treatment that accounts for all the heating modes, including basal heating, internal heating, and combined heating, it is necessary to find in the temperature difference given by equation (4) the portion that is solely caused by internal heating. Letting the inner boundary be adiabatic (i.e., an insulating boundary), in accordance with internal heating, equation (2) yields Substituting equation (5) into equation (4) gives the temperature difference caused by internal heating Therefore, the contribution of internal heating in the temperature difference can be measured by the internal heating fraction defined by with η being the radius ratio defined by The temperature difference given by equation (4) naturally leads to a Rayleigh number definition for combined heating and a partition for the contributions from internal heating and basal heating, that is, where ρ 0 is the fluid density, μ is the viscosity, κ is the thermal diffusivity, α v is the thermal expansion coefficient, and g 0 is the acceleration of gravity at the outer boundary, all being considered as constants; subscripts ch, ih, and bh stand for the combined heating, internal heating, and basal heating, respectively, and b is the thickness of the spherical shell: In the literature, Rayleigh number for internal heating is usually defined by It is straightforward to verify This serves as a direct conversion relationship between the heat production rate based Ra H and the temperature difference based Ra ih , while the latter allows for comparing the internal heating effect with the basal heating effect on a common basis.
Once thermal convection occurs in the spherical shell, it is governed by the conservational laws of mass, momentum, and energy. Based on the Boussinesq approximation (e.g., Bejan, 2006Bejan, , 2013Oertel, 2004;Schubert, 1979), temperature-dependent density variation around a reference temperature T ref is only important in the buoyancy effect. As a result, the perturbation fluid velocity u ! , pressure p, and temperature δT = T − T c for the onset of thermal convection satisfy the following governing equations (Schubert et al., 2001): where e ! r is the unit vector in the radial direction and the local acceleration of gravity is assumed to vary only along the radial coordinate. In consistent with previous studies on convections of constant fluid density (e.g., Schubert et al., 2001), we assume in this paper that the acceleration of gravity is proportional to the radial coordinate and relates to the acceleration of gravity at the outer boundary g 0 by Since the fluid velocity components and the temperature perturbation are all small quantities at the onset of convection, the above momentum and energy equations do not contain nonlinear terms, and the governing equations are all linear.
For the linear stability analysis in this paper, the radial velocity is always assumed to be zero at both the outer and inner boundaries of the spherical shell, resulting in the following two impermeability conditions: As for the thermal boundary conditions, either a constant heat flux or a constant temperature is applied at each boundary, resulting in the following conditions: 1. constant heat flux on the outer boundary (denoted by H) 2. constant temperature on the outer boundary (T) Similarly, the lateral kinematic boundary conditions are assumed to be either free-slip (zero shear stress) or fixed, resulting in the following conditions: 1. free-slip outer boundary (F) The solution procedure for the above defined linear stability problem of thermal convection in a spherical shell is detailed in Appendices A-C, while FORTRAN codes are provided in the supporting information.

Thermal Convections in a Sphere (η = 0)
We first present the results for thermal convections in a sphere, which corresponds to a spherical shell with a radius ratio of η = 0. Since there is no inner boundary for a sphere, all the boundary conditions are specified through the outer boundary and internal heating is the only heating mode to be considered. For four combinations of the boundary conditions, the obtained critical Rayleigh numbers and their associated spherical harmonic degrees are listed in Table 1.
The results listed in Table 1 agree well with previous results available in the literature. Zebib et al. (1983) obtained a critical Rayleigh number of Ra H = 9,276.6 for the TF boundary conditions, while our result ofRa ih = 1,545.596363 can be converted to Ra H = 9,273.578181 according to equation (14). Backus (1955) obtained critical Rayleigh numbers of 3,091.180 ± 0.014 and 8,040.07 ± 0.04, respectively, for  (1955) is exactly twice the value of our definition, the corresponding critical Rayleigh numbers of our results can be converted to be 3,091.192726 and 8,040.046890, respectively, in accordance with his definition. Since our method does not need to calculate the Bessel function, the series solution can be easily constructed to include more high-order terms and provide higher accuracy. Table 1 shows that all the four cases result in a degree-1 thermal convection, with the critical Rayleigh number varying from around 480 for the HF case to over 4,000 for the TX case. The variation of the critical Rayleigh number indicates different constraining strengths posed by the boundary conditions. Table 1 suggests that an isothermal boundary is stronger than a constant heat flux boundary and a fixed boundary is stronger than a free-slip boundary for preventing convection to occur. This is understandable, considering that no lateral temperature gradient is allowed to be built up on an isothermal boundary and that frictional resistance has to be overcome on a fixed boundary, both hindering thermal convection to occur. Furthermore, that the HX case requires a higher critical Rayleigh number than the TF case implies that an isothermal condition is a weaker constraint than a fixed kinematic condition for thermal convection onset in a sphere.
The infinitesimal-amplitude convective flow and perturbation temperature fields for the onset-stage thermal convections in a sphere subject to internal heating are compared in Figure 1 for different boundary conditions. As the conductive temperature (not shown in the figure) is spherically symmetric and decreases monotonically from the spherical center to the boundary, the upward streams traveling from the lower pole to the upper pole are heated through the journey passing the central hot region, resulting in a hotter upper hemisphere and a cooler lower hemisphere. In spite of such a common feature for all the cases shown in Figure 1, the figure clearly shows the distinct effects of different boundary conditions. For the HF case (Figure 1a), both the maximum and minimum perturbation temperature regions appear at the polar regions, while a switch to a fixed boundary (Figure 1b) prevents the convective flow reaching the boundary, causing the extreme perturbation temperature regions inwardly moving a distance from the polar regions. In contrast to the constant heat flux boundary cases (Figures 1a and 1b), the extreme perturbation temperature regions associated with the isothermal boundary cases (Figures 1c and 1d) lie halfway between the spherical center and the poles.

Thermal Convections in Spherical Shells (η > 0)
For spherical shells of η > 0, various combinations of boundary conditions on both the outer and inner boundaries result in 16 cases. Table 2 lists the obtained critical Rayleigh numbers and the associated spherical harmonic degrees for spherical shells of two radius ratios (η = 0.5 and η = 0.95) subject to basal heating (ε = 0) and internal heating (ε = 1) under different boundary conditions.
The critical Rayleigh numbers listed in Table 2 agree well with previous results available in the literature.
For an η = 0.5 spherical shell subject to internal heating under the ThFf boundary conditions, Schubert and Zebib (1980) obtained Ra H = 2,174.5, and Zebib et al. (1983) improved it to 2,174.3, while our result is Ra ih = 724.784653, corresponding toRa H = 2,174.353959 according to equation (14). For the basal heating case of an η = 0.5 spherical shell under the TtFf boundary conditions, the critical Rayleigh number of 978.5 obtained by Zebib et al. (1983) is also in good agreement with our result of 978.558336.
In Table 2, the 16 boundary condition combinations are divided into four groups based on thermal conditions, with each group including four cases associated with different kinematic conditions. It is interesting to note that all the cases in the Hh group are associated with degree-1 thermal convection, suggesting that lateral temperature gradients on both boundaries provide favorable conditions for forming degree-1 convection.  Table 2 indicates that both the critical Rayleigh number and its associated harmonic degree almost always increase along with the intergroup sequence Hh-Ht-Th-Tt, as well as with the intragroup sequence of Ff-Fx-Xf-Xx. Such sequences reflect the increasing constraining strength posed by the boundary conditions. Weaker constraints posed by the Ht and Fx conditions than the Th and Xf conditions reveals the more important role of the outer boundary than the inner boundary, since the outer boundary always has a larger surface area than the inner boundary. The only exceptional case that deviates from the sequences is the ThXf case of the η = 0.95 shell, which is associated with a lower critical Rayleigh number under the basal heating mode and lower harmonic degrees under both the basal heating and internal heating modes, as compared with the counterparts of the ThFx and HtXf cases. For the η = 0.95 shell, the difference in area between the outer boundary and inner boundary is insignificant, while the favorable conditions for basal heating convection and unfavorable conditions for internal heating convection provided by the ThXf case may explain the exceptional behavior associated with this case.

Earth and Space Science
The onset-stage convective flow and perturbation temperature fields in a η = 0.5 shell are shown in Figure 2 for basal heating mode and in Figure 3 for internal heating mode, in terms of axisymmetric flow cases characterized by degree-1 or degree-2 convection. The effects of boundary conditions on the flow and temperature fields are clearly exhibited by the figures.
A comparison between the basal heating columns and the internal heating columns in Table 2 indicates that the critical Rayleigh number associated with basal heating is generally greater than that associated with internal heating for the same shell geometry and boundary conditions. This general trend may be explained by the conductive temperature gradient that causes the onset of thermal convection. According to equation (2), the conductive temperature gradient under basal heating decreases gradually from the inner boundary to the outer boundary, whereas that under internal heating constantly increases outwardly from zero on the inner boundary. Such a difference results in greater buoyancy body forces near the inner boundary for the basal heating mode and near the outer boundary for the internal heating mode. This can explain the major difference between Figures 2 and 3; that is, the extreme perturbation temperature regions are distributed closer to the outer boundary in Figure 3 than those in Figure 2. Considering that the outer boundary has a larger area than the inner boundary, the buoyancy forces in the internal heating mode normally have more profound effects than the basal heating mode, in accordance with smaller critical Rayleigh numbers being required for thermal convections subject to internal heating.
However, such a normal trend may be broken if the boundary conditions are specified to favor inner boundary flow but hinder outer boundary flow. The most remarkable case is the ThXf case, which favors near inner boundary flow by constant heat flux and free-slip conditions and hinders near outer boundary flow by isothermal and fixed conditions. Table 2 indicates that the ThXf case is indeed the most abnormal case, which is associated with the largest difference of the critical Rayleigh number under the internal heating mode over that under the basal heating mode for both radius ratios. It is also evident in Table 2 that the HhXf and ThFf cases are abnormal for both radius ratios, and the ThXx and TtXf cases are abnormal for radius ratio of 0.95. This is understandable, as these four cases are the only cases that share three out of the four conditions with the ThXf case.
For combined heating, the variations of the critical Rayleigh numbers Ra ch , Ra bh , and Ra ih and the associated harmonic degree along with continuously varying internal heating fraction ε are demonstrated in Figure 4, in terms of thermal convections in an η = 0.25 spherical shell under the TtFf boundary conditions as an example. It can be observed that although all the three critical Rayleigh numbers vary monotonically and smoothly with the internal heating fraction, there exists a transition from degree-1 convection to degree-2 convection when increasing the value of ε from 0.451 to 0.452. The vertical line going through the asterisk in Figure 4 represents this transition from l = 1 to l = 2, and the convective flow and perturbation temperature fields before and after the transition are compared in Figure 5, together with the corresponding critical Rayleigh numbers.
In addition to systematically comparing the effects of boundary conditions and heating modes, Table 2 also provides a comparison for the effects of shell geometry. Under the same boundary conditions and heating modes, it is evident in Table 2 that the critical Rayleigh number decreases with the increasing inner-outer radius ratio from 0.5 to 0.95 for all the cases, except for the internal heating cases under the HtFf and HtFx boundary conditions. The general trend of the dependence of the critical Rayleigh

10.1029/2019EA000672
Earth and Space Science number on radius ratio reflects the curvature effects on preventing thermal convection to occur. For a same shell thickness, a spherical shell with a greater curvature has stronger constraint for convective flow, normally requiring a higher temperature difference across the shell thickness, leading to a higher critical Rayleigh number.
In order to thoroughly investigate the shell curvature effects, the variations of the critical Rayleigh numbers and the associated harmonic degrees along with continuously varying radius ratio from 0.01 to 0.95 are demonstrated in Figure 6 for the TtFf case as an example for the normal trend and for the HtFx case as an abnormal example. As shown in Figure 6a, the critical Rayleigh numbers for both the basal heating and internal heating modes decrease monotonically with increasing radius ratio, while the slopes of the curves exhibit discontinuous transitions along with harmonic degree changes, which is remarkably different from the smooth transition with internal heating fraction change shown in Figure 4. The difference between the two curves in Figure 6a tends to decrease with increasing radius ratio, reflecting that the difference between basal heating and internal heating is more profound in highly curved shells.
The abnormal HtFx case shown in Figure 6b exhibits distinctly features, as compared with Figure 6a. Although both the two curves in Figure 6b are not strictly monotonic, the basal heating curve shows a general trend that is similar to the normal trend shown in Figure 6a, while the internal heating curve shows a rather complicated trend, with the critical Rayleigh number even increasing with the radius ratio at η > 0.6. Such an abnormal behavior is a consequence of the specified boundary conditions. Contrary to the ThXf case, the HtFx case favors near outer boundary flow and hinders near inner boundary flow, so it is especially suitable for developing internal heating convection, which accumulates buoyancy forces near the outer boundary region. Along with an increase in radius ratio, the resistance due to the shell curvature to hinder convection decreases, and the driving forces from the accumulative buoyancy to develop convection also decrease. The competition between these two effects results in the complicated and abnormal trend in the internal heating curve shown in Figure 6b.

Conclusions
In this paper, we propose a unified solution to linear stability analysis for thermal convections in spherical shells, and systematically investigate the effects of boundary conditions, heating modes, and shell geometry. The general trends for dependence of critical Rayleigh numbers and their associated harmonic degrees upon boundary condition type, internal heating fraction, and shell radius ratio reveal the fundamental physical mechanisms affecting thermal convections, while the remarkable exceptional cases deviating from the general trends clearly demonstrate that thermal convections are consequences of complicated interplay between multiple competing effects.
The explicit formulas given in Appendix C for calculating the elements of the vanishing determinant provide a convenient tool for linear stability analysis of thermal convections in spheres and spherical shells. Although only extreme boundary conditions (i.e., constant heat flux or isothermal for thermal conditions and free-slip or fixed for kinematic conditions) are discussed in this paper, there is no significant difficulty in principle in modifying the formulas in Appendix C for solving problems subject to more realistic boundary constraints (e.g., convective boundary condition).

Appendix A: Solution Procedure of Linear Stability Analysis
We start from the linear equations (15)-(17). The continuity equation for incompressible fluid, equation (15), can be automatically satisfied if the solenoidal velocity field is represented by the sum of poloidal and toroidal vector fields (Chandrasekhar, 1961). In the problem of thermal convection onset, the toroidal function is identically zero as buoyancy forces only induce a poloidal velocity field (Schubert et al., 2001;Zebib et al., 1985). As a result, the velocity components in spherical coordinates can be written only in terms of a poloidal function Φ by where the operator L 2 is defined as and its eigenfunctions are spherical harmonics Y m l θ; φ ð Þ of degree l and order m, that is, The Laplacian is related to the operator L 2 by By taking the curl of equation (16), the pressure term in the momentum equation can be eliminated and results in Figure 6. Variations of the critical Rayleigh numbers and harmonic degrees along with the radius ratio η of spherical shells for thermal convections due to different heating modes under the (a) TtFf and (b) HtFx boundary conditions. The Ra bh and Ra ih curves correspond to the basal heating (ε = 0) and internal heating (ε = 1) cases, respectively, while the combined heating cases lie between the two extreme cases. Solid lines connecting two symbols show transitions for the critical harmonic degrees, that is, the line connecting the two asterisks (labeled 1-2) represents the transition from l = 1 to l = 2.

10.1029/2019EA000672
Earth and Space Science The ϕ component of this equation yields Substituting equation (A1) into equation (17), the energy equation reduces to The problem of convective instability onset in a sphere or a spherical shell reduces to the solution of equations (A3b) and (A3c). The solution can be most readily obtained by an expansion of the perturbation temperature field and the perturbation poloidal velocity field in terms of spherical harmonics, that is, Substituting equations (A4) and (A5) into equations (A3b) and (A3c) yields Combining equations (A6a) and (A6b) and invoking equations (18) and (2), respectively, to evaluate the acceleration of gravity and the conductive temperature gradient results in a single ordinary differential equation of six orders for the poloidal function and The ratio between these two Rayleigh numbers is 10.1029/2019EA000672

Earth and Space Science
which depends only upon the internal heating fraction defined in equation (7) and the radius ratio defined in equation (8). Invoking equation (4), the classical Rayleigh number for combined heating can be directly related to Ra 1 and Ra 2 by The central task for the linear stability analysis is solving equation (A9) subject to different boundary conditions. The boundary conditions proposed in equations (19)-(24a) reduce to the following conditions for b All these conditions are homogenous. Although the thermal boundary conditions proposed in equations (22a)-(23b) are nonhomogenous, the corresponding conditions for perturbation temperature and its gradient are homogenous.
Recognizing that b r ¼ 0 is a regular singular point (Asmar, 2004) of equation (A9), we can expand the nondimensional poloidal function into Frobenius series (Asmar, 2004, p. A53), that is, where s represents the lowest power in the series and z n are expansion coefficients, and substituting it into equation (A9) gives The coefficient of the leading term b r s−6 in equation (A14) must vanish, resulting in According to the Frobenius method (Asmar, 2004, pp. A51-A62; Courant & Hilbert, 1989, pp. 492-501), the six roots for s in equation (A15) define six linearly independent fundamental solutions of equation (A9), which can be written by where The general solution of equation (A9) is a combination of the six fundamental solutions, that is, Substituting equation (A21) into the six boundary conditions given by equations (A12a)-(A12f) results six algebraic equations that can be written in matrix form as where B is a 6 × 6 square matrix and its elements are provided in Appendix C.
Since equation (A9) and boundary conditions (A12a)-(A12f) are all homogeneous, nontrivial solution for equation (A22) exists only when the determinant of matrix B vanishes. Therefore, solving the linear stability problem reduces to search for the minimum Rayleigh number (defined as the critical Rayleigh number) that makes the determinant of matrix B be zero by varying l value for specified radius ratio and internal heating fraction under various boundary conditions. In order to ensure that all the involved computations are sufficient to provide critical Rayleigh numbers accurate to at least six decimal places, we use quad-precision (about 36 significant digits) computations for spheres and spherical shells with a radius ratio less than 0.7, and invoke a thread-safe arbitrary precision package MPFUN-Fort (Bailey, 2017;Bailey et al., 2002) to employ up to 1,000 significant digits in computations for thin shells. The relevant FORTRAN codes are provided in the supporting information.
It should be noted that the critical Rayleigh number is determined together with its associated spherical harmonic degree l, while the harmonic order m for the thermal convection cannot be determined from linear stability analysis, since equation (9) is independent of m. The determination of preferred value of m requires a nonlinear stability analysis of the convective problem. It has been proven (Busse, 1975(Busse, , 1978Busse & Riahi, 1982) that for degree-1 and degree-2 thermal convections at the onset stage, stable solutions corresponding to different m values are identical and axisymmetric, whereas the most stable solutions for l > 2 convections are in general three-dimensional and asymmetric.