A search for dark matter particles is performed by looking for events with large transverse momentum imbalance and a recoiling Higgs boson decaying to either a pair of photons or a pair of $\tau$ leptons. The search is based on proton-proton collision data at a center-of-mass energy of 13 TeV collected at the CERN LHC in 2016 and corresponding to an integrated luminosity of 35.9 fb$^{-1}$. No significant excess over the expected standard model background is observed. Upper limits at 95% confidence level are presented for the product of the production cross section and branching fraction in the context of two benchmark simplified models. For the Z'-two-Higgs-doublet model (where Z' is a new massive boson mediator) with an intermediate heavy pseudoscalar particle of mass $m_\mathrm{A} =$ 300 GeV and $m_\mathrm{DM} =$ 100 GeV, Z' masses from 550 GeV up to 1265 GeV are excluded. For a baryonic Z' model, with $m_\mathrm{DM} =$ 1 GeV, Z' masses up to 615 GeV are excluded. Results are also presented for the spin-independent cross section for the dark matter-nucleon interaction as a function of the mass of the dark matter particle. This is the first search for dark matter particles produced in association with a Higgs boson decaying to two $\tau$ leptons.