  We report on the investigation of height distributions (HDs) and spatial
covariances of two-dimensional surfaces obtained from extensive numerical
simulations of the celebrated Clarke-Vvedensky (CV) model for homoepitaxial
thin film growth. In this model, the effect of temperature, deposition flux,
and strengths of atom-atom interactions are encoded in two parameters: the
diffusion to deposition ratio $R=D/F$ and $\varepsilon$, which is related to
the probability of an adatom "breaking" a lateral bond. We demonstrate that the
HDs present a strong dependence on both $R$ and $\varepsilon$, and even after
the deposition of $10^5$ monolayers (MLs) they are still far from the
asymptotics in some cases. For instance, the temporal evolution of the HDs'
skewness (kurtosis) displays a pronounced minimum (maximum), for small $R$ and
$\varepsilon$, and only at long times it passes to increase (decrease) toward
its asymptotic value. However, it is hard to determine whether they converge to
a single value or different nonuniversal ones. For large $R$ and/or
$\varepsilon$, on the other hand, these quantities clearly converge to the
values expected for the Villain-Lai-Das Sarma (VLDS) universality class. A
similar behavior is observed in the spatial covariances, but with weaker
finite-time effects, so that rescaled curves of them collapse quite well with
the one for the VLDS class at long times. Simulations of a model with limited
mobility of particles, which captures some essential features of the CV model
in the limit of irreversible aggregation ($\varepsilon=0$), reveal a similar
scenario. Overall, these results point out that the study of fluctuations in
homoepitaxial thin films' surfaces can be a very difficult task and shall be
performed very carefully, once typical experimental films have $\lesssim 10^4$
MLs, so that their HDs and covariances can be in the realm of transient
regimes.
