{'fraction_non_alphanumeric': 0.031161966717974037, 'fraction_numerical': 0.008393111259547528, 'mean_word_length': 4.799898938858009, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Classical field model for arrays of photon condensates

We introduce a classical phasor model for the description of multimode photon condensates that thermalize through repeated absorptions and reemissions by dye molecules. Thermal equilibrium is expressed through the fluctuation-dissipation relation that connects the energy damping to spontaneous emission fluctuations. We apply our model to a photonic Josephson junction (two coupled wells) and to one- and two-dimensional arrays of photon condensates. In the limit of zero pumping and cavity losses, we recover the thermal equilibrium result, but in the weakly driven-dissipative case in the canonical regime, we find suppressed density and phase fluctuations with respect to the ideal Bose gas.


I. INTRODUCTION
Bose-Einstein condensates (BEC) of photons can be realized by embedding dye molecules in a high-quality optical cavity [1]. When the rate of absorption and reemission of photons by the dye molecules is much higher than the photon losses, the photons are brought to thermal equilibrium with the molecular rovibrational states [2][3][4][5]. These are in turn thermalized by collisions with solvent molecules. The photon gas then assumes the solvent temperature and features Bose-Einstein condensation above the saturation density [1].
Since the ideal Bose gas is one of the simplest systems from the theoretical point of view [6], one may wonder whether this system presents any theoretical challenges or is just a nice platform to perform some demonstration experiments of elementary textbook physics. When the photon losses are fully absent, the latter is the case [7,8], since then the system is guaranteed to relax to the thermal state. When on the other hand losses cannot be fully neglected, as is the case in current experiments [1,[9][10][11][12], the physics becomes richer and our understanding of the interplay between pumping, losses, external potentials and thermalization is still not complete. Most experiments are performed with harmonically trapped photons, but more recently there has been experimental progress in the creation of double well and periodic potentials [12].
The simplest theoretical description of out of equilibrium photon condensation consists of rate equations for the occupation of the single particle energy levels [13]. The ingredient that is missing here is the coherence between the single particle states, which is necessary to form localized photon wave packets [14]. Such a state is for example formed when the photon condensate is pumped with a finite size pumping spot [15]. An extension of the rate equation model to take the spatial distribution of the molecules into account was developed by Hesten et al. [16].
On the other size of the complexity spectrum are the quantum optics based approaches where the master equation for the open quantum system [17] consisting of the coupled molecular and photonic system is solved [14,18,19]. Apart from the needed computational re-sources, a disadvantage of this master equation approach is that it is hard to include correlation between photons and molecules. It has been shown that the correlation between photons and molecules can affect the density fluctuations [20,21]. The underlying physics is elementary: when a large number of molecules is present, they form a bath for particle exchange as in the grandcanonical statistical ensemble. In this so-called grandcanonical regime, the number fluctuations are large, due to the unrestricted exchange of particles between system and bath. When on the other hand the number of molecules is small, it becomes unlikely that all photons are simultaneously absorbed; hence number fluctuations are reduced. In order to describe this physics correctly, it is essential to include the correlations between the number of photons and the number of excited molecules. From the quantum optics side, these classical molecule-photon correlations can be included most easily in a quantum trajectory approach. For single mode photon condensates, this was implemented in Ref. [22]. It was also shown in this work that the quantum optical description is well reproduced by a classical field model.
For the description of close to ideal Bose gases, classical field theory has proved to be an indispensable tool [23]. Most experiments with weakly interacting ultracold Bose gases are excellently described by the Gross-Pitaevskii equation (GPE) [24]. For exciton-polariton condensates, a system closely related to photon condensates, many experiments are modeled with a generalized Gross-Pitaevskii equation that includes pumping and losses [25]. Exciton-polaritons are hybrid light-matter quasiparticles that interact with each other thanks to their excitonic component, setting them apart from noninteracting photons. In practice however, the interaction energy in experimentally realized polariton condensates is quite small due the relatively small value of the interaction constant [26]. In parallel to photon condensates, cavity losses make it necessary to pump the system in order to reach a steady state. When cavity losses are small, thermal equilibrium is closely approached [27].
Classical field theories exist at various levels of complexity. The most elementary version is the standard GPE, containing only kinetic and interaction energy, that is suitable for the description of zero temperature weakly interacting bosons [24]. Particle exchange with a reservoir can be added by including an imaginary term [25]. Energy exchange between the bosons and their environment can be modeled by making the prefactor of the time derivative complex. It was originally introduced to model the friction between the superfluid and normal components of liquid Helium [28], but has also been employed in the description of ultracold atoms [23,29] and polariton condensates [30]. The standard Gross-Pitaevskii equation assumes perfect coherence of the bosons. Decoherence can be incorporated by including some stochasticity. This can for example be derived in the truncated Wigner approximation [23,[31][32][33].
For photon condensates, a classical field description has been used to model their phase coherence [34]. The 'phasor model' version of the classical field description was developed in the context of laser physics [35,36] and compares favorably with a quantum trajectory description for single-mode photon condensates [22]. It is the purpose of this paper to extend the phasor model to multimode photon condensates.
We show in Sec. II that the thermalization of the photon gas by the molecules is described by adding the same term that models the friction between superfluid and normal components in atomic condensates. The fluctuations in the phasor model are shown to be related to this friction through a fluctuation-dissipation relation. In Sec. III, we recapitulate the physics of single mode photon condensates, with specific attention to the regimes of small and large density fluctuations, the 'canonical' and 'grandcanonical' regimes respectively. We then analyze in detail the case of two coupled photon traps in Sec. IV, a photonic Josephson junction (PJJ). For this system, we analytically compute the density and phase fluctuations in the linearized Bogoliubov approximation. In the limit of zero losses, we recover the equilibrium correlators in the classical regime, justifying our model as an adequate description of photon condensates. Finally, in Sec. V, we extend our analysis to 1D and 2D lattices of photon condensates. We find that the spatial coherence of photon condensates is much better in the canonical regime as compared to the grandcanonical regime.

A. Kennard-Stepanov relation and energy relaxation
For dye molecules that interact sufficiently strongly with their solvent, the emission and absorption coefficients are related by the Kennard-Stepanov (KS) law, which reads at the inverse temperature β = 1/(k B T ) [2][3][4] Here B 12 (B 21 ) is the Einstein coefficient for absorption (emission) of a photon, ω 0 is the molecular transition frequency and ω is the photon frequency. Here we absorbed possible degeneracy factors in a renormalization of the molecular transition frequency. A sketch of an absorption-emission spectrum that satisfies the KS relation is shown in Fig. 1.
The KS law brings the photon gas to thermal equilibrium, as can be seen from the steady state of the kinetic equation for the photon number in a single mode cavity: Setting dn/dt = 0, one finds n = (B 12 /B 21 − 1) −1 , which reduces to the Bose-Einstein distribution when the KS relation (1) is used. For a single mode of noninteracting photons, the KS relation can be straightforwardly implemented in a theoretical model, but for multimode systems, the photonic frequency ω is not a priori known. When the photonic frequency is still close to a certain cavity frequency ω c , one can proceed by writing the KS relation as where ∆ = ω c −ω 0 is the cavity-molecule detuning. To be specific, we will assume that all the energy dependence is in the absorption coefficient. This is valid close to the maximum of the emission. However, we do not expect that our results will be significantly altered when the energy dependence is moved to the emission or distributed between absorption and emission. Choosing ω c as the zero of energy and replacing ω → i∂ t , one obtains In a classical field description, the photon dynamics is described by a generalized Gross-Pitaevskii equation (gGPE), setting = 1, Here,T formally represents the kinetic energy and M 1 (2) are the number of ground state (excited) molecules. The cavity loss rate is denoted by γ. The wave function is position dependent (explicitly: ψ = ψ(j), where j labels a lattice cite), as well as the number of ground state and excited molecules, that satisfy at all times M 1 (j) + M 2 (j) = M , where M is the number of dye molecules at each lattice site. Eq. (5) does not contain interaction energy, which is quite negligible in current experiments [37], except for a slow thermo-optical nonlinearity [38], inclusion of which would be a straightforward extension of our model. With Eq. (4), the gGPE becomes where is the energy relaxation rate. For typical photon condensates the relative fluctuations in the number of ground state molecules is small, such that κ can be approximated by a constant. The above derivation made use of the formal substitution ω → i ∂ ∂t in the KS relation, whose validity may be questioned. In order to further justify this approach, we show in Appendix A that the same equation can be rigorously derived for the case where the energy dependent absorption is due to coupling with a lossy bosonic mode.
For κ ≪ 1, which is satisfied if the absorption of photons is slow on the thermal time scale β, the gGPE can be approximated by where iκT ψ forms an imaginary tunneling term, while corrections due to κ in the emission-absorption term were neglected.

B. Fluctuations
A classical field model without fluctuations only captures absorption and stimulated emission. In order to describe spontaneous emissions, fluctuations have to be introduced. One possibility is the heuristic phasor model, where spontaneous emissions are modeled by adding a unit length phasor with random angle to the photonic field [35]: with θ a random phase. This noise should be added at random times, with probability p = dt B 21 M ↑2 (j) in a time interval dt. For the single mode case, this model was demonstrated to show excellent agreement with a full quantum trajectory description [22], motivating its use to describe fluctuations in the multimode regime.

C. Molecule dynamics
The equation of motion (8) has to be coupled to the dynamics of the number of excited molecules. The absorption and stimulated emission dynamics are local and lead to a local change in the number of excited molecules of the form From the energy relaxation term (proportional to κ), we have the contribution where we have neglected relaxation of the molecules in other modes than the cavity mode. Finally, the spontaneous emission (9) is accompanied by a change Note that this change can be positive or negative and is in a given realization not equal to one (this is only true on average). The terminology 'spontaneous emission' may therefore be a bit confusing for this term. The crucial physics that it does capture are the phase diffusion and density fluctuations in the photon condensate [22]. In order to compensate for the excitations, which are lost through the cavity mirrors, and reach a steady state, the system has to be continuously pumped. This is modeled by the following term in the equation of motion for the excited molecules wheren is the targeted steady state number of photons. The Gaussian white term with autocorrelation ξ p (t)ξ p (t ′ ) = δ(t − t ′ ) comes from the shot noise in the excitation of the molecules. We include it here for completeness, but it will turn out that its effect is much smaller than that of the spontaneous emission noise (12).

III. SINGLE MODE PHYSICS: CANONICAL AND GRANDCANONICAL REGIMES
The noninteracting Bose gas in the grandcanonical ensemble has large number fluctuations. When the photons of a photon condensate are coupled to a large number of molecules, the molecules form a reservoir and photon number fluctuations are large. With less molecules in the cavity, the photon number fluctuations are reduced.
The dynamical analysis of the density fluctuations for a single-mode photon condensate as a function of reservoir size and detuning was performed in Ref. [22]. The deviations of the number of photons n and total number of excitations X = M 2 + n from their equilibrium values n andX evolve in linear approximation as Here, the number fluctuation decay rate is given by [39] where the effective reservoir size equals This reduces to the form from Ref. [21] when γ = 0. The Gaussian white noises ξ p,n have zero mean and variance equal to The stochastic term ξ n originates from the fluctuations due to spontaneous emission and can be derived in the diffusion approximation to the phasor model by considering the effect of the spontaneous emission on the variance of the density. In a spontaneous emission ψ = ψ + e iθ , the density variance increases by With the spontaneous emission rate being B 21 M 2 , one arrives at the noise term in Eq. (15). Analogously, the noise term in Eq. (14) originates from the shot noise deviations of the pumping, which is needed to compensate for the photon losses, from its average rate γn. According to Eq. (15), the density fluctuations in the absence of losses are given by [22] Typical experimental photon losses do not significantly alter the density fluctuations [21]. In the limit of a large reservoir (grandcanonical regime) (n 2 ≪ M eff ), one obtains large number fluctuations σ 2 n =n 2 , where in the opposite limit of a small reservoir (canonical regime), one obtains small number fluctuations σ 2 n = M eff ≪n 2 . In the grandcanonical limit, phase jumps occur when the density goes to zero [34], but the overall phase coherence time is still of the Schawlow-Townes form [36,40].

IV. THE PHOTONIC JOSEPHSON JUNCTION
The simplest system to illustrate our model for a lattice of photon condensates is a photonic Josephson juntion (PJJ), which consists of two coupled sites (L and R) [12]. For this example, we will write explicitly the gGPE and relaxation contribution to the molecular dynamics for the left site; the equations for the right site can be obtained by the replacement L ↔ R.

A. Equations of motion
The deterministic part of the equations of motion for the field amplitude on the left hand site reads from Eq.
The relaxation coefficient κ leads to a larger linewidth for the antisymmetric state as compared to the symmetric state. It leads to the following change in the number of excited molecules, cf. Eq. (11), For the remaining equations of motion, we refer to Eqs. (9), (10) and (12).

B. Fluctuations
As in the single mode case, further analytical insight can be obtained by linearizing the equations of motion for small phase and density difference, writing ψ j = n + δn j e iθj . For the phase difference ∆θ = θ R − θ L , one then obtains where the phase diffusion originates from spontaneous emission. The phase diffusion constant is Eq. (22) shows that the energy relaxation parameter κ drives the system to zero relative phase. With Eq. (7), the phase damping and noise are seen to obey the fluctuation-dissipation relation For the photon and excitation number density difference ∆n = δn R − δn L and ∆X = δX R − δX L , one obtains where Γ and σ 2 n are still given by Eqs. (16) and (19). With Ito calculus, equations of motion for the correlation functions can be constructed and the correlators can be evaluated analytically. The full expressions are cumbersome, but both in the limit for large and for small tunneling, the variance of the phase difference takes the simple expression where η ≥ 1 is a noise enhancement factor that describes the increase of fluctuations due to losses: In practice, this factor becomes appreciably larger than one only for large negative detuning. Inserting in Eq. (27) the fluctuation-dissipation relation (24), one finds For the relative density fluctuations, one obtains in the limit of large J In the absence of losses (γ = 0), the results (29) and (30) agree with those for noninteracting bosons in the weak fluctuation regime at thermal equilibrium (see appendix B) . Note that this limiting equilibrium expression is ensured by the fluctuation-dissipation relation, which originates from the Kennard-Stepanov relation. From the density and phase fluctuations, also the first order coherence can be computed in linearized approximation: In the limit of large J one obtains which also reduces to the equilibrium expression when η = 1. The dependence of density fluctuations on the tunneling rate in the presence of losses is shown in Fig. 2(a) for photon condensates both in the grandcanonical (black squares), canonical (green triangles) and intermediate (red dots) regimes. These points were obtained from direct numerical simulations of the equations of motion.
At sufficiently large tunneling rate, the density fluctuations in Bogoliubov approximation (30) are recovered for all photon and molecule numbers. The Bogoliubov approximation breaks down when density fluctuations are too large (at small J), but in equilibrium, the density fluctuations can be computed analytically at any T /J (see appendix B). This expression is shown in Fig. 2(a) with the black dashed line and corresponds very well to the numerical simulations in the deep grandcanonical regime for all tunneling rates.
In the canonical regime, losses affect the density fluctuations significantly. This is seen both in the numerics (red and green symbols) and in the nonequilibrium Bogoliubov (NEB) expression (red and green lines). The density fluctuations in the canonical regime can be understood from the interplay between tunneling and losses. By coupling the two wells, the relative density can fluctuate because of particle exchange. When the tunneling becomes large, the condensate is almost entirely in the symmetric state, such that the relative density fluctuations become small. On the other hand, when the effect of tunneling is much smaller than that of losses, particle exchange is barely possible and relative density fluctuations are suppressed. Consequently, there is a nonmonotoneous dependence of the relative density fluctuations on the tunneling rate. By combining the small and large J expansions of the density fluctuations, the position of the maximal density fluctuations can be estimated to be at the tunneling rate J = 1 2 (T γB 21 n 2 /M 2 ) 1/3 . The first order coherence of the PJJ is shown in Fig.  2(b). At large tunneling, it shows the same behavior as the density fluctuations. Also in analogy with density fluctuations, the coherence for grandcanonical condensates is almost unaffected by experimentally relevant losses. In the canonical regime, on the other hand, we find a again a nonmonotonous behavior. At the tunneling strength where density fluctuations become suppressed, also the first order coherence improves. Upon further decreasing the tunneling strength, the coherence becomes optimal and then worsens again. This worsening at small J is entirely due to phase fluctuations. Indeed, in the limit of zero tunneling, the phases between the two condensates become uncorrelated, such that there is no coherence between the two condensates. From the small J expansion of the coherence, the maximum is found to be at the tunneling rate J = [B 21 nγ/(4 √ 2)] 1/2 . In some numerical simulations (not shown here), we have also observed a long lived antisymmetric state with π phase difference between the two wells. This is not the lowest energy state, but for small tunneling rate the energy relaxation is too weak in order to cause a dynamical instability of the antisymmetric state. The metastability of the antibonding state is relevant for the potential use of photonic Bose-Einstein condensates as analog simulators for (classical) minimization problems [41]: already for the simple two-site problem, there appears the possibility for the system to be stuck in an excited state. A more detailed analysis of metastable states of photon condensates is beyond the scope of this paper and will be deferred to another study.
The previous discussion assumed perfect symmetry of the system. In Fig. 3, we show the sensitivity of the density fluctuations and coherence to the detuning of the wells in the canonical regime. Where the spatial coherence is very good for small energy offset between the wells, it quickly deteriorates when the detuning is in- creased [see Fig. 3(b)]. For larger losses, the coherence is lost for smaller values of the detuning. The reason is that the system goes to a desynchronized state, where the condensates in the two wells have different frequencies and hence no phase coherence. When losses are larger, less particles can be transferred between the two wells and synchronization is lost for smaller detuning. Also the density fluctuations are sensitive to a detuning between the wells [see Fig. 3(a)]. When the system loses synchronisation, the density fluctuations increase, but for larger detuning, the density fluctuations again decrease, the exchange of photons being suppressed when the detuning is much larger than the hopping.

V. LATTICES OF PHOTON CONDENSATES
With the physics of the two-site PJJ understood, we now turn to the study of one-and two-dimensional arrays. From the analytical side, we will again study the dynamics in Bogoliubov approximation, which al-lows to compute the spatial coherence. We complete these calculations with numerical studies of one-and twodimensional systems.
The linear analysis for larger lattices can be performed along the same lines as for the PJJ. The density and phase variables are again introduced as ψ(j) = n + δn(j)e iδθ(j) . The Fourier components of the phase fluctuations are defined by where L is the number of lattice sites, and analogous for δn(j) and δX(j). As observables, we will consider the momentum distribution where ψ k is the Fourier transform of the field ψ(j). It is worth stressing that the momentum distribution N k is different from the Fourier transform of the density δn k . We will also consider the normalized static structure factor In the linear Bogoliubov approximation to Eqs. (8)- (12), the fluctuations obey the equations of motion For a tight binding hamiltonian with hopping amplitude J, the single particle dispersion equals ǫ k = 2J[2 − cos(k x ) − cos(k y )]. The white noise terms have zero average and variance ξ (α) . The diffusion constants are D θ = B 21 M 2 /4n, D n = B 21 M 2n and D X = γn. As in the two-cavity case, the last one appears to play a negligible role. These equations are identical to the ones in the case of the PJJ with the replacement 2J → ǫ k . The momentum distribution can be computed analogously to Eq. (31) as Numerically, we can check the validity of the linear approximation through the static structure factor, shown in Fig. 4(a) and the momentum distribution, shown in Fig.  4(b). In order to stress the similarity with the PJJ, we plot the fluctuations S kx,0 and N kx,0 against the dispersion ǫ k . The parameter values used in the figures correspond to a 1D grandcanonical (blue squares and triangles), a 1D canonical (red circles) and 2D canonical (red stars) regimes.
In order to span the five orders in magnitude in energy, simulations on systems with different tunneling were combined. As can be seen from Eqs. (37)- (39), in the linear regime, the fluctuations only depend on the energy ǫ k and not on the specific value of J. We verified that in the overlapping energy regions, numerical simulations with different J gave the same results. This universality breaks down in the regime of strong fluctuations, where the momentum distribution and static structure factor explicitly depend on J at small energy (compare upright and inverted blue triangles).
In the grandcanonical regime at low energy, the fluctuations become large and the linear approximation breaks down. In analogy with the PJJ, one can use the grandcanonical equilibrium ensemble in order to describe the condensate fluctuations in the grandcanonical regime (see Appendix C). The corresponding blue dashed and dotted lines closely reproduce the numerical results at all energies.
In the large energy limit, we obtain from the nonequilibrium Bogoliubov approximation the limiting expressions In order to gain insight in the behavior of correlations at very large distances, one can start from the low energy dependence in Bogoliubov approximation of the phasephase correlator: The divergence of the phase fluctuations at small energy is reflected in the momentum distribution ( Fig. 4(b)). It appears that at the lowest energies, there is a discrepancy between the results of the numerical simulations and the NEB prediction. Moreover, the numerical results appear to depend on the dimensionality (see the discrepancy between the 1D and 2D results), where the linear theory is dimension independent. The 2D simulations appear to follow the analytical curve somewhat longer, but in both cases, the fluctuations in the numerical simulations are smaller than the analytical prediction. When density fluctuations are small, the phase correlator (43) results in an exponential decay of the first order spatial coherence at large distances, ψ † (x)ψ(x ′ ) ∼ exp(−|x − x ′ |/ℓ c ), with coherence length ℓ c = 4nJ/k B T . This expression coincides with the correlation length of the 1D interacting Bose gas [24]. In equilibrium systems, the condition of small density fluctuations is sat- isfied thanks to interactions, in the present case of noninteracting bosons out of equilibrium, the density fluctuations can be suppressed by losses. This was illustrated in Fig. 4 (red circles, stars and dash-dotted line). If the main contribution to the momentum distribution comes from the large momentum tail (for not too large systems and not too small J), the decay of the first order spatial correlation function is exponential. For a 1D lattice in the limit wherenJ ≫ k B T , the coherence length (in units of the lattice spacing) equals ℓ c = 2nJ/(ηk B T ). In the equilibrium limit η = 1, it reduces to the correlation length of the ideal Bose gas [42] (see Appendix C). From the above analytical considerations, one expects an increase by a factor of two in the correlation length when going from the grandcanonical to the canonical regime. The spatial coherence in real space, obtained from numerical simulations in a 1D array, is shown in Fig. 5(a). From these results obtained for various values of the de- tuning ∆, it is seen that the spatial coherence improves when going from the grandcanonical (red squares) to the canonical (blue triangles) regime. It is also clear from this figure that the increase in the spatial coherence length between grandcanonical (red squares) and canonical (blue triangles) regimes is much larger than the analytically predicted factor of two. The larger than expected coherence reflects the fact that the numerical momentum distribution in the canonical regime in Fig. 4(b) lies below the analytical prediction. The theoretical understanding of the coherence length in the canonical regime requires further investigation. Moreover, in our numerical simulations with small tunneling rates, we have observed that, in analogy with the two-well case, the system can spend a long time in states with large phase difference between neighboring wells. We also defer a study of the occurrence of states with large phase difference to a further study.
In one dimension the suppression of density fluctuations only quantitatively affects the spatial coherence, where in two dimensions the difference is qualitative. Where the ideal 2D bose gas (which has large density fluctuations) does not feature a phase transition to a phase coherent state at finite temperature (the decay of coherence is exponential at all nonzero temperatures), the interacting Bose gas (with suppressed density fluctuations) features a Berezinskii-Kosterlitz-Thouless transition.
In the nonequilibrium case, we have shown that density fluctuations can remain small for noninteracting photons in the canonical regime when losses are present (see Fig.  4(a)). The phase correlator (43) then leads to an algebraic decay of the spatial coherence [24] with the exponent This prediction for the long distance spatial coherence reduces to the equilibrium one when η → 1. Fig. 5(b) shows a comparison of the spatial coherence between a one and a two dimensional array. As predicted by our theoretical analysis, the spatial coherence is better in the two-dimensional case.

VI. CONCLUSIONS AND OUTLOOK
We have introduced a classical model to describe Bose-Einstein condensates of photons in coupled cavities. The model consists of a photon field that is coupled to the molecular states by emission and absorption. The transfer of photons between the cavities is described by the usual nearest neighbor tunneling term. The interplay between tunneling and thermalization through the repeated absorption and emission processes is modeled with an imaginary tunneling term. This term was derived in the approximation that the temperature is much larger than the tunneling energy. If this is not the case, a higher order expansions of the Kennard-Stepanov relation has to be made.
We have studied numerically the full nonlinear equations of motion and obtained analytical expressions from the linearization of the model around a homogeneous steady state, both for two coupled wells and for lattices. From our analytical solutions, we have recovered the equilibrium expressions in the limit of zero cavity losses.
In the grandcanonical ensemble with large density fluctuations, our numerical simulations coincide with the thermal equilibrium results, even when losses are present. In the canonical regime with small density fluctuations, on the other hand, the cavities decouple from each other when the tunneling rate is reduced. This results in a nonmonotoneous dependence of the relative density fluctuations on the tunneling rate.
For one-and two-dimensional lattices of photon condensates in the grandcanonical regime, the spatial coherence reduces to that of the ideal bose gas at thermal equilibrium. In the canonical regime, the spatial coherence markedly improves. In the one-dimensional case, the decay is still exponential, but with a longer coherence length. Where the analytical Bogoliubov analysis predicts an enhancement of the coherence length by a factor of two, the numerical simulations show a much larger coherence length. Further theoretical work will be needed to understand the spatial coherence in the canonical regime. In 2D lattices in the canonical regime, the coherence decays much slower than in 1D. From our Bogoliubov analysis, we predict a power law decay, but the numerics is not conclusive on this point because of the limited lattice sizes that were accessible in our simulations. If the power law coherence exists, this would open the way to the observation of the Berezinksii-Kosterlitz-Thouless transition in (noninteracting) photon condensates, that is stabilized by driving and dissipation [43][44][45][46].
For the PJJ, we have shown that a detuning between the wells is detrimental for the spatial coherence. It will be interesting to analyze the role of disorder in extended lattices as well and study the interplay between Anderson localization and driving and dissipation A further outlook concerns the study of photon BECs in lattices with complex tunneling phases, where the engineering of artificical gauge fields [47] could be possible.

VII. ACKNOWLEDGEMENTS
Stimulating discussions with M. Weitz, J. Klaers, F. Oztürk and W. Verstraelen are warmly acknowledged. This work was financially supported by the grant UA-BOF-FFB150168.

Appendix A: Energy relaxation from coupled bosonic modes
In order to motivate further the introduction of an energy relaxation term for a system with energy-dependent losses, we derive it also for the case of two coupled bosonic modes, where one (ψ) is conserved and one (χ) is dissipative. Due to its Lorentzian spectrum, the losses in the dissipative mode introduces losses in the conserved one that are energy dependent. At the classical field level, this system is described by For a constantT (no kinetic energy) and strong detuning (|ǫ −T | ≫ g, Γ) it has an eigenvalue showing an energy (T )-dependent absorption rate. Expanding aroundT = 0, one has from which one sees that the damping rate has the energy-dependent form γ = γ 0 + 2κT with γ 0 = g 2 Γ ǫ 2 and κ = − Let us now show that this system can be well approximated by a gGPE with energy relaxation parameter κ. The equation of motion for χ can be formally solved as Substituting in Eq. (A1), one obtains After expansion to first order of the denominator in ∂ t and subsequently in Γ, the equation for ψ becomes where in the last line the definitions (A5) of γ 0 and κ were used. We recover here the same relation between the energy dependence of the loss rate and the gGPE as in our derivation based on the KS relation.