  A search for a light pseudoscalar Higgs boson (a) decaying from the 125 GeV
(or a heavier) scalar Higgs boson (H) is performed using the 2016 LHC
proton-proton collision data at $\sqrt{s} =$ 13 TeV, corresponding to an
integrated luminosity of 35.9 fb$^{-1}$, collected by the CMS experiment. The
analysis considers gluon fusion and vector boson fusion production of the H,
followed by the decay H $\to$ aa $\to$ $\mu\mu\tau\tau$, and considers
pseudoscalar masses in the range 3.6 $\lt$ $m_\mathrm{a}$ $\lt$ 21 GeV. Because
of the large mass difference between the H and the a bosons and the small
masses of the a boson decay products, both the $\mu\mu$ and the $\tau\tau$
pairs have high Lorentz boost and are collimated. The $\tau\tau$ reconstruction
efficiency is increased by modifying the standard technique for hadronic $\tau$
lepton decay reconstruction to account for a nearby muon. No significant signal
is observed. Model-independent limits are set at 95% confidence level, as a
function of $m_\mathrm{a}$, on the branching fraction ($\mathcal{B}$) for H
$\to$ aa $\to$ $\mu\mu\tau\tau$, down to 1.5 (2.0) $\times$ 10$^{-4}$ for
$m_\mathrm{H} =$ 125 (300) GeV. Model-dependent limits on $\mathcal{B}$(H $\to$
aa) are set within the context of two Higgs doublets plus singlet models, with
the most stringent results obtained for Type-III models. These results extend
current LHC searches for heavier a bosons that decay to resolved lepton pairs
and provide the first such bounds for an H boson with a mass above 125 GeV.
