which depends explicitly on the nuclear wavefunction, on the electronic wavefunction via the presence of the TDVP, and acts on the parametric dependence of Φ(r, t; R) as a spatial derivative. The ambiguity of the product form of the molecular wavefunction in Eq. ( ) is partially eliminated by imposing the normalization condition |Φ(r, t; R)| 2 dr = 1 ∀R, t, which allows to identify |Φ(r, t; R)| 2 as a conditional probability density and |χ(R, t)| 2 = |Ψ(r, R, t)| 2 dr as a marginal probability density that yields the nuclear density from the full wavefunction. Equation ( ) is thus unique up to a gauge encoded in a phase factor e (i/ )θ(R,t) (with θ(R, t) a real function): multiplying the nuclear wavefunction by this factor and the electronic term by its complex conjugate, Eq. (3) remains unaffected and the time-dependent potentials transform as standard gauge potentials. Therefore, this ambiguity has to be eliminated by imposing a choice of gauge as discussed in Section II A.