Note that in the original model of Ref. [ 62], I 1 was a function of r . However, in our calculations, we used the constant value of I 1 at r 1 , i.e., at the equilibrium OH distance of the electronic ground state. In the above expressions, α is the COH angle, r CC is the CC distance and r CH is the CH distance at the equilibrium ground state geometry. All the parameters of the model are listed in Table ; the reported values of α, r CC and r CH have been determined by optimizing the phenol molecule using density functional theory at the B3LYP/6-311G(d,p) level of theory with Gaussian16 75 . Diagonalization of the electronic Hamiltonian (A1) yields the adiabatic PESs, shown in Fig. . The PESs and corresponding adiabatic states are labelled S 0 (blue surface), S 1 (orange surface) and S 2 (red surface), and present two conical intersections that are clearly identified in Fig. a.t.u. (about 1.2 fs). In this Lanczos procedure, the Krylov space is expanded until the norm of a new vector is less than 10 -10 . The wavepacket is represented in the diabatic representation, and in each electronic state the wavepacket is expanded in a two-dimensional direct product grid or basis set. The primitive basis sets associated with the angle θ and with the r distance are a real Fourier series (64 basis functions) and a particle-in-a-box in the range [0.5, 60] bohr with 2048 basis functions, respectively. This very large range along r is needed to avoid reflection before the end of the propagation time. Since the quantum dynamics is performed in the diabatic representation, the corresponding initial wavepacket is slightly different from the adiabatic one used for the trajectories. In the diabatic representation, the initial wavepacket is given by Eq. (B1) in the third diabatic state, which corresponds to a population in S 1 of more than 97% instead of the 100% for the trajectories.