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.