Efficient and exact numerical approach for many multi-level systems in open system CQED

In this paper we generalize the numerically exact expansion scheme for many two-level systems coupled to a bosonic (cavity) mode under open system conditions to the case of multi-level systems. For the two-level systems the method is derived from physical and intuitive arguments and the efficient numerical modeling up to hundreds of two-level systems is described. We thereby perform an analysis that allows a natural generalization to multi-level systems. The focus of the paper lies on a comprehensible view on the underlying physics. This facilitates direct application of the method and the use of additional generalizations and approximations.

Omitting the restriction to two levels opens even wider ranges of application: e.g. three-and general multilevel systems are used to study noise induced coherences [13,14], electrically induced transparency [15], and the entangled photon creation via the quantum dot biexciton cascade is usually described by a four-level system [16][17][18][19][20][21]. Three-and four-level systems are also used for more realistic and/or lower threshold laser theories [22,23].
In many applications and experiments ensembles of similar or identical systems are present-for instance more emitters in a laser gain medium lead to an increased output power or irradiating a sample with a diffraction limited laser spot may result in measuring multiple systems at a time. The physical behavior of these ensembles may change drastically compared to the individual systems when quantum correlations become important: e.g. ensembles of cold atoms [24] or quantum dots [25] can exhibit superradiance or ensembles of identical two-level systems interacting with a cavity mode may form a single giant quantum oscillator obeying bose statistics [26], which was experimentally validated in e.g. [27]. In recent years genuine multipartite entanglement has attracted a lot of attention [28,29], which is also closely tied to the occurrence of quantum correlations in ensembles.
In situations where neither one nor many-but few-individual systems are present, standard techniques are not feasible: the full microscopic theory is too elaborate and the approximations invented for many particle systems are generally poor in the few particle case 1 . In addition these systems generally interact with their environment, making an open system description desirable. Mathematically such systems are often treated in some limit (e.g. weak versus strong excitation and also one versus many individual systems), however for a comprehensible understanding it is crucial to be able to treat all-including intermediate-operating regimes.
This creates the need for having a consistent theoretical tool that allows to study single, few and many systems under open system conditions on an equal footing.
In a recent publication a numerically exact and feasible solution for the Markovian quantum master equation of a set of identical two-level emitters was introduced [12]. The solution is based on an expansion of the density matrix in direct product states. Direct product states are products of single system states. Our method provides a significant reduction of the number of degrees of freedom from 4 N to , i.e. from an exponential to a (low) polynomial scaling, without using any approximations [12]. This directly implies the wide range applicability of the method. However, naive implementations of this method lead to issues of numerical convergence, limiting the number of two-level systems to  N 20. We will show that our implementation resolves this problem with minor effort and we are able to treat several hundreds of two-level systems on a single processor.
In the present work we focus on a more exhaustive discussion of the underlying symmetries. We start by giving a detailed introduction into the method for two-level systems. We use the laser as an example to illustrate the mechanisms and physical quantities arising. The concepts introduced in this discussion then facilitate an intuitive generalization to multi-level systems. This generalization is discussed in detail for three-and four-level systems using again the laser example.
In section 2 we will describe the method in detail using N two-level emitters coupled to one photon mode. We start by giving an introduction to the associated Hamiltonian and to Lindblad equation theory for the open system dynamics. Thereafter we construct the basis states, illustrate the simplification procedure and ensure numerical stability. The method is illustrated using the two-level laser example. In section 3 we generalize the method to the case of three-and multi-level systems. Especially in the multi-level emitter case the use of direct product states becomes beneficial since it enables us to identify further reduction of the degrees of freedom and thus further computational speedup. As in the two-level case the method will be illustrated using laser setups.

Two-level systems
In this section we explain the numerical expansion scheme for the Lindblad quantum master equation of a set of N two-level systems in detail [12]. We will show that the exponential scaling of the naive solution 4 N is drastically reduced to~N 3 . The only requirement is that the parameters of the theory-e.g. coupling constants and lifetimes-are identical for the individual two-level systems. Even if this requirement does restrict the application to rather ordered N emitter setups, a solution that is feasible from one up to hundreds of emitters is a significant step, since it allows to study a multitude of different interaction regimes consistently.
Furthermore the method gives equations of motions in a compact form, allowing for an easy implementation even for high emitter numbers.

