  A search for exotic decays of the Higgs boson (H) with a mass of 125 GeV to a
pair of light pseudoscalars $\mathrm{a}_1$ is performed in final states where
one pseudoscalar decays to two b quarks and the other to a pair of muons or
$\tau$ leptons. A data sample of proton-proton collisions at $\sqrt{s}$ = 13
TeV corresponding to an integrated luminosity of 138 fb$^{-1}$ recorded with
the CMS detector is analyzed. No statistically significant excess is observed
over the standard model backgrounds. Upper limits are set at 95% confidence
level (CL) on the Higgs boson branching fraction to $\mu\mu$bb and to
$\tau\tau$bb, via a pair of $\mathrm{a}_1$s. The limits depend on the
pseudoscalar mass $m_{\mathrm{a}_1}$ and are observed to be in the range
(0.17-3.3) $\times$ 10$^{-4}$ and (1.7-7.7) $\times$ 10$^{-2}$ in the
$\mu\mu$bb and $\tau\tau$bb final states, respectively. In the framework of
models with two Higgs doublets and a complex scalar singlet (2HDM+S), the
results of the two final states are combined to determine model-independent
upper limits on the branching fraction $\mathcal{B}$(H $\to$
$\mathrm{a}_1\mathrm{a}_1$ $\to$ $\ell\ell$bb) at 95% CL, with $\ell$ being a
muon or a $\tau$ lepton. For different types of 2HDM+S, upper bounds on the
branching fraction $\mathcal{B}$(H $\to$ $\mathrm{a}_1\mathrm{a}_1$) are
extracted from the combination of the two channels. In most of the Type II
2HDM+S parameter space, $\mathcal{B}($H $\to$ $\mathrm{a}_1\mathrm{a}_1$)
values above 0.23 are excluded at 95% CL for $m_{\mathrm{a}_1}$ values between
15 and 60 GeV.
