Heat operator approach to quantum stochastic thermodynamics in the strong-coupling regime

Heat exchanged between an open quantum system and its environment exhibits fluctuations that carry crucial signatures of the underlying dynamics. Within the well-established two-point measurement scheme, we identify a “heat operator,” whose moments with respect to the vacuum state of a thermofield-doubled Hilbert space correspond to the stochastic moments of the heat exchanged with a bath. This recasts heat statistics as a unitary time evolution problem, which we solve by combining chain-mapped reservoirs with tensor-network propagation. In a multibath setup, all total and bath-resolved heat moments then follow from a single pure state evolution. We employ this approach to compute transient and steady-state heat fluctuations in Ohmic spin-boson models in and out of equilibrium, accessing the challenging low-temperature and long memory time regimes of the environment. In the nonequilibrium case, we show a crossover in the Fano factor from super-Poissonian to nearly Poissonian statistics under strong coupling asymmetry, corresponding to thermal rectification behavior. The method applies to noninteracting (bosonic or fermionic) nonequilibrium environments with arbitrary spectral densities, offering a powerful, nonperturbative framework for understanding heat transfer in open quantum systems.