Photon condensation in circuit QED by engineered dissipation

Photon condensation in circuit QED by engineered dissipation

D Marcos, A Tomadin, S Diehl and P Rabl Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, A-6020 Innsbruck, Austria Institute for Theoretical Physics, University of Innsbruck, A-6020 Innsbruck, Austria david.marcos@uibk.ac.at
July 15, 2019
Abstract

We study photon condensation phenomena in a driven and dissipative array of superconducting microwave resonators. Specifically, we show that by using an appropriately designed coupling of microwave photons to superconducting qubits, an effective dissipative mechanism can be engineered, which scatters photons towards low-momentum states while conserving their number. This mimics a tunable coupling of bosons to a low temperature bath, and leads to the formation of a stationary photon condensate in the presence of losses and under continuous-driving conditions. Here we propose a realistic experimental setup to observe this effect in two or multiple coupled cavities, and study the characteristics of such an out-of-equilibrium condensate, which arise from the competition between pumping and dissipation processes.

1 Introduction

Massive particles and spin systems have been at the center of investigations in quantum many-body physics since the inception of the field. Photonic systems, however, have been the object of interest in this field only much more recently. This is, on the one hand, due to the fact that photons are intrinsically non-interacting, and on the other hand, due to the difficulty to confine and experiment with photons before they decay. One important exception in this context are ideas to produce a Bose-Einstein condensate (BEC) of weakly interacting polaritons [1, 2, 3, 4, 5], whose large photonic component is essential to achieve a low effective mass and therefore transition temperatures which are accessible in solid-state experiments. While current cold-atom experiments offer a well controlled setting to study bosonic quantum gases [6], the investigation of photonic condensates is still an active area of research, stimulated by new experimental directions in the field [7, 8, 9] and controversial discussions on the exact nature of such condensates [10, 11, 12].

