This paper presents a search for decays of the Higgs boson with a mass of 125 GeV to a pair of new spin-0 particles, $H\rightarrow aa$, where one $a$-boson decays to a $b$-quark pair and the other to a muon pair. The search uses 139 fb$^{-1}$ of the LHC proton--proton collision data at a center-of-mass energy of $\sqrt{s}$=13 TeV recorded by the ATLAS experiment between 2015 and 2018. A narrow dimuon resonance is searched for in the invariant mass spectrum between 16 GeV and 62 GeV. The largest excess of events above the Standard Model backgrounds is observed at an invariant dimuon mass of 52 GeV and corresponds to a local (global) significance of $3.3 \sigma$ ($1.7 \sigma$). Upper limits at 95% confidence level are placed on the branching ratio of the Higgs boson to the $bb\mu\mu$ final state, $ \mathcal{B}(H\rightarrow aa\rightarrow bb\mu\mu)$, and are in the range $(0.22-4.0) \times10^{-4}$, depending on the signal mass hypothesis. The events selected for the analysis are required to have two muons of opposite charge with the leading and sub-leading muons satisfying either $p_T^{\text{leading}}>27$ GeV and $p_T^{\text{subleading}}>5$ GeV, and the event being triggered with the single muon triggers, or with both muons having $p_T>15$ GeV, and the event being triggered with the dimuon triggers. The dimuon invariant mass is required to be between 15 GeV and 65 GeV. Furthermore, the events have to contain exactly two $b$-tagged jets with a $p_T$ above 20 GeV.

Post-fit number of background events in all SR bins without applying the BDT cuts that are tested for the presence of signal. The bins are 2 GeV (3 GeV) wide in mµµ for $m_a$ ≤ 45 GeV ($m_a$ > 45 GeV). Events in neighbouring bins partially overlap. Discontinuities in the background predictions appear when the BDT discriminant used for the selection changes from the one trained in the lower mass range to the one trained in the higher mass range.