  We present an application of the Extended Stochastic Liouville-von Neumann
equations (ESLN) method introduced earlier [PRB 95, 125124 (2017); PRB 97,
224310 (2018)] which describes the dynamics of an exactly thermalised open
quantum system reduced density matrix coupled to a non-Markovian harmonic
environment. Critically, the combined system of the open system fully coupled
to its environment is thermalised at finite temperature using an imaginary time
evolution procedure before the application of real time evolution. This
initialises the combined system in the correct canonical equilibrium state
rather than being initially decoupled. We apply our theory to the spin-boson
Hamiltonian and develop a number of competing ESLN variants designed to reduce
the numerical divergence of the trace of the open system density matrix. We
find that a careful choice of the driving noises is essential for improving
numerical stability. We also investigate the effect of applying higher order
numerical schemes for solving stochastic differential equations, such as the
Stratonovich-Heun scheme, and conclude that stochastic sampling dominates
convergence with the improvement associated with the numerical scheme being
less important for short times but required for late times. To verify the
method and its numerical implementation, we consider evolution under a fixed
Hamiltonian and show that the system either remains in, or approaches, the
correct canonical equilibrium state at long times. Additionally, evolution of
the open system under non-equilibrium Landau-Zener (LZ) driving is considered
and the asymptotic convergence to the LZ limit was observed for vanishing
system-environment coupling and temperature. When coupling and temperature are
non-zero, initially thermalising the combined system at a finite time in the
past was found to be a better approximation of the true LZ initial state than a
pure state.