More recently, the connection between many-body physics and photons has been addressed again from the perspective of photonic quantum simulators [13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Here, the general goal is to implement strongly interacting many-body systems in a controllable way, to simulate the properties of non-trivial condensed-matter models. This can be achieved, for example, using ideas from cavity quantum electrodynamics (QED) [23], where effective photon-photon interactions can be obtained through the coupling to an intermediary system, such as atoms [24, 25] or solid-state systems [26, 27, 28]. Based on this principle, various schemes for implementing bosonic Hubbard models for photons on a lattice [13, 14], photonic quantum Hall systems [29, 30, 31, 32], and strongly interacting photons in a 1D continuum [33, 34], have been theoretically investigated. Although the experimental implementation of these ideas is still challenging, the development of scalable cavity-QED systems in on-chip devices [35] is rapidly progressing, and analogous (‘circuit QED’) systems in the microwave regime [36, 37, 38, 39, 40, 41] already approach the stage where photonic lattice models can be realized.

A central and still open question in the field of photonic quantum simulators is how to prepare and probe quantum many-body states of light. Because of unavoidable losses and the absence of a chemical potential (arising from an equilibrium particle reservoir), photonic many-body systems must necessarily be studied under driven and non-equilibrium conditions. Therefore, familiar concepts from condensed-matter physics, such as ground-state or equilibrium phase diagrams, are a priori not accessible in photonic many-body systems. In previous works it has been suggested to study, for example, the transient dynamics of an initially pumped system [15, 17, 42], or to use excitations spectroscopy and photon-correlation measurements of a weakly driven system to infer certain properties of the interacting photonic system [34, 43, 44, 45, 46].

In this work we consider a different scenario and study the dissipative dynamics of a strongly driven cavity array in the presence of engineered dissipation [47, 48, 49, 50, 51, 52, 53, 54]. More precisely, we will show how in a circuit QED setting, the coupling of microwave photons to superconducting qubits can be used to scatter photons from high to low momentum states, eventually accumulating photons in a zero-momentum condensate. Previously, a similar scheme has been explored as a dissipative way to prepare a Bose Einstein condensate (BEC) of atoms [48]. In the case of photonics systems, it is essential that this mechanism implements a number-conserving coupling of photons to an effective low-temperature bath, and therefore provides a new approach for the preparation of quasi-equilibrium many-body states in open and driven photonic systems. Here we propose and analyze a proof-of-principle experiment to study condensation of microwave photons in a dissipative cavity array, and describe the properties of such a non-equilibrium BEC, which arise from the interplay between driving, decay, and thermalization processes. Such experiment would realize a controlled setting for the simulation of non-equilibrium condensation phenomena, and more generally, would offer a new route towards preparing stationary states of strongly interacting photons. Moreover, our work shows that the flexibility provided by circuit QED, which so far has mainly been employed to engineer strong coherent interactions, can equally well be used for the design of various non-trivial dissipative couplings, and be applied to simulate driven quantum systems in unconventional environments.

2 Coupled cavity arrays

Figure 1: Schematics of a coupled array of photonic cavities (represented in grey). The photons in each cavity are tunnel-coupled to neighboring sites with amplitude and decay with rate . The cavities are driven with an external coherent field of strength . This external driving compensates losses and ensures a finite stationary photon population.

Fig. 1 illustrates the basic setup for an array of coupled cavities, where each cavity is represented by a single photonic mode of frequency and bosonic annihilation operator . Photons can tunnel between neighboring cavities with hopping amplitude and, in addition, local interactions with two- or many-level systems can be used to induce Kerr-type nonlinearities with an effective photon-photon interaction strength . A generic model for this system is then given by the Bose-Hubbard Hamiltonian (see e.g. Ref. [13])

(1)

Various generalizations of this model have been discussed in the literature and several implementations have been proposed using optical or microwave cavities. The photonic Hubbard model given in Eq. (1) describes a gas of interacting bosons on a lattice, where the hopping competes with on-site interaction . However, in contrast to the cold-atom physics scenario, the bosons here can decay, and under experimentally-relevant conditions (where is the temperature and the Boltzmann constant), the equilibrium state of this model is simply the vacuum state. Therefore, in photonic many-body systems we are mainly interested in the out-of-equilibrium dynamics of in the presence of losses and external driving fields. In particular, in this work we model the resulting dissipative dynamics for the system density operator by a master equation (ME) of the form,

(2)

where . In Eq. (2) the Hamiltonian describes an external driving field of frequency which is used to excite the system, and the second term accounts for photon losses in each cavity with a field decay rate . While a finite driving field is required to counteract the losses, it will in general also compete with and, for strong driving fields, even dominate the system dynamics. Therefore, in previous works it has been suggested to study either the transient dynamics of an initially prepared photonic state [15, 17, 42] (where for times [42] or to use weak excitation spectroscopy [34, 43, 44, 45, 46] () to probe the many-body spectrum of Hamiltonian .

In this work we are interested in the opposite regime of a strongly and continuously driven cavity array. We study the dynamics of this system in the presence of an additional artificial thermalization mechanism, denoted by in Eq. (2). More precisely, we will show below how a non-local coupling of photons to superconducting qubits can be engineered in an array of microwave cavities to implement a dissipative photon scattering process of the form

(3)

The interpretation of this term can be seen best in the case of just two cavities. Then, the first term in Eq. (3) describes the scattering of photons from the asymmetric (energetically higher) mode into the symmetric (energetically lower) mode , while preserving the total photon number. The second term in Eq. (3) describes the reverse process. In the absence of losses, the action of would, for two cavities, drive arbitrary initial photon states into a thermal equilibrium distribution, characterized by a detailed balance between symmetric and antisymmetric modes, with the ratio defining an effective temperature . Similarly, for a whole cavity array, the processes scatter photons to lower momenta and – as previously shown for an analogous system of cold atoms [48] – drive the system towards a condensate of photons in the zero momentum mode, , where is the quasimomentum and the total number of particles. Also in this case the opposite processes can be roughly interpreted as a finite temperature effect, although, as it will be shown below, for and the stationary state of does no longer represent a true thermal equilibrium state.

The competition of with a Kerr-type nonlinearity in the Hamiltonian has already been discussed in the atomic, number-conserving setup [55, 56]. However, in the photonic case, under realistic conditions, competes with external driving fields and intrinsic photon losses as described by the full ME (2). To study this competition more explicitly, we focus in the following on the limit , and consider the experimental scenario where neighboring cavities are driven by coherent fields of alternating amplitude (see Fig. 1). In this case, the accumulation of photons in the mode can be interpreted as a clear signature for photon condensation induced by the dissipative photon-photon interactions in .

3 Physical implementation

Figure 2: Dissipation engineering with circuit QED. a) Implementation of the Hamiltonian in Eq. (4). Two stripline resonators are coupled through a system of two mutually interacting Cooper pair boxes. b) Four-level diagram corresponding to the two coupled qubits, whose dynamics is described by Eq. (6).

