  We study the halo mass function (HMF) in modified gravity (MG) models using a
set of large $N$-body simulations -- the ELEPHANT suite. We consider two
popular beyond-general relativity scenarios: the Hu-Sawicki chameleon $f(R)$
model and the normal branch of the Dvali-Gabadadze-Porrati braneworld (nDGP).
We show that in MG, analytic formulation based on the Press-Schechter framework
offers a grossly inaccurate description of the HMF. We find, however, that once
the HMF is expressed in terms of the dimensionless multiplicity function, it
approximately assumes a redshift-independent universal character for all the
models. Exploiting this property, we propose universal fits for the MG HMF in
terms of their fractional departures from the $\mathrm{\Lambda \text{CDM}}$
case. We find two enclosed formulas, one for $f(R)$ and another for nDGP, that
provide a reliable description of the HMF over the mass range covered by the
simulations. These are accurate to a few percent with respect to the $N$-body
data. We test the extrapolation potential of our fits against separate
simulations with a different cosmological background and mass resolution, and
find very good accuracy, within $\sim 10\%$. A particularly interesting finding
from our analysis is a Gaussian-like shape of the HMF deviation that seems to
appear universally across the whole $f(R)$ family, peaking at a mass variance
scale characteristic for each $f(R)$ variant. We attribute this behavior to the
specific physics of the environmentally-dependent chameleon screening models.
