We report our investigation on the doubly virtual TFFs $F_{{rm P}gamma^*}(Q^2_1,Q^2_2)$ for the ${rm P}togamma^*(q_1)gamma^*(q_2) ;({rm P}=pi^0,eta,eta)$ transitions using the light-front quark model (LFQM). Performing a LF calculation in the exactly solvable manifestly covariant Bethe-Salpeter (BS) model as the first illustration, we used $q^+_1=0$ frame and found that both LF and manifestly covariant calculations produce exactly the same results for $F_{{rm P}gamma^*}(Q^2_1,Q^2_2)$. This confirms the absence of the LF zero mode in the doubly virtual TFFs. We then mapped this covariant BS model to the standard LFQM using the more phenomenologically accessible Gaussian wave function provided by the LFQM analysis of meson mass spectra. For the numerical analyses of $F_{{rm P}gamma^*}(Q^2_1,Q^2_2)$, we compared our LFQM results with the available experimental data and the perturbative QCD (pQCD) and the vector meson dominance (VMD) model predictions. As $(Q^2_1, Q^2_2)toinfty$, our LFQM result for doubly virtual TFF is consistent with the pQCD prediction, i.e. $F_{{rm P}gamma^*}(Q^2_1, Q^2_2)sim 1/(Q^2_1 + Q^2_2)$, while it differs far from the result of VMD model which behaves $F^{rm VMD}_{{rm P}gamma^*}(Q^2_1, Q^2_2)sim 1/(Q^2_1 Q^2_2)$. Our LFQM prediction for $F_{etagamma^*}(Q^2_1,Q^2_2)$ shows an agreement with the very recent experimental data obtained from the BaBar collaboration for the ranges of $2< Q^2_1, Q^2_1 <60$ GeV$^2$.