In this section we show how the photon scattering mechanism described by Eq. (3) can be implemented in an array of superconducting microwave resonators. The basic setting is illustrated in Fig. 2a for two cavities, each coupled to a nonlinear element, for example a superconducting Cooper pair box (‘charge qubit’). The qubits are placed next to each other to obtain strong electrostatic or magnetic interactions. As we discuss now, this configuration allows us to engineer both nonlinear as well as non-local dissipation processes for photons.

For simplicity we restrict the following analysis to a single block of only two cavities, but the generalization to a whole array of linked cavities is straightforward. The Hamiltonian for this system can be derived from the corresponding equivalent circuit model schematically shown in Fig. 2a. By restricting each resonator to a single mode , , and by approximating each Cooper pair box by a two-level system with ground state and excited state , we obtain

(4)

where . The first part is the free cavity Hamiltonian, while the second term describes the dipole interaction between photons and qubits with strength . Finally, is the Hamiltonian of the two coupled qubits, which we assume of the generic form 111The explicit derivation using circuit theory is omitted here. For a comprehensive review on the topic c.f. [38, 57, 58].

(5)

Here, are the Rabi frequencies of an applied microwave field, which is used to drive the qubits at frequency . The coupling constants and are the strengths of the diagonal and off-diagonal qubit-qubit interactions, respectively. The Hamiltonian (4) can be rewritten in terms of the symmetric and antisymmetric cavity modes, and the qubit states , , , and . In this new basis, choosing , and changing into a frame rotating with , we obtain

(6)

where , , . Apart from the coherent dynamics described by Eq. (6), we include dissipation in the form of intrinsic cavity loss with rate and qubit decay with rate . The full system dynamics is then described by the ME

(7)

and a summary of the energy levels and the most relevant transitions is shown in Fig. 2b.

Our goal in the following is to eliminate the qubit dynamics and derive an effective ME for the cavity modes. Specifically, we are interested in the dissipative two-photon process, where the qubits change from to by absorbing a photon from the antisymmetric cavity mode () and emitting a photon into the symmetric cavity mode () (see Fig. 2b). Since the overall process is in the photon-qubit coupling strength, a general derivation is quite involved, and to make the following discussion more transparent we will now focus explicitly on the hierarchy of energy scales drawn in Fig. 2b. In particular, we assume that , , and are much larger than all the other frequency scales . In this limit, none of the single-photon processes is resonant and we can use a Schrieffer-Wolff transformation to derive an effective two-photon Hamiltonian. We perform the unitary transformation with , and make the ansatz

(8)

We define and as the first and last lines of the Hamiltonian (6), respectively, and choose the coefficients such that . This can be achieved by setting

(9)

In view of the assumptions discussed above, we will use the approximate results , and define the parameter . After this transformation the Hamiltonian (6) reads

(10)

Here is the qubit Hamiltonian written above, with small frequency shifts absorbed into a redefinition of the detunings . The modified cavity Hamiltonian is

(11)

where , , , ,

(12)
(13)

and the average is taken with respect to the stationary qubit state in the limit . For the parameter regime of interest, , and corresponds to the free cavity Hamiltonian with a qubit-mediated tunneling amplitude . Finally, the effective coupling between photons and qubits can be written as

(14)

where . Note that in we have only kept resonant two-photon processes and already omitted terms like and . Both of these terms oscillate with the detuning and will only give small corrections to the results presented below.

The transformation induces a weak mixing between photon and qubit states, which apart from generating new effective interactions also modifies the dissipative couplings. Since we consider , we find that the main result of this mixing is an enhancement of the cavity decay rates . After a transformation of the jump operators in the ME (7), and averaging over the qubit states, we find that for the parameter regime of interest, these rates are

(15)
(16)

In summary, we find that the system dynamics in the new dressed state representation is well described by the ME

(17)

where

(18)
(19)
(20)

In the limit , qubits and cavities are decoupled, and the system relaxes on a timescale into the state , where . For finite , we can use a perturbative projection operator technique to eliminate the qubit degrees of freedom and derive an effective ME for the reduced cavity density operator [59],

(21)

Evaluating this expression we obtain

(22)

where and is an effective photon-photon interaction

(23)

In Eq. (22) we have defined , , and the interactions and , given in terms of the qubit correlation spectra

