  We describe a numerical scheme for exactly simulating the heat current
behavior in a quantum harmonic chain with self-consistent reservoirs.
Numerically-exact results are compared to classical simulations and to the
quantum behavior under the linear response approximation. In the classical
limit or for small temperature biases our results coincide with previous
calculations. At large bias and for low temperatures the quantum dynamics of
the system fundamentally differs from the close-to-equilibrium behavior,
revealing in particular the effect of thermal rectification for asymmetric
chains. Since this effect is absent in the classical analog of our model, we
conclude that in the quantum model studied here thermal rectification is a
purely quantum phenomenon, rooted in the quantum statistics.
