ABSTRACT The galaxy size–stellar mass and central surface density–stellar mass relationships are fundamental observational constraints on galaxy formation models. However, inferring the physical size of a galaxy from observed stellar emission is non-trivial due to various observational effects, such as the mass-to-light ratio variations that can be caused by non-uniform stellar ages, metallicities, and dust attenuation. Consequently, forward-modelling light-based sizes from simulations is desirable. In this work, we use the skirt dust radiative transfer code to generate synthetic observations of massive galaxies ($M_{*}∼ 10^{11}\, \rm {M_{⊙ }}$ at z = 2, hosted by haloes of mass $M_{\rm {halo}}∼ 10^{12.5}\, \rm {M_{⊙ }}$) from high-resolution cosmological zoom-in simulations that form part of the Feedback In Realistic Environments project. The simulations used in this paper include explicit stellar feedback but no active galactic nucleus (AGN) feedback. From each mock observation, we infer the effective radius (Re), as well as the stellar mass surface density within this radius and within $1\, \rm {kpc}$ (Σe and Σ1, respectively). We first investigate how well the intrinsic half-mass radius and stellar mass surface density can be inferred from observables. The majority of predicted sizes and surface densities are within a factor of 2 of the intrinsic values. We then compare our predictions to the observed size–mass relationship and the Σ1−M⋆ and Σe−M⋆ relationships. At z ≳ 2, the simulated massive galaxies are in general agreement with observational scaling relations. At z ≲ 2, they evolve to become too compact but still star forming, in the stellar mass and redshift regime where many of them should be quenched. Our results suggest that some additional source of feedback, such as AGN-driven outflows, is necessary in order to decrease the central densities of the simulated massive galaxies to bring them into agreement with observations at z ≲ 2.