Open system Dicke model
where we have introduced the bosonic creation and annihilation operators † a , a, which can e.g. describe a single cavity mode. Please note that the theory of the present paper is also valid for more (including offresonant) photon modes. The operators describing the state of the two-level systems are, see figure 1 and [31] where ñ |0 i and ñ |1 i refer to the ground and excited state of the ith two-level system. Thus s i 00 and s i 11 are the projection operators onto the lower and upper state and s s = ( ) † i i 01 10 are the flip operators connected to quantum transitions, see figure 1. We choose this representation as it allows for a clear view on the physical processes 2 . Please note that in (2.1) the energy gap w  1 and the coupling constant g are assumed to be identical for every two-level system. This property is essential for the method introduced in this work.
A common approximation to (2.1) is the rotating wave approximation (RWA) [33] å s This Hamiltonian is usually referred to as Tavis-Cummings Hamiltonian, which is the many emitter generalization of the Jaynes-Cummings model [34]. Our method does not require the RWA, but it can always be included for further simplification. However, in situations where the RWA fails, i.e. in the (ultra) strong coupling regime, one needs be cautious to treat the system-environment coupling correctly [35]. Most equations and sketches in this work are valid for both RWA and Non-RWA-instances where this is not the case will be discussed regarding the respective differences.
In the present work we are interested in the open system dynamics of this model. A popular way to include the effects of environment into the system dynamics is by using a quantum master equation in Born-Markov approximation or rather the Lindblad equation [36], which will be used to introduce the method. Please note that the approach introduced in the following can also be applied to non-Markovian open system theories like Nakajima-Zwanzig [36] and time convolution less theory [36,37].
In open system theories the time evolution of the density matrix instead of state vectors is considered, since the interaction with the environment generally leads to mixed states. The Lindblad equation has the form [36] where the first term on the rhs represents the closed-system, von-Neumann time evolution. The second term has the general form where the A k , † A k are general operators acting on the reduced Hilbert space of the system. These terms are sometimes called Lindblad dissipators and describe the dynamics arising from the Markovian coupling to the environment: i.e. external pumping, dephasing and energy relaxation [36]. The radiative spontaneous emission Liouvillian is given by [38]  where the superscript k on the lhs just indicates that the two-level systems decay individually. Again the inverse lifetime due to spontaneous decay γ should not depend on the index i of the individual two-level systems for the method to be applicable. To construct a theory of optical emitters one could further include a pumping process for the two-level systems and/or a pure dephasing term. Incoherent pumping of the two-level systems can be described by å s r s rs s r rs . In appendix A we give an overview over the Lindblad dissipators that are compatible with the theoretical scheme 3 .
We will introduce the method using a cavity light emission setup (laser) as an example. This setup is chosen because it is well understood, has a multitude of different operating regimes, and was formulated for two-, three-and four-level systems. Our exact solution scheme can be applied in all these cases, illustrating the versatility of the method. Nonetheless we wish to make clear that the current work has a much wider range of application than just laser theory.
Apart from dissipators acting on the two-level systems this theory needs a Liouvillian that describes the cavity photon lifetime. This term is given by where k 1 is the cavity photon lifetime. The total Lindblad equation for a two-level laser theory in the interaction picture reads r r s s We will discuss the solution of the full dynamics within our formalism in section 2.5.
It is beneficial to work with direct product states for the open system dynamics, which will be introduced in the next section. The closed system Dicke model is best described using Dicke states (see appendix B), which are the cause for processes like super-and subradiance [33]. The connection of our method to the Dicke states will be briefly discussed in section 2.6.
Hereafter we interpret the index i as unique label for the individual two-level systems.