(24)

which can be calculated using the quantum regression theorem [60].

3.1 Discussion

By ignoring coherent photon-photon interactions for the moment, we see that the first line of the effective cavity dynamics given in Eq. (22) represents the dissipatively coupled cavities as introduced in Eq. (2). In particular, we can write

(25)

By comparison of Eqs. (22) and (3) we identify , with the parameter determining the minimal effective temperature of the engineered two-photon process. For typical parameters, GHz, GHz and GHz, we obtain , and hence the limit considered in most parts of the paper can be implemented to a very good approximation. However, we point out that the ratio can always be increased artificially, for example, by adding an additional coherent or incoherent driving field to populate the state . Our derivation also shows the existence of coherences between and , which lead to additional squeezing terms of order , and that in Eq. (25) are summarized by . For , and in the presence of cavity losses, we expect the influence of these coherences to be negligible, while the relevant effect on the populations is already captured by a finite . Indeed, numerical simulations similar to those discussed in Sec. 4 show no significant differences between the exact ME (22) and our model Liouvillian (3), even up to values of . However, by choosing and , it would be possible to extend our model and study physical effects directly linked to the parameter , such as number-conserving squeezing as previously analyzed in [61].

Figure 3: a) Single and two-photon rates which appear in the effective model of Eq. (22). The photon scattering rate of interest, , becomes significantly larger than the others in a certain region of the parameter space. Here we have chosen MHz, GHz, MHz, MHz, MHz, kHz. b) Figure of merit for the validity of our model for high-momentum modes, , as a function of qubit driving and detuning . Rest of parameters as in (a). c) The same plot as in (b) but for low-momentum modes, .

Apart from the two-photon scattering processes of interest, the coupling to the qubits introduces various imperfections. On the one hand, these are the enhanced cavity decay rates introduced above, and on the other hand, we obtain additional dephasing terms due to fluctuating Stark shifts when the qubit jumps incoherently between states and . Since the ratio scales as , single photon losses can be suppressed in the strong coupling regime, as long as also exceeds the bare decay rate . The ratio is independent of , and depends on the details of the correlation functions (24) and therefore on the parameters , . In Fig. 3a, is compared with the other decoherence rates, and we see that an experimentally feasible parameter range exists, where is the dominant rate. The main corrections arise from the enhanced decay and the dephasing of the antisymmetric mode, while the corresponding rates for the symmetric mode are much smaller. In Fig. 3b and 3c we plot the two quantities , , as a function of the control parameters , . These quantities represent a simple figure of merit for characterizing the validity of our model for low () and high () momentum modes. We see that within a large parameter regime, the number conserving photon scattering rate can exceed all decoherence processes discussed here, and therefore the proposed model in Eq. (2) provides a reliable description of the system dynamics. In particular, this is true for the low momentum (symmetric) modes, where condensation occurs.

Finally, let us briefly comment on the effective photon-photon interactions described by . For a resonant two-photon process, , we find . Moreover, for the typical parameters considered above, and in particular these interactions vanish for . Therefore, in this work we focus exclusively on the purely dissipative cavity dynamics. However, by setting instead, the regime of strong coherent photon-photon coupling can be engineered as well, having .

In summary, we find that the ME (2) can be implemented with superconducting microwave cavities under realistic experimental conditions. Our analysis, which has been detailed here for generic coupled two level systems, can be thus adapted to various nonlinear Josephson devices, where the parameter plays the role of the nonlinearity.

4 Condensation of photons: two cavities

4.1 Semiclassical approximation

In the previous section we have shown how to implement the effective ME (2). Now, we discuss how the scattering between photons induced by the Liouvillian leads to condensation of photons in an open and driven cavity array. In this section we first illustrate the basic process for the minimal case with cavities only. As already mentioned above, we assume , and consider the case where only the antisymmetric mode is excited, i.e. . In the rotating frame with respect to the drive frequency the equations of motion (EOM) for the field expectation values and of the symmetric and antisymmetric cavity modes are

(26)

where is the detuning of the antisymmetric mode with respect to the driving field. For the moment, we assume in Eq. (3) and postpone the discussion of this term to the end of the present section. The EOM for the field expectation values couple to third-order correlation functions, starting a hierarchy of equations that cannot be solved analytically. To truncate the hierarchy, we resort here to a semiclassical approximation in which the state of the system is assumed to be coherent. With this assumption, each normal-ordered correlation function can be readily substituted by the corresponding product of field expectation values. The EOM read

