A search for the decay of the Higgs boson to a $Z$ boson and a light, pseudoscalar particle, $a$, decaying respectively to two leptons and to two photons is reported. The search uses the full LHC Run 2 proton-proton collision data at $\sqrt{s}=13$ TeV, corresponding to 139 fb$^{-1}$ collected by the ATLAS detector. This is one of the first searches for this specific decay mode of the Higgs boson, and it probes unexplored parameter space in models with axion-like particles (ALPs) and extended scalar sectors. The mass of the $a$ particle is assumed to be in the range 0.1-33 GeV. The data are analysed in two categories: a merged category where the photons from the $a$ decay are reconstructed in the ATLAS calorimeter as a single cluster, and a resolved category in which two separate photons are detected. The main background processes are from Standard Model $Z$ boson production in association with photons or jets. The data are in agreement with the background predictions, and upper limits on the branching ratio of the Higgs boson decay to $Za$ times the branching ratio $a\to\gamma\gamma$ are derived at the 95% confidence level and they range from 0.08% to 2% depending on the mass of the $a$ particle. The results are also interpreted in the context of ALP models.