Abstract: With the rapid advancement of cutting-edge technologies such as hypersonic vehicles operating in the upper atmosphere, micro-electro-mechanical systems (MEMS), and extreme ultraviolet (EUV) lithography machines, increasing attention has been directed toward the analysis of flow stability under the coupled influence of rarefied gas effects and thermal convection. In such high-tech applications, traditional flow stability analysis methods based on the Navier–Stokes equations fall short, as they are not suitable for accurately describing flows at higher levels of rarefaction. Meanwhile, the direct simulation Monte Carlo (DSMC) method, although capable of handling rarefied flows, suffers from high computational costs and significant statistical fluctuations, which make it unsuitable for accurately tracking the evolution of small disturbances. Furthermore, the widely used Boltzmann-BGK (Bhatnagar-Gross-Krook) model equation is limited by its assumption of a fixed Prandtl number equal to one, which restricts its applicability in real-world situations where thermophysical properties vary with conditions. To overcome these challenges, a linear stability analysis framework based on the Boltzmann-Shakhov model equation is proposed. This model accommodates a tunable Prandtl number and thus better captures the physical characteristics of different flow regimes. Based on the assumption of small perturbations and using a modal analysis approach, a full-domain linear stability equation is derived. This enables a comprehensive assessment of flow stability that accounts for both rarefied gas dynamics and thermal effects. Applying this method to the Rayleigh–Bénard convection problem, the study shows that the flow becomes highly unstable when the density distribution between the upper and lower walls shifts from a configuration of “lighter above, heavier below” to “heavier above, lighter below.” A noticeable turning point in the disturbance growth rate appears near a specific wavenumber, likely due to the eigenfunction of the disturbance changing from having one extremum to having two. Moreover, as rarefaction increases, the Prandtl number corresponding to the most unstable mode gradually decreases. These findings shed light on the intricate interplay between rarefied gas behavior and thermal convection in determining flow stability.