(27)

where and for brevity. The semiclassical approximation, in general, is valid for large photon numbers but it offers a satisfactory explanation of the dynamics of the system for smaller occupation as well, as we verify below with the numerical integration of the EOM (26).

Figure 4: Results (31) of the semiclassical approximation to the EOM (26). a) Population fraction in the symmetric mode as a function of cavity pumping and dissipation strength . b) Population in the symmetric (solid line) and antisymmetric (dashed line) mode for as a function of the pumping strength . The threshold (28) is .

The exact equations (26) become nonlinear -number equations for the expectation values and . Nonlinear equations may have multiple, stable and unstable steady-state solutions. A paradigmatic example of this behavior are the classical equations of the laser [62], where a solution which lacks coherent emission is possible at any pumping strength, but is unstable above a certain threshold. Here we find a comparable result, with a threshold value

(28)

for the pumping strength . Below the threshold, the only stable steady-state solution to Eq. (27) is (c.f. also A)

(29)

and only the pumped antisymmetric mode is populated. For driving strengths above the threshold we obtain the stable solution

(30)

where the symmetric mode features nonvanishing occupation. Note that these equations depend on only via the detuning , which has the only effect to reduce the driving strength. Therefore, from now on we will restrict ourselves to and . The stationary solution to the semiclassical EOM then simplifies to

(31)

The behavior of the populations and is plotted in Fig. 4 as a function of the pumping intensity and the effective decay . Below threshold, the symmetric mode is empty and the antisymmetric mode occupation grows quadratically with the pumping intensity. For , the antisymmetric mode population saturates, and the symmetric mode occupation grows linearly. Thus, if is sufficiently large, most of the photon population is transferred to the symmetric mode. This transition to a state in which photons ‘macroscopically’ occupy the symmetric superposition of two cavity modes is a direct consequence of the engineered nonlinear Liouvillian (3).

The discussion of Eqs. (31) illustrates both, similarities and differences between the present setup and a laser [63, 64, 65, 66]. In a single- or multi-mode laser, the light field is coupled to a pumped reservoir and the lasing threshold corresponds to the point where the total light intensity starts to grow. In our case, since the antisymmetric mode is driven directly by a classical field, the total photon population is always finite and grows monotonically. The threshold then marks the point where the many-body scattering mechanism contained in dominates, and a macroscopic transfer of photons from the antisymmetric to the symmetric mode occurs. In this sense, one can rather think of the antisymmetric mode as a photon reservoir from which a condensation of photons (‘lasing’) into the symmetric mode occurs. However, as we will see below, this analogy between lasing and condensation [10, 11, 12, 67, 68] holds in our system only for certain parameter regimes, and in particular striking differences can be observed in the regime of low photon numbers.

Finally, note that both for as well as for the condensation is suppressed, and that the phase diagram in Fig. 4 is symmetric with respect to the ratio . For small , not enough photons are scattered into the symmetric mode, while for large the antisymmetric mode is highly overdamped and the population of the whole system is small as in the first place.

4.2 Exact diagonalization

Figure 5: Numerical solution to the exact equations of motion for cavities. The occupation numbers and of the symmetric and antisymmetric mode are plotted as a function of the coherent driving strength for a) and b) . In both plots , and the dashed lines indicate the corresponding semiclassical predictions given by Eq. (31). The dotted lines show the two-photon correlation function of the symmetric mode .

To support the results of the semiclassical approximation discussed above, we resort to the exact diagonalization of the system composed by two cavities. The time-evolution of the system for typical parameters features an initial transient in which the photonic population accumulates in the antisymmetric mode, which is directly pumped, and a later stage in which photons scatter into the symmetric mode. Finally, a steady state is reached in which a dynamical equilibrium takes place between the populations in the two modes. We remark that the population in the modes is continuously subjected to losses.

The steady-state values of the populations and in the symmetric and antisymmetric mode, respectively, are shown in Fig. 5 as a function of the pumping strength and for two different values of . In the limit , where we expect the transition to occur at large mode occupation numbers , we find that the exact solution matches qualitatively and quantitatively very well the semiclassical predictions. In particular, we see that above the population grows linearly with the driving strength, while the antisymmetric mode saturates at a value . Differently from the prediction of the semiclassical analysis, however, does not vanish exactly below the threshold but it is smaller than . A study of the equal-time two-photon correlation reveals that, when crossing , the photon statistics of the symmetric mode changes from a thermal () to a Poissonian () distribution. This aspect is in agreement with the standard lasing transition.

