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 multi-bath 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, non-perturbative framework for understanding heat transfer in open quantum systems.