I've only looked cross reading on the above cited paper by Verstraelen et al. The key point is that you deal with an explicitly time-dependent Hamiltonian and thus energy is not conserved here. That's also the case in classical mechanics: If you have forces derived from a time-dependent potential, you can still use the Hamilton principle but total energy is not conserved.
This is easily seen in the quantum case using the fundamental rule that the time-derivative of an observable (represented by a self-adjoint operator ##\hat{O}##) is represented by the operator
$$\mathring{\hat{O}}=\frac{1}{\mathrm{i} \hbar} [\hat{O},\hat{H}]+\partial_t \hat{O},$$
where the partial time derivative refers to the explicit time-dependence of ##\hat{O}## only (i.e., the only time-dependence at all in the Schrödinger picture of time evolution). Now set ##\hat{O}=\hat{H}##. Then you have
$$\mathring{\hat{H}}=\partial_t \hat{H},$$
i.e., energy is conserved if and only if the Hamiltonian is not explicitly time-dependent.
Of course the time evolution is still unitary. In the Schrödinger picture the states evolve in time with the unitary operator ##\hat{U}(t,t_0)## definied by the "operator Schrödinger equation",
$$\mathrm{i} \hbar \partial_t \hat{U}(t,t_0)=\hat{H}(t) \hat{U}(t,t_0).$$
This time-evolution operator now is a function of ##t## and ##t_0## separately, not only on ##(t-t_0)## as in the case of a time-independent Hamiltonian. The formal solution of the operator Schrödinger equation is given, using the time-ordering symbol, by
$$\hat{U}(t,t')=\mathcal{T}_c \exp \left [-\frac{\mathrm{i}}{\hbar} \int_{t_0}^t \mathrm{d} t' \hat{H}(t') \right].$$
The timeordering symbol tells you that in the expansion of the operator exponential you have to time order the Hamiltonians in the occurring multiple integrals always such the the time arguments are ordered as increasing from right to left.
Now the paper deals with the interesting situation that you of course also do not have the usual equilibrium situation since the canonical (or grand canonical) stat. op. ##\exp[-\hat{H}/(k_B T)]## is not stationary since ##\hat{H}## depends explicitly ontime, but that you still can define a generalized "Gibbs ensemble" by using the maximum-entropy method under appropriate constraints.
