  We introduce a method for calculating the stationary state of a translation
invariant array of weakly coupled cavities in the presence of dissipation and
coherent as well as incoherent drives. Instead of computing the full density
matrix our method directly calculates the correlation functions which are
relevant for obtaining all local quantities of interest. It considers an
expansion of the correlation functions and their equations of motion in powers
of the photon tunneling rate between adjacent cavities, leading to an exact
second order solution for any number of cavities. Our method provides a
controllable approximation for weak tunneling rates applicable to the strongly
correlated regime that is dominated by nonlinearities in the cavities and thus
of high interest.
