  We develop a $ P $-multigrid solver to simulate locally preconditioned
unsteady compressible Navier-Stokes equations at low Mach numbers with implicit
high-order methods. Specifically, the high-order flux reconstruction/correction
procedure via reconstruction (FR/CPR) method is employed for spatial
discretization and the high-order time integration is conducted by means of the
explicit first stage, singly diagonally implicit Runge-Kutta (ESDIRK) method.
Local preconditioning is used to alleviate the stiffness of the compressible
Navier-Stokes equations at low Mach numbers and is only conducted in pseudo
transient continuation to ensure the high-order accuracy of ESDIRK methods. We
employ the element Jacobi smoother to update the solutions at different $ P
$-levels in the $ P $-multigrid solver. High-order spatiotemporal accuracy of
the new solver for low-Mach-number flow simulation is verified with the
isentropic vortex propagation when the Mach (Ma) number of the free stream is
0.005. The impact of the hierarchy of polynomial degrees on the convergence
speed of the $ P $-multigrid method is studied via several numerical
experiments, including two dimensional (2D) inviscid and viscous flows over a
NACA0012 airfoil at $\text{Ma} = 0.001$, and a three dimensional (3D) inviscid
flow over a sphere at $\text{Ma} = 0.001$. The $ P $-multigrid solver is then
applied to coarse resolution simulation of the transitional flows over an
SD7003 wing at $ 8^\circ $ angle of attack when the Reynolds number is 60000
and the Mach number is 0.1 or 0.01.
