We consider exponential Lawson multistep methods for the time integration of the equations of motion associated with the multi-configuration time-dependent Hartree-Fock (MCTDHF) approximation for high-dimensional quantum dynamics. These provide high-order approximations at a minimum of evaluations of the computationally expensive nonlocal potential terms, and have been found to enable stable long-time integration. In this work, we prove convergence of the numerical approximation on finite time intervals under minimal regularity assumptions on the exact solution. A numerical illustration shows adaptive time propagation based on our methods.Mathematics Subject Classification. 65L05, 65L06, 65M12, 65M15, 65M20, 81-08. on the Hilbert space L 2 . Here, A : D ⊆ L 2 → L 2 is a self-adjoint differential operator and B a generally unbounded nonlinear operator. We will denote the flow of −iĤ by E −iĤ , that is, ψ(t) = E −iĤ (t, ψ 0 ) is the solution of ∂ t ψ = −iĤ(ψ), ψ(0) = ψ 0 , and similarly for E A and E B .Equations of this type commonly arise from model reductions of high-dimensional quantum dynamical systems serving to make many-particle quantum simulations computationally tractable. Examples of such model reductions are the Gross-Pitaevskii equation (GPE) for Bose-Einstein condensates (BEC) [41], which constitutes a meanfield theory for ultracold dilute bosonic gases, the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method for ultracold bosonic atoms [1], the multiconfigurational time-dependent Hartree-Fock (MCTDHF) method [27,32,48] and its most current extension, the time-dependent complete-active-space self-consistent-field (TD-CASSCF) method for electron dynamics in atoms and small molecules [44], and timedependent density functional theory (TDDFT) for extended systems such as large molecules, nanostructures