  We present a new analytic expression for the heat current through a general
harmonic network coupled with Ohmic reservoirs. We use a new method that
enables us to express the stationary state of the network in terms of the
eigenvectors and eigenvalues of a generalized cubic eigenvalue problem. In this
way, we obtain exact formulae for the heat current and the local temperature
inside the network. Our method does not rely in the usual assumptions of weak
coupling to the environments or on the existence of an infinite cutoff in the
environmental spectral densities. We use this method to study non-equilibrium
processes without the weak coupling and Markovian approximations. As a first
application of our method, we revisit the problem of heat conduction in 2D and
3D crystals with binary mass disorder. We complement previous results showing
that for small systems the scaling of the heat current with the system size
greatly depends on the strength of the interaction between system and
reservoirs. This somewhat counter-intuitive result seems not to have been
noticed before.