Excited state basis
The direct product basis used in the following is usually referred to as excited state basis. The excited state basis of the N two-level system Hilbert space is composed of direct products of the ground and excited states of the individual two-level systems ñ |0 i and ñ |1 i (see figure 1). The two-level systems are assumed to only interact via the cavity mode, i.e. there is no direct coupling between individual two-level systems. For modest two-level system -photon interaction strengths [35] the ground state of the collective system is which is just the direct product of all individual ground states. Higher excited states can be constructed by repeated application of raising operators s j  contains all labels of the excited two-level systems. We use the term label instead of index, since these labels uniquely define one two-level system (see figure 2), a meaning that might be overlooked if we were using the word index. The set containing the labels of all N two-level systems is figure 2. Of course each label can occur only once. Thus the state ñ |n u , n is only one possible realization of a state containing n excited two-level systems. The total number of quantum states with n two-level systems being excited is i.e. the number of possibilities of taking n elements out of a set of N elements. Since  matrix (and thus the Lindblad equation) in this basis is given by , which is the brute force complexity of this problem as stated in the introduction. We will continue our discussion using the basis elements (2.17) and the associated expansion of the density matrix in the form of (2.16), as opposed to (2.15), as it allows for a more compact notation and a more transparent view on the mathematics.
In (2.17) every two-level system is represented by one ket and one bra, i.e. by an operator ñ á | | i j k k for two level system k. There are four possible single two level system operators: ñ á . Based on these basis elements we can make a distinction: let n 11 be the number of two-level systems that are represented by ñ á | | 1 1 k k in (2.17) and let u 11 be the set of the respective labels. Completely analogous we introduce n 10 , n 01 , n 00 and u 10 , u 01 , u 00 for the other three cases. These quantities obey the relations = + + + ( ) N n n n n 2.18 11  Please note that n 00 and u 00 do not enter this expression: as in the definition of the excited state (vector) basis (2.12), we can omit the two-level systems that are in the ground state, because the total number of two-level systems is fixed and thus the information about the ground state two-level systems can be recovered (from (2.18) and (2.19)). An important quantity is the number of basis elements for fixed n 11 , n 10 , n 01 but variable sets u 11 , u 10 , u 01 . It is given by the multinomial coefficient where we have used (2.18) in the last step. The completeness relation in this notation is given by where AE = {} is the empty set and å all comb. runs over all possible sets u. In the next section we look at the time evolution of the density matrix elements associated to the basis in (2.20) and identify symmetries in the action of the quantum master equation.

