  We develop and test a computational framework to study heat exchange in
interacting, nonequilibrium open quantum systems. Our iterative full counting
statistics path integral (iFCSPI) approach extends a previously
well-established influence functional path integral method, by going beyond
reduced system dynamics to provide the cumulant generating function of heat
exchange. The method is straightforward; we implement it for the nonequilibrium
spin boson model to calculate transient and long-time observables, focusing on
the steady-state heat current flowing through the system under a temperature
difference. Results are compared to perturbative treatments and demonstrate
good agreement in the appropriate limits. The challenge of converging
nonequilibrium quantities, currents and high order cumulants, is discussed in
detail. The iFCSPI, a numerically exact technique, naturally captures strong
system-bath coupling and non-Markovian effects of the environment. As such, it
is a promising tool for probing fundamental questions in quantum transport and
quantum thermodynamics.
