  The production of dark matter in association with Higgs bosons is predicted
in several extensions of the Standard Model. An exploration of such scenarios
is presented, considering final states with missing transverse momentum and
$b$-tagged jets consistent with a Higgs boson. The analysis uses proton-proton
collision data at a centre-of-mass energy of 13 TeV recorded by the ATLAS
experiment at the LHC during Run 2, amounting to an integrated luminosity of
139 fb$^{-1}$. The analysis, when compared with previous searches, benefits
from a larger dataset, but also has further improvements providing sensitivity
to a wider spectrum of signal scenarios. These improvements include both an
optimised event selection and advances in the object identification, such as
the use of the likelihood-based significance of the missing transverse momentum
and variable-radius track-jets. No significant deviation from Standard Model
expectations is observed. Limits are set, at 95% confidence level, in two
benchmark models with two Higgs doublets extended by either a heavy vector
boson $Z'$ or a pseudoscalar singlet $a$ and which both provide a dark matter
candidate $\chi$. In the case of the two-Higgs-doublet model with an additional
vector boson $Z'$, the observed limits extend up to a $Z'$ mass of 3 TeV for a
mass of 100 GeV for the dark matter candidate. The two-Higgs-doublet model with
a dark matter particle mass of 10 GeV and an additional pseudoscalar $a$ is
excluded for masses of the $a$ up to 520 GeV and 240 GeV for $\tan \beta = 1$
and $\tan \beta = 10$ respectively. Limits on the visible cross-sections are
set and range from 0.05 fb to 3.26 fb, depending on the missing transverse
momentum and $b$-quark jet multiplicity requirements.
