  This paper presents a measurement of the production cross-section of a $Z$
boson in association with $b$- or $c$-jets, in proton-proton collisions at
$\sqrt{s} = 13$ TeV with the ATLAS experiment at the Large Hadron Collider
using data corresponding to an integrated luminosity of 140 fb$^{-1}$.
Inclusive and differential cross-sections are measured for events containing a
$Z$ boson decaying into electrons or muons and produced in association with at
least one $b$-jet, at least one $c$-jet, or at least two $b$-jets with
transverse momentum $p_\textrm{T} > 20$ GeV and rapidity $|y| < 2.5$.
Predictions from several Monte Carlo generators based on next-to-leading-order
matrix elements interfaced with a parton-shower simulation, with different
choices of flavour schemes for initial-state partons, are compared with the
measured cross-sections. The results are also compared with novel predictions,
based on infrared and collinear safe jet flavour dressing algorithms. Selected
$Z + \ge 1 c$-jet observables, optimized for sensitivity to intrinsic-charm,
are compared with benchmark models with different intrinsic-charm fractions.
