  Stellar convection poses two main gargantuan challenges for astrophysical
fluid solvers: low-Mach number flows and minuscule perturbations over steeply
stratified hydrostatic equilibria. Most methods exhibit excessive numerical
diffusion and are unable to capture the correct solution due to large
truncation errors. In this paper, we analyze the performance of the Spectral
Difference (SD) method under these extreme conditions using an arbitrarily
high-order shock capturing scheme with a posteriori limiting. We include both a
modification to the HLLC Riemann solver adapted to low Mach number flows
(L-HLLC) and a well-balanced scheme to properly evolve perturbations over steep
equilibrium solutions. We evaluate the performance of our method using a series
of test tailored specifically for stellar convection. We observe that our
high-order SD method is capable of dealing with very subsonic flows without
necessarily using the modified Riemann solver. We find however that the
well-balanced framework is unavoidable if one wants to capture accurately small
amplitude convective and acoustic modes. Analyzing the temporal and spatial
evolution of the turbulent kinetic energy, we show that our fourth-order SD
scheme seems to emerge as an optimal variant to solve this difficult numerical
problem.
