  We derive the formula of energy flux for the hierarchical equations of motion
(HEOM) method with the help of stochastic decoupling technique. The resulting
expression is a combination of the terms in the first two layers of the
hierarchy. The formula is applied to a three-level "heat engine" coupled
respectively to three baths, of which two for heat sources and one for work
dump. All the three baths are modeled by collections of harmonic oscillators.
We illustrate the proper parameterizing to converge the third bath to the
"work-dump" limit. As an example, the effect of the engine parameters on
working efficiency is studied.
