  We calculate the next-to-next-to-leading order QCD corrections to the
leptonic invariant mass ($q^2$) spectrum of semileptonic $b \to c$ inclusive
decays, taking into account the mass of the charm quark and the charged lepton
in the final state. We obtain analytic results in terms of generalized
polylogarithms and present numerical studies of the $\mathcal{O}(\alpha_s^2)$
corrections to the $q^2$ spectrum of $b \to c \ell \bar \nu_\ell$ decays, for
$\ell =e, \mu$ and $\tau$, in the kinetic scheme. Our computation can be used
to incorporate the recent measurements of $q^2$ moments by Belle and Belle II
into global fits of inclusive semileptonic $B$-decays.
