  The time-dependent numerical renormalization group (td-NRG) [Anders et al.
Phys. Rev. Lett. {\bf 95}, 196801 (2006)] offers the prospect of investigating
in a non-perturbative manner the time-dependence of local observables of
interacting quantum impurity models at all time scales following a quantum
quench. We present a generalization of this method to arbitrary finite
temperature by making use of the full density matrix [Weichselbaum et al. Phys.
Rev. Lett. {\bf 99}, 076402 (2007)]. We show that all terms in the projected
density matrix $\rho^{i\to f} = \rho^{++} + \rho^{--} + \rho^{+-}+\rho^{-+}$ in
the time-evolution of a local observable may be evaluated in closed form, with
$\rho^{+-}=\rho^{-+}=0$. The expression for $\rho^{--}$ is shown to be finite
at finite temperature, becoming negligible only in the limit of vanishing
temperatures. We prove that this approach recovers the short-time limit for the
expectation value of a local observable exactly at arbitrary temperatures. In
contrast, the corresponding long-time limit is only recovered exactly in the
limit $\Lambda\rightarrow 1^{+}$. This limit is impractical within NRG. We
suggest how to overcome this problem by noting that the long-time limit becomes
increasingly more accurate on reducing the size of the quantum quench. This
suggests an improved generalized td-NRG approach in which the system is
time-evolved between the initial and final states via a sequence of small
quantum quenches within a finite time interval instead of a single large
quantum quench. The formalism for this is provided, generalizing the td-NRG to
multiple quantum quenches, periodic switching, and general pulses. This
formalism, like our finite temperature generalization of the single-quench
case, rests only on the NRG approximation. The results are illustrated by
application to the Anderson impurity model.