In the opposite regime, , Eq. (31) predicts that the population of the driven mode is always less than one, and we do not expect the semiclassical analysis to give accurate results. Indeed, Fig. 5b shows that in this limit the transition is completely washed out and exceeds for all pumping strengths. Nevertheless, above we still observe the characteristic saturation of and – apart from an additional offset – the correct linear scaling of the population in the symmetric mode. However, in strong contrast to the semiclassical regime, we find that the photon statistics of the symmetric mode is sub-Poissonian () at low driving strengths and approaches the classical limit from below. This anti-bunching effect can be understood from the fact that the effective damping of the mode is given by . Therefore, the scattering of the first photon into the symmetric mode changes the damping significantly and suppresses the following repopulation of the ‘reservoir mode’ . The limit and would then result in a ‘single photon condensate’ with and a non-classical statistics with .

4.3 Effective temperature

Finally, let us briefly comment on the effect of the reverse photon scattering process in Eq. (3). As we have discussed in Sec. 2, in the case of two isolated cavities the presence of the reverse photon scattering term can be interpreted as a finite temperature effect, and we can set . By including this term, the EOM are

(32)

These equations are of the same form of (26), with the replacements and . Therefore, in the semiclassical regime we obtain the same physics, but with a renormalized pumping threshold

(33)

In particular, for any fixed driving strength we can use Eq. (33) to define a critical effective temperature above which photon condensation does not occur.

5 Photon condensation in a cavity array

We turn now to the solution of Eq. (2) in a lattice with cavities. The central goal of this section is to show that the many-body Liouvillian (3) induces condensation of photons into the homogeneous lattice mode with zero momentum . Following the scheme introduced in the previous section in the case of cavities only, we consider the scenario in which each cavity mode in the array, with , is coherently driven by a classical field with staggered amplitude . In this way, the coherent driving acts on the edge of the Brillouin zone, while dissipation-induced condensation is expected to take place at the center.

For simplicity, we assume periodic boundary conditions for the lattice. We rewrite Eq. (2) in terms of the annihilation operators of the photonic modes in momentum space, where takes values (given here in units of the inverse lattice spacing) in the discretized Brillouin zone. In the rotating frame with respect to the coherent driving, the Hamiltonian is then given by

(34)

where and is the detuning of the mode from the driving frequency. The Liouville operators are

(35)

and

(36)

In the latter equation, the function plays the role of a scattering kernel, which has the general form

(37)

In the following analysis, we will focus on the case . Furthermore, as discussed above, the photon condensation effect does not depend on the tunneling amplitude and we can assume in the Hamiltonian.

Solving the ME (2) in an extended array is in general a very demanding task, and exact diagonalization can only be provided for very limited system size (see our results in Fig. 6b for cavities). For this reason, we resort in the following to two complementary approximations, which allow us to treat the full dynamics of the photonic population in the whole Brillouin zone. First we consider the equation of motion for the field expectation values within the semiclassical approximation, which was discussed in Sec. 4 for the case , and showed reliable results for large photon numbers. Then we treat the photonic population directly and derive a Boltzmann-like equation starting from the exact EOM for the two-point correlation functions. Within the latter approximation we are able to estimate the finite correlation length of the photonic field in the array. In this respect, this technique fills the gap between the semiclassical analysis of an infinitely extended array and the numerical solution to a small-size system. We emphasize that all techniques agree on the central results that we present here, i.e. the existence of a macroscopic photonic population at the center of the Brillouin zone.

5.1 Semiclassical approximation

Figure 6: a) Numerical solution of the semiclassical equations of motion (38). momentum states are used in the integration and the population is shown in the positive half of the Brillouin zone. b) One-particle density matrix in momentum space, obtained with the numerical integration of the exact equations of motion for cavities, integrated until . In both panels we use and .

For many cavities, the EOM for the cavity fields (analogous to Eq. (27)) read

(38)

where . The numerical solution of this system of nonlinear coupled equations is straightforward. Fig. 6a shows the population on the discretized Brillouin zone for cavities. We see that, in the initial stage of the dynamics, the coherently pumped mode at is populated on the time scale of . Later, photonic population is transferred to the mode