Time evolution-simplification
We will now expand the quantum master equation in the excited state basis by using where we setÔ to the basis elements defined in (2.20). We will start by discussing the spontaneous emission contribution (see (2.6)) in the quantum master equation (2.23). The whole set of equations arising from the laser example will be discussed in section 2.5. Here we omit the cavity mode as it is not necessary for spontaneous emission into vacuum.
Using  We will start by discussing the mathematical properties of this equation and come back to the physics afterwards. The Liouvillian s [ ]  k 01 has three contributions: in the first term on the rhs of (2.24) the single system operators act on both sides (bra and ket) of the basis element simultaneously, in the remaining two terms they act on only one side each (bra or ket). Due to this the Liouvillian makes a distinction regarding the sets of excited two-level systems: the first term on the rhs of (2.24) is non-zero only if the two-level system k is in the ground state on both sides of the basis element. Similarly, the remaining two terms are non-zero only if the two-level system is excited on the right or on the left side respectively. However this is the same distinction we made earlier when introducing the different sets of (un)excited two-level systems u ,... 11 (see (2.18)-(2.20)). We proceed by writing  is the set u 11 including the additional element k. The rhs of (2.25) consists of three sums over two-level system labels i. These sums have n 00 , + n n 11 10 , and + n n 11 01 summands regardless of the specific sets of unique two-level system labels u ,... 11 . In fact the whole equation does not depend on the specific sets as long as the parameter of the Liouvillian γ is the same for all two-level systems. This holds for all contributions arising from the Lindblad equation. Upon further requiring that at some initial time all density matrix entries of the form á ñá ñ | | n n n n u u u u , , , , 10 11 11 01 01 with the same numbers n 11 , n 10 , n 01 but different sets u ,... 11 are equal, the information about the sets is redundant. This requirement is fulfilled if the system starts in the ground state or a thermal equilibrium state. Other states such as e.g. one specific two-level system excited and the rest unexcited, which is a popular setup for entanglement creation [42,43], may require a more cautious treatment. If the information about the sets is redundant we may omit it and introduce which contains all information about the system. The action of all Liouvillians on the density matrix elements defined by (2.26) (n xy fixed, u xy variable) is the same. The total number of these elements (for n xy fixed) is given by (2.21).
The number of different r [ ] n n n , , 11 10 01 , i.e. the total number degrees of freedom of the system is The same scaling was independently observed in other publications: to our knowledge the first published mention of this scaling was reported by Sarkar and Satchell in 1987 [44], however without providing a concise explanation. In reference to their work, Carmichael provides an explanation of this scaling in his book based on an operator expectation value hierarchy expansion [38]. In 2012 Hartmann found this scaling behavior by exploiting a SU(4) symmetry of the quantum master equation using abstract group theory [32]. This method was also used to calculate the emission from a coupled metal-nanoparticle-N quantum emitter setup [45] and in the context of conditional Ramsey spectroscopy [40]. Our derivation comes from a simple and accessible argument and we believe that the direct product state representation gives a clear view on the underlying physics.
Furthermore the equations of motion are easily derived and we can provide a straightforward generalization to multi-level systems (below). The action of this Liouvillian on the density matrix element r [ ] n n n , , 11 10 01 , containing in-and outscattering contributions, can be understood as follows: the in-scattering contribution (first term) stems from a density matrix element that has a higher number of emitters in the excited state, i.e. + n 1 11 . The out-scattering contribution (second term) is proportional to the same density matrix element r [ ] n n n , , 11 10 01 . Such a process, where a state with a higher excitation number scatters into a state with a lower excitation number, describes an exponential decay of excitation. The part proportional to + n n 10 01 describes the dephasing. Both contributions together describe spontaneous radiative decay. The decrease in excited two-level systems is reflected in the diagonal elements ( = = n n 0 10 01 ) and the loss of coherence in the off-diagonals ( ¹ ¹ n n 0 10 01 ). A sketch in figure 3 illustrates this mechanism. The arrow pointing from the n 11 to the n 00 circle depicts the process in which states with higher numbers n 11 of excited two-level systems (and thus lower n 00 ) decay to the states with higher n 00 (and thus lower n 11 ), which results in a reduction of the excited state population. The arrows pointing from the n 10 and n 01 circles to the outside depict the destruction of quantum coherence.
In principle we could stop here and implement a simulation based on the quantity r [ ] n n n , , 11 10 01 . However at this level there are issues concerning the numerical stability and for the optical emission example (see (2.10)) we also need to include photon number states.

Preventing numerical instability
The density matrix contains all accessible information about a quantum system [46]. However the this information is in general only indirectly accessible to an experiment through measurable observables. Thus in order to make meaningful predictions we need to calculate expectation values of operators. As an example we could calculate the excited state population expectation value s å i i 11 , using (2.22) Figure 3. Sketch illustrating the action of the spontaneous emission Liouvillian. Arrows and corresponding terms have the same color. The green (black) arrow depicts the loss of excitation, states with higher numbers n 11 decay into states with higher numbers n 00 until all probability has decayed to the ground state (i.e. = = n n N 0, 11 00 ). There are two terms responsible for this process-one in-and one out-scattering term-since the total probability has to be conserved. The yellow and purple (gray) arrows depict the dephasing. . This reduces numerical accuracy or makes computation impossible altogether: the larger binomial coefficients grow faster than exponentially in N, thus quickly leaving reasonable number storage formats. Furthermore since the trace of the density matrix is conserved, i.e.
it is clear that the magnitudes of those elements r [ ] n, 0, 0 that are associated with large binomial coefficients become (less than) exponentially small. Multiplying larger than exponentially large numbers with smaller than exponentially small numbers on a finite precision machine results in enormous numerical errors.
Fortunately, this problem can conveniently be circumvented by accounting for the multinomial number of r [ ] ... at the level of the equations of motion. Since all density matrix entries of the form r [ ] n n n , , 11 10 01 are equal, we can simply multiply the respective equations (like (2.28)) by the constant ( )  n n n , , 11 10 01 (see (2.21)). By defining  n n n n n n n n n , , , This expression requires arguments of comparable magnitude only, which greatly improves numerical convergence and accuracy. The physical interpretation of [ ]  n n n , , 11 10 01 is even more accessible than r [ ] n n n , , 11 10 01 : [ ]  n, 0, 0 gives the full probability of finding the system in a state with n two-levels excited and N−n unexcited and for ¹ n n , 0 [ ]  n n n , , 11 10 01 gives the full polarization/transition probability between adjacent states.
In the next section we give the full equations of motion for the two-level laser example (2.10).

