  We provide a systematic renormalization group formalism for the mass effects
in the relation of the pole mass $m_Q^{\rm pole}$ and short-distance masses
such as the $\overline{\rm MS}$ mass $\overline{m}_Q$ of a heavy quark $Q$,
coming from virtual loop insertions of massive quarks lighter than $Q$. The
formalism reflects the constraints from heavy quark symmetry and entails a
combined matching and evolution procedure that allows to disentangle and
successively integrate out the corrections coming from the lighter massive
quarks and the momentum regions between them and to precisely control the large
order asymptotic behavior. With the formalism we systematically sum logarithms
of ratios of the lighter quark masses and $m_Q$, relate the QCD corrections for
different external heavy quarks to each other, predict the ${\cal
O}(\alpha_s^4)$ virtual quark mass corrections in the pole-$\overline{\rm MS}$
mass relation, calculate the pole mass differences for the top, bottom and
charm quarks with a precision of around $20$ MeV and analyze the decoupling of
the lighter massive quark flavors at large orders. The summation of logarithms
is most relevant for the top quark pole mass $m_t^{\rm pole}$, where the
hierarchy to the bottom and charm quarks is large. We determine the ambiguity
of the pole mass for top, bottom and charm quarks in different scenarios with
massive or massless bottom and charm quarks in a way consistent with heavy
quark symmetry, and we find that it is $250$ MeV. The ambiguity is larger than
current projections for the precision of top quark mass measurements in the
high-luminosity phase of the LHC.
