  We investigate the full-counting statistics (FCS) of energy transport carried
by electrons in molecular junctions for the Anderson-Holstein model in the
polaronic regime. Using two-time quantum measurement scheme, generating
function (GF) for the energy transport is derived and expressed as a Fredholm
determinant in terms of Keldysh nonequilibrium Green's function in the time
domain. Dressed tunneling approximation is used in decoupling the phonon cloud
operator in the polaronic regime. This formalism enables us to analyze the time
evolution of energy transport dynamics after a sudden switch-on of the coupling
between the dot and the leads towards the stationary state. The steady state
energy current cumulant GF in the long time limit is obtained in the energy
domain as well. Universal relations for steady state energy current FCS are
derived under finite temperature gradient with zero bias and this enables us to
express the equilibrium energy current cumulant by a linear combination of
lower order cumulants. Behaviors of energy current cumulants in steady state
under temperature gradient and external bias are numerically studied and
explained. Transient dynamics of energy current cumulants is numerically
calculated and analyzed. The universal scaling of normalized transient energy
cumulants is found under both temperature gradient and external bias.
