Density-based solvers, originally designed for high Mach number flows, suffer from a decrease in solution quality and convergence rate as the Mach number decreases. Low Mach preconditioning has been shown by previous investigations [1] to improve the efficiency of density-based solvers at low Mach numbers in steady-state computations. Applying a similar procedure, this study presents the implementation of low Mach preconditioning in harmonic balance solvers to enhance the solution quality and optimize the convergence rate of unsteady frequency-domain computations in the low Mach regime. The convergence rate is improved by multiplying the time derivative of the equations with a preconditioning matrix, to equalizes the acoustic and convective propagation velocities. However, for unsteady equations, the preconditioner introduces a frequency dependence in the acoustic propagation velocity. Thus, to optimize the convergence rate in the Harmonic Balance solver, each harmonic requires its own frequency dependent preconditioner. Inaccuracies at low Mach numbers occur due to the excessive numerical dissipation of the roe solver. In this study we utilize a preconditioned Roe solver for unsteady computations [2], which ensures the correct Mach number scaling in the artificial diffusion, significantly improving the solution quality. Using the low Mach preconditioned Harmonic Balance solver, we computed a laminar cylinder and found that an inappropriate preconditioner can reduce the convergence rate. Optimal convergence is attained by employing the steady-state preconditioner for the zeroth harmonic and applying the frequency-dependent preconditioner for higher harmonics. In addition, the study highlights the shortcomings of the standard Roe solver in predicting the solution’s unsteadiness for decreasing Mach numbers. In contrast, the preconditioned Roe solver produces Mach number independent solutions, establishing its superior predictive capability.