  We present a wavefunction methodology to account for finite temperature
initial conditions in the quantum Rabi model. The approach is based on the
Davydov-Ansatz together with a statistical sampling of the canonical harmonic
oscillator initial density matrix. Equations of motion are gained from a
variational principle and numerical results are compared to those of the
thermal Hamiltonian approach. For a system consisting of a single spin and a
single oscillator and for moderate coupling strength, we compare our new
results with full quantum ones as well as with other Davydov-type results based
on alternative sampling/summation strategies. All of these perform better than
the ones based on the thermal Hamiltonian approach. The best agreement is shown
by a Boltzmann weighting of individual eigenstate propagations. Extending this
to a bath of many oscillators will, however, be very demanding numerically. The
use of any one of the investigated stochastic sampling approaches will then be
favorable.
