Abstract: 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$$ s = 13 TeV with the ATLAS experiment at the Large Hadron Collider using data corresponding to an integrated luminosity of 140 $$\hbox {fb}^{-1}$$ 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_\text {T} > 20$$ p T > 20 GeV and rapidity $$|y| < 2.5$$ | 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$$ Z + ≥ 1 c -jet observables, optimized for sensitivity to intrinsic-charm, are compared with benchmark models with different intrinsic-charm fractions.