  We present two diagrammatic Monte Carlo methods for quantum systems coupled
with harmonic baths, whose dynamics are described by integro-differential
equations. The first approach can be considered as a reformulation of Dyson
series, and the second one, called "bold-thin-bold diagrammatic Monte Carlo",
is based on resummation of the diagrams in the Dyson series to accelerate its
convergence. The underlying mechanism of the governing equations associated
with the two methods lies in the recurrence relation of the path integrals,
which is the most costly part in the numerical methods. The proposed algorithms
give an extension to the work ["Fast algorithms of bath calculations in
simulations of quantum system-bath dynamics", Computer Physics Communications,
to appear], where the algorithms are designed based on reusing the previous
calculations of bath influence functionals. Compared with the algorithms
therein, our methods further include the reuse of system associated functionals
and show better performance in terms of computational efficiency and memory
cost. We demonstrate the two methods in the framework of spin-boson model, and
numerical experiments are carried out to verify the validity of the methods.
