Abstract Cross-sections for the production of a Z boson in association with two photons are measured in proton–proton collisions at a centre-of-mass energy of 13 TeV. The data used correspond to an integrated luminosity of 139 fb $$^{-1}$$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msup> <mml:mrow /> <mml:mrow> <mml:mo>-</mml:mo> <mml:mn>1</mml:mn> </mml:mrow> </mml:msup> </mml:math> recorded by the ATLAS experiment during Run 2 of the LHC. The measurements use the electron and muon decay channels of the Z boson, and a fiducial phase-space region where the photons are not radiated from the leptons. The integrated $$Z(\rightarrow \ell \ell )\gamma \gamma $$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mrow> <mml:mi>Z</mml:mi> <mml:mo>(</mml:mo> <mml:mo>→</mml:mo> <mml:mi>ℓ</mml:mi> <mml:mi>ℓ</mml:mi> <mml:mo>)</mml:mo> <mml:mi>γ</mml:mi> <mml:mi>γ</mml:mi> </mml:mrow> </mml:math> cross-section is measured with a precision of 12% and differential cross-sections are measured as a function of six kinematic variables of the $$Z\gamma \gamma $$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mrow> <mml:mi>Z</mml:mi> <mml:mi>γ</mml:mi> <mml:mi>γ</mml:mi> </mml:mrow> </mml:math> system. The data are compared with predictions from MC event generators which are accurate to up to next-to-leading order in QCD. The cross-section measurements are used to set limits on the coupling strengths of dimension-8 operators in the framework of an effective field theory.