7 in an infinite 3D system there would be no outside for photons to escape to, and they would instead mix with the atomic excitations to form two polariton modes [28]. (In the infinite system limit there are also retardation effects neglected by our master equation, but for prac- tical experimental sizes the boundary is reached before these become significant.) In the following, we therefore restrict ourselves to one-dimensional lattices. dipole interaction becomes the dispersion c ( p ) of collec- tive electronic excitations, which allows this term to be treated non-perturbatively. The first and second order terms d$_{α}$ ( p ), e$_{αβ}$ ( p ), which we do still have to treat per- turbatively, describe interaction processes where an elec- tronic excitation emits or absorbs one or two phonons. Because of the discrete translational symmetry, wavevec- tor (quasimomentum) is conserved mod the Brillouin zone, with the laser contributing its wavevector k$_{L}$ to each excitation it creates. In the large system limit N → ∞ , these Fourier sums converge in one dimension (except for a logarithmic di- vergence at the points p = ± k$_{0}$ ) but diverge in two or more dimensions. While they can still be calculated for finite sizes, in the presence of such strong finite size ef- fects we cannot neglect the difference between periodic and open boundary conditions. Physically, this is be- cause in 2D and 3D photons have a significant probabil- ity of being reabsorbed after travelling a large distance; Similarly, the reset operator R ( ρ ) in Eq. (20) can now be expanded as with C. The reset operator R ( ρ ) = R $^{(0)}$( ρ ) + R $^{(1)}$( ρ ) + R $^{(2)}$( ρ ) , (26) R $^{(0)}$( ρ ) = ∑ p Re c ( p ) s$_{p}$ρs † $_{p}$, R $^{(1)}$( ρ ) = ∑ p , q ,α i √ N $^{[}$Im d$_{α}$ ( q ) s$_{q}$ρs † p + $_{q}$( b$_{p}$$_{α}$ + b † − p $_{α}$) − Im d$_{α}$ ( p + q ) ( b$_{p}$$_{α}$ + b † − p $_{α}$) s$_{q}$ρs † p + q ] , R $^{(2)}$( ρ ) = ∑ p , q , r ,α,β 1 N $^{[}$Re e$_{αβ}$ ( r + q )( b$_{p}$$_{α}$ + b † − p $_{α}$)( b$_{q}$$_{β}$ + b † − q $_{β}$) s$_{r}$$_{−}$$_{p}$ρs † r + q − 2Re e$_{αβ}$ ( r )( b$_{p}$$_{α}$ + b † − p $_{α}$) s$_{r}$$_{−}$$_{p}$ρs † r + $_{q}$( b$_{q}$$_{β}$ + b † − q $_{β}$) + Re e$_{αβ}$ ( r − p ) s$_{r}$$_{−}$$_{p}$ρs † r + $_{q}$( b$_{p}$$_{α}$ + b † − p $_{α}$)( b$_{q}$$_{β}$ + b † − q $_{β}$) ] . (27) Because e ipx + e − ipx = 2 cos px is real and e ipx − e − ipx = 2i sin px is pure imaginary, even functions such as C$_{ij}$, E αβ ij transform real part to real part (Re C$_{ij}$ → Re c ( p )) while odd functions such as Re D α ij transform real part to imaginary part (Re D α ij → i Im d$_{α}$ ( p )). The zeroth order term describes the spontaneous decay of the electronic excitations, with a wavevector-dependent rate because of dipole-dipole interaction. The first and sec- ond order terms are decays that also create or destroy phonons. In this subsection, we describe a modified form of Feyn- man diagrams that can be used on open systems de- scribed by a time independent Markovian master equa- tion. Our diagrams are similar to those of Keldysh [29, 30], but our starting point is a master equation rather than a system plus bath Hamiltonian. Our diagrams divide the Liouvillean superoperator L ( ρ ) = ˙ ρ into a propagation part and a vertex part, and treat the vertex part as a perturbation, in the same way as conventional Feynman diagrams divide the Hamil- tonian operator [31]. Propagation terms are those in H$_{cond}$ with a matching creation and annihilation oper- ator, i.e. − i 2 $_{ℏ}$c ( p ) s † $_{p}$s$_{p}$ for an electronic excitation and $_{ℏ}$ν$_{α}$ b † p $_{α}$b$_{p}$$_{α}$ for a phonon. (These lines represent propaga- tion through time as well as space, so it makes sense for a non-propagating particle such as our phonon to have a propagation line.) The remaining terms in H$_{cond}$ and all terms in R ( ρ ) are vertices. Fig. 4(a) shows the vertices appearing in our master equation. Propagation lines are labelled with whatever indices Under the master equation (11), the density matrix converges to a steady state, and for cooling we would like it to converge within a reasonable time to a state with the lowest possible temperature, or equivalently phonon number. In this section, we calculate this steady state using an open system version of Feynman diagrams. IV. STEADY STATE DENSITY MATRIX A. Open system Feynman diagrams