A search for exotic decays of the Higgs boson ($$\text {H}$$$\text{H}$) with a mass of 125$$\,\text {Ge}\hspace{-.08em}\text {V}$$$\text{Ge}\text{V}$to a pair of light pseudoscalars$$\text {a}_{1} $$${\text{a}}_1$is performed in final states where one pseudoscalar decays to two$${\textrm{b}}$$$\text{b}$quarks and the other to a pair of muons or$$\tau $$$\tau$leptons. A data sample of proton–proton collisions at$$\sqrt{s}=13\,\text {Te}\hspace{-.08em}\text {V} $$$\sqrt{s}=13\text{Te}\text{V}$corresponding to an integrated luminosity of 138$$\,\text {fb}^{-1}$$${\text{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 ($$\text {CL}$$$\text{CL}$) on the Higgs boson branching fraction to$$\upmu \upmu \text{ b } \text{ b } $$$\mu \mu \text{b}\text{b}$and to$$\uptau \uptau \text{ b } \text{ b },$$$\tau \tau \text{b}\text{b},$via a pair of$$\text {a}_{1} $$${\text{a}}_1$s. The limits depend on the pseudoscalar mass$$m_{\text {a}_{1}}$$$m_{{\text{a}}_1}$and are observed to be in the range (0.17–3.3) $$\times 10^{-4}$$$\times {10}^{-4}$and (1.7–7.7) $$\times 10^{-2}$$$\times {10}^{-2}$in the$$\upmu \upmu \text{ b } \text{ b } $$$\mu \mu \text{b}\text{b}$and$$\uptau \uptau \text{ b } \text{ b } $$$\tau \tau \text{b}\text{b}$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 upper limits on the branching fraction$${\mathcal {B}}(\text {H} \rightarrow \text {a}_{1} \text {a}_{1} \rightarrow \ell \ell \text{ b } \text{ b})$$$B(\text{H}\to {\text{a}}_1{\text{a}}_1\to \ell \ell \text{b}\text{b})$at 95%$$\text {CL}$$$\text{CL}$, with$$\ell $$$\ell$being a muon or a$$\uptau $$$\tau$lepton. For different types of 2HDM+S, upper bounds on the branching fraction$${\mathcal {B}}(\text {H} \rightarrow \text {a}_{1} \text {a}_{1} )$$$B(\text{H}\to {\text{a}}_1{\text{a}}_1)$are extracted from the combination of the two channels. In most of the Type II 2HDM+S parameter space,$${\mathcal {B}}(\text {H} \rightarrow \text {a}_{1} \text {a}_{1} )$$$B(\text{H}\to {\text{a}}_1{\text{a}}_1)$values above 0.23 are excluded at 95%$$\text {CL}$$$\text{CL}$for$$m_{\text {a}_{1}}$$$m_{{\text{a}}_1}$values between 15 and 60$$\,\text {Ge}\hspace{-.08em}\text {V}$$$\text{Ge}\text{V}$.