Two-level laser theory
In order to discuss the full time evolution of the two-level laser example (see (2.10)) we need to introduce a basis for the photon mode. These are the usual Fock photon number states ñ |k and the basis for the joint system is constructed via , . For reasons of brevity we will use the RWA in this section. Please note that the RWA has no influence on the scaling behavior of the method. The Non-RWA equations simply contain more terms and may require finer time discretization in numerical integration algorithms. Also the ground state may change in systems where the RWA fails-this does not affect the applicability of the method but one needs to be cautious to treat the system environment coupling correctly [35].  where we have used = --n N n n n 00 11 10 01 for visual clarity. The terms are all associated to processes where one excitation in the two-level systems is destroyed (created) while a photon is created (destroyed).
Omitting the RWA would result in twice the number of terms. A sketch of the action of the electron-photon interaction on the two-level systems is given in figure 4 (b): every arrow corresponds to one term in (2.38). We see that density is exchanged between the two levels 0 and 1 via the build up of quantum coherence, as opposed to the action of the Lindblad dissipators, which are inherently incoherent contributions, see figure 4 (a).

Symmetric subspaces-Dicke model
In section 2.1 we stated that the closed Dicke model is conveniently described by the (entangled) Dicke basis states. For open system conditions we develop here a concise solution using (unentangled) direct product states. The interaction with an environment generally destroys entanglement/quantum coherence within the system [36,46]. In this section we briefly discuss the link between the two approaches and elucidate how the  where the sum in the last line again runs over all possible sets u ,... 11 . This expression produces a totally symmetric superposition of all the operators defined by fixed n 11 , n 10 , n 01 .
We now consider the superradiant projector Hence we have found a direct and accessible link between our method and the Dicke basis. The binomial prefactor stems from the normalization of the Dicke states. It can get very small and thus illustrates that the dimension of the symmetric subspace spanned by (2.40) is small compared to the total dimension of the Liouville space.
Interactions with the environment destroy entanglement and Dicke states are entangled states. With this in mind we can look at the action of e.g. the pure dephasing operator on such a superradiant projector , where p 2 is the number specifying how far away from the main diagonal the entry is. A more exhaustive discussion on the connection to the Dicke basis states can be found in [47], which is based on the group theoretic considerations of [32].

Multi-level systems
The dimension of the Liouville space associated to N n-level systems is n N 2 , which makes computation impossible even for small N and n.
In the previous sections we have considered two-level systems. The action of the quantum master equation on certain groups of elements of the density matrix was found to be identical and the number of degrees of freedom was greatly reduced. The arguments we used are quite general and hold for arbitrary multi-level systems. Therefore also for multi-level systems the exponential complexity can be reduced to a polynomial one.
In this section we start by introducing the formal theory for the three-and multi-level system case. Thereafter the application of this method is discussed using three-and four-level laser schemes. In analogy to the notation of the previous section the ground state of the individual systems shall always be the 0 state, i.e. ñ |0 i for the individual unexcited system and thus the collective ground state is given by (2.11).

