The exact nuclear time-dependent potential energy surface arises from the exact decomposition of the electronic and nuclear motions, recently presented in [A. Abedi, N.T. Maitra, and E.K.U. Gross, Phys. Rev. Lett. 105, 123002 (2010)]. Such time-dependent potential drives nuclear motion and fully accounts for the coupling to the electronic subsystem. We investigate the features of the potential in the context of electronic non-adiabatic processes and employ it to study the performance of the classical approximation on nuclear dynamics. We observe that the potential, after the nuclear wave packet splits at an avoided crossing, develops dynamical steps connecting different regions, along the nuclear coordinate, in which it has the same slope as one or the other adiabatic surface. A detailed analysis of these steps is presented for systems with different non-adiabatic coupling strength. The exact factorisation of the electron-nuclear wave function is at the basis of the decomposition. In particular, the nuclear part is the true nuclear wave function, solution of a time-dependent Schrödinger equation and leading to the exact many-body density and current density. As a consequence, the Ehrenfest theorem can be extended to the nuclear subsystem and Hamiltonian, as discussed here with an analytical derivation and numerical results.