Light-front wave functions play a central role in the program of understanding the structure of hadrons as QCD bound states. Using continuum Schwinger methods, based on Dyson–Schwinger and Bethe–Salpeter equations, they can be computed directly within a framework connected to QCD. For light pseudoscalar mesons, previous studies revealed an approximate separability of longitudinal and transverse lightcone momentum dependences in the LFWFs, leading to a simple relation between distribution functions and amplitudes. In this work, we extend those previous studies to the case of pseudoscalar charmonium and bottomonium, using the fictitious _s meson as a benchmark. Motivated by the observed deviations, we propose a modified non-separable LFWF ansatz that successfully reproduces the properties of heavy pseudoscalar quarkonia and allows the calculation of zero-skewness generalised parton distribution functions, electromagnetic and gravitational form factors, and transverse charge and mass distributions.