Three-level systems
The simplest extension of the theory is to consider three-level systems. Three level systems are often divided into the three categories V, Λ, Ξ according to the relative energies of the levels [31]. Three level systems are used for realistic laser theories [22], noise induced coherences [13,14], electrically induced transparency [15], or for more realistic quantum dot models including e.g. the trion state [48,49]. Following the spirit of the previous section we can define the 9 basis elements for the three-level system based on the basis states ñ where the three-level systems in u n 2 and u n 1 are in state ñ |2 and ñ |1 respectively. The two sets u n 2 and u n 1 are disjoint. The corresponding Hamiltonians are where we have included multiple electromagnetic (cavity) modes. Possible Lindblad dissipators are again e.g. spontaneous radiative decay between individual levels where again the parameter needs to be identical for all individual three-level systems. For instance let us consider a Ξ system, where  2 1 and  1 0 are dipole allowed and  2 0is not dipole allowed and the energy eigenvalues for the states are ascending in the order ñ |0 , ñ |1 , and ñ |2 , see g 12 and g 01 in figures 6(a) and 7(b Here n 2 gives the number of three-level systems that are in the Ket state ñ |2 k in above expression, the associated labels are contained in set u n 2 . The associated Bra state can either be á | 2 k , á | 1 k , or á | 0 k . Analogously to the previous section we distinguish these cases and write For completeness we also define the contribution for the ground state population s i 00 , even though we can omit it in the equations of motion as we did in the two-level system case. n 00 and u n 00 can be defined as  The interaction of levels 0 and 1 with a single cavity mode is described by the H 10 Hamiltonian (see (3.3b)) å s 3 . 1 2 k k k 10 10 10 01 In figure 5 the action the of the H 10 Hamiltonian on the three-level systems is sketched. We have additional degrees of freedom compared to the two-level system case and get an additional exchange of quantum coherence. Please note that the additional exchange of coherence is completely decoupled from the other dynamics. This will become important in section 3.3, where we discuss the application of the many multi-level system expansion scheme by looking at a three-(and four-) level laser scheme.

Multi-level systems
There is a large application range for multi-level systems: e.g. four-level systems describe gain in optical devices with reduced thresholds or are model systems to study the quantum-dot biexciton cascade. Generally more realistic descriptions of quantum optical systems such as (coupled) quantum dots require multi-level systems [49,50], which leads to increasingly complex quantum dynamics [51,52]. Figure 5. Sketch of the action of the H 10 three-level system Dicke Hamiltonian. It is apparent that we get, additionally to the two-level system case, contributions exchanging quantum coherence between n 20 , n 21 and n 02 , n 12 , which are completely separated from the other dynamics of the system. We have simplified the representation with respect to figure 4, i.e. we draw single arrows with two heads.
In the RWA this corresponds to two terms in the equations of motion (see (2.38) and figure 4(b)), without RWA there are four. The dots indicate the transition to higher multi-level system dynamics.
For a general + ( ) l 1 -level system we can introduce the basis matrices   In the next section we will discuss this expansion scheme by looking at three-and four-level laser schemes. In particular we will show that there are instances where the dynamics decouple and the complexity of the solution is reduced further without making any approximations.

Three-and four-level optical emitters
Again-as in the two-level case-optical emitters like the laser serve as an illustration. In this section we discuss two different setups for each three-and four-level systems. We omit the cavity loss Liouvillian and the full equations of motion, as they are lengthy and do not benefit the understanding. The full equations are given in appendix C. The discussion in this section will be centered around the sketches introduced above, see figures 3 and 4. Please note that at the level of the sketches there is no difference between RWA and Non-RWA treatments, thus all conclusion of this section are valid in both cases.
We do not discuss some spontaneous emission contributions and all conceivable pure dephasings for reasons of brevity. The scalings/complexities of the solutions of the respective quantum master equation do not change if these contributions were included. The only difference in the results are more terms in the equations of motion, but it does not affect the numerical scaling. We do include the spontaneous emission of the lasing transition into non-lasing modes, as this rate is crucial for a definition of the β factor, a central parameter in laser theory [12,53,54].
We discuss four different setups: coherently and incoherently driven three-and four-level systems, respectively. The density matrix in the four corresponding equations will be labeled with a roman number index to emphasize that they describe different systems. The quantum master equation for an optically pumped threelevel laser scheme (system I) would look like Here H 10 is the interaction with the cavity mode (see (3.12)), H 20 is the optical pumping of the system (see (3.12)   The incoherent setup corresponds to replacing the arrows belonging to the pumping Hamiltonian g 20 in figure  6 (a) with a single arrow pointing upwards (belonging to the pumping Liouvillian). This emphasizes that densities are directly transferred between the levels by the Lindblad dissipator and not via the build up of polarizations as in the H 20 case.
In figure 8 (a) the action of this Liouvillian is sketched. Upon combining all occurring processes in a single sketch we observe that four polarization degrees of freedom (i.e. n n n n , , , 21 ) are completely decoupled from the other system dynamics, see figure 8 (b). This greatly reduces the numerical complexity and the scaling with the system size: coherences can only be driven by coupling to densities. If we start in a state that has no quantum coherences the density matrix entries for ¹ n n n n , , , remain zero throughout the whole Figure 6. Sketches of the processes occurring in the individual multi-level system for (a) three level and (b) four level optical emitters (see (3.25) and (3.30)).
time evolution. Furthermore, regardless of the initial state these elements will be zero in the steady state since they are only dephased. Therefore this quantum master equation for three-level systems is exactly described by the quantity  This scaling is considerably lower than that of the coherently driven three level systems (see (3.25),~N 8 ) and of course also the naive brute force solution (~3 N 2 ). Please bear in mind that this does not correspond to a Hilbert/ Liouville space truncation-we solely observed that some density matrix elements are not connected to the other elements and thus remain zero and we therefore do not need to compute them. We still can calculate all observables and expectation values without any approximation.
With the help of the diagrams we have access to the exact solution and can identify the simplifications without deriving a single equation of motion, emphasizing the usefulness of the diagrams.
The four-level laser theory would be constructed with the master equation for an incoherent pump (system IV). The level scheme of the coherently driven four-level laser is shown in figure 6 (b). The incoherent setup corresponds to replacing the arrows belonging to pumping Hamiltonian g 30 with a single arrow pointing upwards (belonging to an incoherent pump Liouvillian) in figure 6 (b). The corresponding contributions in the diagrams are the purple arrows in figures 9 (a) and 10 (a).
In figures 9 (a) and (b) the processes included in the Hamiltonians and Lindblad dissipators of the four-level laser with coherent pump are sketched. Again we observe a decoupling of the polarization degrees of freedom. The corresponding elements ( )  n ,... ij are zero, see figures 9 (c) and (d) and thus do not need to be computed. Consequently the complexity of this solution does not scale with N 15 but with N 7 instead. Interestingly this is even lower than in the three-level system case: the symmetry of the four levels allows for a better decoupling of polarization degrees of freedom as opposed to the three-level systems, where all polarization degrees of freedom are connected to densities.
Again a further reduction is achieved via an incoherent pump term (see (3.31)). The solution of this setup scales with~N 5 , see figure 10.
In this section we have illustrated how the use of diagrams facilitates the application of the many identical n-level system method. The sketches for the individual contributions in the master equation are easy to draw and sketches for the full quantum master equation provide an easy and quick tool to identify simplifications and to calculate the scaling behavior for a specific setup. Even if the diagrams for the full quantum master equation dynamics may appear complicated at first sight, their construction is straightforward.
Please bear in mind that, from a mathematical point of view, all presented solutions are exact.

Summary and conclusion
In summary we have presented a theory that allows the numerically exact treatment N multi-level systems in CQED setups.
We introduced the method by discussing two-level systems, i.e. the open system version of the Dicke model. We were able to derive the numerically feasible method completely from intuitive and physical arguments. This allows for a clearer view on the physics and also the limitations of the method. By rescaling the equations of motion we were able to considerably increase the numerical stability of the method.
Our approach also allowed for a straightforward generalization to arbitrary multi-level systems. We developed diagrams that intuitively explain the occurring processes and give a clear view on possible, additional simplifications. Especially in the three-and four-level system examples drastic reductions in complexity were found by just drawing these diagrams, a process of a few minutes. This makes the use of the diagrams appealing as one can state the numerical effort for a specific system very fast and without deriving a single equation of motion.
In this section we give a brief introduction to the Dicke basis. For a more thorough discussion we refer to the existing literature e.g. [33].
As stated above the proper choice of basis for the closed system, which is solely governed by the Dicke Hamiltonian (2.1), is the Dicke basis. In the discussion of the Dicke states it is customary to introduce the collective spin operators 11, 10, 01, 00. B.1 The closed system can be expressed solely in terms of the collective operators and the Dicke basis sates emerge as natural basis of the collective spin operators. The Lindblad dissipators can only nontrivially be written solely in terms of collective operators [32], which is why the simple Dicke state description breaks down in the open system case.
For the Dicke basis it is customary to further introduce the total inversion operator å s s 00 11 00 and the total pseudo-spin operator [33] = + + ( ) i.e. has the full complexity of a many three-level system problem, since all polarization degrees of freedom are coupled to density degrees of freedom, see figure 7. We will use the RWA for reasons of brevity.