  The standard approach for path integral Monte Carlo simulations of open
quantum systems is extended as an efficient tool to monitor the time evolution
of coherences (off-diagonal elements of the reduced density matrix) also for
strong coupling to environments. Specific simulations are performed for two
level systems embedded in Ohmic and sub-Ohmic reservoirs in the domains of
coherent and incoherent dynamics of the polarization. In the latter regime, the
notorious difficulty to access the long time regime is overcome by combining
simulations on moderate time scales with iteratively calculated initial
densities. This allows to extract relaxation rates for sub-Ohmic environments
at finite temperatures and over a broad range of couplings and to compare them
to analytical predictions. The time evolution of the von Neumann entropy
provides insight into the quantum phase transition at thermal equilibrium from
a delocalized to a localized state at zero temperature.
