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

The standard example for quantum optical systems consisting of N two-level systems is the Dicke model or Dicke Hamiltonian. It is the theoretical foundation for a large variety of phenomena including lasers, optical parametric amplifiers, superradiance, etc. The Dicke Hamiltonian reads

$$\begin{eqnarray}&&H=\underset{{H}_{0}}{\underbrace{{\rm{\hslash }}{\omega }_{0}{a}^{\dagger }a+{\rm{\hslash }}{\omega }_{1}\displaystyle \sum _{i=1}^{N}{\sigma }_{11}^{i}}}+\underset{{H}_{10}}{\underbrace{{\rm{\hslash }}g(a+{a}^{\dagger })\displaystyle \sum _{i=1}^{N}({\sigma }_{01}^{i}+{\sigma }_{10}^{i})}},\end{eqnarray} \tag{ 2.1 }$$

where we have introduced the bosonic creation and annihilation operators ${a}^{\dagger }$, 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]

$$\begin{eqnarray}{\sigma }_{11}^{i} & = & {| 1\rangle }_{i}{\langle 1| }_{i},\quad {\sigma }_{10}^{i}={| 1\rangle }_{i}{\langle 0| }_{i},\\ {\sigma }_{01}^{i} & = & {| 0\rangle }_{i}{\langle 1| }_{i},\quad {\sigma }_{00}^{i}={| 0\rangle }_{i}{\langle 0| }_{i},\end{eqnarray} \tag{ 2.2 }$$

where ${| 0\rangle }_{i}$ and ${| 1\rangle }_{i}$ refer to the ground and excited state of the ith two-level system. Thus ${\sigma }_{00}^{i}$ and ${\sigma }_{11}^{i}$ are the projection operators onto the lower and upper state and ${\sigma }_{01}^{i}={({\sigma }_{10}^{i})}^{\dagger }$ 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² . Please note that in (2.1) the energy gap ${\rm{\hslash }}{\omega }_{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.

Zoom In Zoom Out Reset image size

A common approximation to (2.1) is the rotating wave approximation (RWA) [33]

$$\begin{eqnarray}&&{H}_{10}\to {H}_{10}^{{\rm{RWA}}}={\rm{\hslash }}g\displaystyle \sum _{i=1}^{N}(a{\sigma }_{10}^{i}+{a}^{\dagger }{\sigma }_{01}^{i}).\end{eqnarray} \tag{ 2.3 }$$

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]

$$\begin{eqnarray}&&\dot{\rho }=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[\rho ,H]+\displaystyle \sum _{k}{ \mathcal D }[{A}_{k}]\rho ,\end{eqnarray} \tag{ 2.4 }$$

where the first term on the rhs represents the closed-system, von-Neumann time evolution. The second term has the general form

$$\begin{eqnarray}&&{ \mathcal D }[{A}_{k}]\rho =\displaystyle \frac{{\gamma }_{k}}{2}(2{A}_{k}\rho {A}_{k}^{\dagger }-{A}_{k}^{\dagger }{A}_{k}\rho -\rho {A}_{k}^{\dagger }{A}_{k}),\end{eqnarray} \tag{ 2.5 }$$

where the A_k, ${A}_{k}^{\dagger }$ 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]

$$\begin{eqnarray}&&{ \mathcal D }[{\sigma }_{01}^{k}]\rho =\displaystyle \frac{\gamma }{2}\displaystyle \sum _{i=1}^{N}(2{\sigma }_{01}^{i}\rho {\sigma }_{10}^{i}-{\sigma }_{11}^{i}\rho -\rho {\sigma }_{11}^{i}),\end{eqnarray} \tag{ 2.6 }$$

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

$$\begin{eqnarray}&&{ \mathcal D }[{\sigma }_{10}^{k}]\rho =\displaystyle \frac{P}{2}\displaystyle \sum _{i=1}^{N}(2{\sigma }_{10}^{i}\rho {\sigma }_{01}^{i}-{\sigma }_{00}^{i}\rho -\rho {\sigma }_{00}^{i}),\end{eqnarray} \tag{ 2.7 }$$

with the generally variable pump rate P. The (phenomenological) pure dephasing contribution is given by

$$\begin{eqnarray}&&{ \mathcal D }[{\sigma }_{z}^{k}]\rho ={\gamma }^{\prime }\displaystyle \sum _{i=1}^{N}({\sigma }_{z}^{i}\rho {\sigma }_{z}^{i}-\rho ),\end{eqnarray} \tag{ 2.8 }$$

with ${\sigma }_{z}={\sigma }_{11}-{\sigma }_{00}$ and ${\sigma }_{z}^{2}=1$. In appendix A we give an overview over the Lindblad dissipators that are compatible with the theoretical scheme³ .

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

$$\begin{eqnarray}&&{ \mathcal D }[b]\rho =\displaystyle \frac{\kappa }{2}(2b\rho {b}^{\dagger }-{b}^{\dagger }b\rho -\rho {b}^{\dagger }b),\end{eqnarray} \tag{ 2.9 }$$

where $1/\kappa$ is the cavity photon lifetime. The total Lindblad equation for a two-level laser theory in the interaction picture reads

$$\begin{eqnarray}&&\dot{\rho }=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[\rho ,{H}_{10}]+{ \mathcal D }[{\sigma }_{01}^{k}]+{ \mathcal D }[{\sigma }_{10}^{k}]+{ \mathcal D }[{\sigma }_{z}^{k}]+{ \mathcal D }[b].\end{eqnarray} \tag{ 2.10 }$$

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.

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\rangle }_{i}$ and ${| 1\rangle }_{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

$$\begin{eqnarray}&&{| 0\rangle }_{N}=\displaystyle \prod _{i=1}^{N}{| 0\rangle }_{i},\end{eqnarray} \tag{ 2.11 }$$

which is just the direct product of all individual ground states. Higher excited states can be constructed by repeated application of raising operators ${\sigma }_{10}^{j}$

$$\begin{eqnarray}&&(\displaystyle \prod _{j\in {{\bf{u}}}_{n}}{\sigma }_{10}^{j})\quad {| 0\rangle }_{N}=\displaystyle \prod _{j\in {{\bf{u}}}_{n}}{| 1\rangle }_{j}\displaystyle \prod _{k\notin {{\bf{u}}}_{n}}{| 0\rangle }_{k}=| n,{{\bf{u}}}_{n}\rangle .\end{eqnarray} \tag{ 2.12 }$$

Here the set ${{\bf{u}}}_{n}=\{{i}_{1},{i}_{2}...{i}_{n}\}$ 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 ${{\bf{u}}}_{N}=\{1,2,...N\}$, see figure 2. Of course each label can occur only once. Thus the state $| n,{{\bf{u}}}_{n}\rangle$ 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

$$\begin{eqnarray}&&{ \mathcal C }(n)=\displaystyle \left(\genfrac{}{}{0em}{}{N}{n}\right),\end{eqnarray} \tag{ 2.13 }$$

i.e. the number of possibilities of taking n elements out of a set of N elements. Since ${H}_{0}| n,{{\bf{u}}}_{n}\rangle ={\rm{\hslash }}{\omega }_{1}n\quad | n,{{\bf{u}}}_{n}\rangle$ we see that ${ \mathcal C }(n)$ is in fact the degeneracy of this state. Accordingly, the total dimension of the Hilbert space is

$$\begin{eqnarray}&&\displaystyle \sum _{n=0}^{N}\displaystyle \left(\genfrac{}{}{0em}{}{N}{n}\right)={2}^{N}.\end{eqnarray} \tag{ 2.14 }$$

Having introduced a vector basis for the Hilbert space of a set of N identical two-level systems (see (2.12)), we can construct the basis elements for the Liouville space of the open system dynamics. The expansion of the density matrix (and thus the Lindblad equation) in this basis is given by

$$\begin{eqnarray}&&\langle m,{{\bf{u}}}_{m}| \rho | n,{{\bf{u}}}_{n}\rangle ,\end{eqnarray} \tag{ 2.15 }$$

which can be rewritten as

$$\begin{eqnarray}\langle m,{{\bf{u}}}_{m}| \rho | n,{{\bf{u}}}_{n}\rangle & = & \mathrm{tr}[| n,{{\bf{u}}}_{n}\rangle \langle m,{{\bf{u}}}_{m}| \rho ]\\ & \equiv & \langle | n,{{\bf{u}}}_{n}\rangle \langle m,{{\bf{u}}}_{m}| \rangle .\end{eqnarray} \tag{ 2.16 }$$

Here the operators

$$\begin{eqnarray}&&| n,{{\bf{u}}}_{n}\rangle \langle m,{{\bf{u}}}_{m}| \end{eqnarray} \tag{ 2.17 }$$

form a complete basis in the Liouville space of the set of N two-level systems. The total number of these operators is ${({2}^{N})}^{2}={4}^{N}$, 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.

Zoom In Zoom Out Reset image size

In (2.17) every two-level system is represented by one ket and one bra, i.e. by an operator ${| i\rangle }_{k}{\langle j| }_{k}$ for two level system k. There are four possible single two level system operators: ${| 1\rangle }_{k}{\langle 1| }_{k}$, ${| 1\rangle }_{k}{\langle 0| }_{k}$, ${| 0\rangle }_{k}{\langle 1| }_{k}$, and ${| 0\rangle }_{k}{\langle 0| }_{k}$, see (2.2). Based on these basis elements we can make a distinction: let n₁₁ be the number of two-level systems that are represented by ${| 1\rangle }_{k}{\langle 1| }_{k}$ in (2.17) and let ${{\bf{u}}}_{11}$ be the set of the respective labels. Completely analogous we introduce n₁₀, n₀₁, n₀₀ and ${{\bf{u}}}_{10}$, ${{\bf{u}}}_{01}$, ${{\bf{u}}}_{00}$ for the other three cases. These quantities obey the relations

$$\begin{eqnarray}&&N={n}_{11}+{n}_{10}+{n}_{01}+{n}_{00}\end{eqnarray} \tag{ 2.18 }$$

and

$$\begin{eqnarray}&&{{\bf{u}}}_{N}={{\bf{u}}}_{11}\cup {{\bf{u}}}_{10}\cup {{\bf{u}}}_{01}\cup {{\bf{u}}}_{00}.\end{eqnarray} \tag{ 2.19 }$$

With this we can write the basis element (2.17) as

$$\begin{eqnarray}| n,{{\bf{u}}}_{n}\rangle \langle m,{{\bf{u}}}_{m}| & = & | {n}_{11}+{n}_{10},{{\bf{u}}}_{11}\cup {{\bf{u}}}_{10}\rangle \langle {n}_{11}+{n}_{01},{{\bf{u}}}_{11}\cup {{\bf{u}}}_{01}| \\ & \equiv & | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| .\end{eqnarray} \tag{ 2.20 }$$

Please note that n₀₀ and ${{\bf{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₁₁, n₁₀, n₀₁ but variable sets ${{\bf{u}}}_{11}$, ${{\bf{u}}}_{10}$, ${{\bf{u}}}_{01}$. It is given by the multinomial coefficient [41]

$$\begin{eqnarray}{ \mathcal C }({n}_{11},{n}_{10},{n}_{01}) & = & \displaystyle \left(\genfrac{}{}{0em}{}{N}{{n}_{11},{n}_{10},{n}_{01},{n}_{00}}\right)\\ & = & \displaystyle \frac{N!}{{n}_{11}!{n}_{10}!{n}_{01}!(N-{n}_{11}-{n}_{10}-{n}_{01})!},\end{eqnarray} \tag{ 2.21 }$$

where we have used (2.18) in the last step. The completeness relation in this notation is given by

$$\begin{eqnarray}&&1=\displaystyle \sum _{n=0}^{N}\displaystyle \sum _{\genfrac{}{}{0em}{}{{\rm{all}}}{{\rm{comb.}}}}| n,{\bf{u}},0,\varnothing \rangle \langle n,{\bf{u}},0,\varnothing | ,\end{eqnarray} \tag{ 2.22 }$$

where $\varnothing =\{\}$ is the empty set and ${\sum }_{\genfrac{}{}{0em}{}{{\rm{all}}}{{\rm{comb.}}}}$ runs over all possible sets ${\bf{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.

We will now expand the quantum master equation in the excited state basis by using

$$\begin{eqnarray}&&{\partial }_{t}\langle \hat{O}\rangle ={\partial }_{t}\mathrm{tr}[\hat{O}\rho ]=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}\mathrm{tr}[\hat{O}[H,\rho ]]+\mathrm{tr}[\hat{O}\displaystyle \sum _{i}{ \mathcal D }[{A}_{i}]\phantom{[\hat{O}}\phantom{\rule{-1.35em}{0ex}}],\end{eqnarray} \tag{ 2.23 }$$

where we set $\hat{O}$ 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 (2.6), (2.20) and (2.23) we write

$$\begin{eqnarray}&&{\partial }_{t}\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle {| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\gamma \displaystyle \sum _{k}[\phantom{\rule[-0000]{0ex}{3.55ex}}\langle {\sigma }_{10}^{k}| {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| {\sigma }_{01}^{k}\rangle \\ &&\qquad -\displaystyle \frac{1}{2}\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| {\sigma }_{11}^{k}\rangle \\ &&\qquad -\displaystyle \frac{1}{2}\langle {\sigma }_{11}^{k}| {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle ]\phantom{\rule[-0000]{0ex}{3.55ex}}.\end{eqnarray} \tag{ 2.24 }$$

We will start by discussing the mathematical properties of this equation and come back to the physics afterwards. The Liouvillian ${ \mathcal D }[{\sigma }_{01}^{k}]$ 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 ${{\bf{u}}}_{11},...$ (see (2.18)–(2.20)).

We proceed by writing

$$\begin{eqnarray}&&{\partial }_{t}\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle {| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\gamma [\phantom{}\displaystyle \sum _{k\in {{\bf{u}}}_{00}}\langle | {n}_{11}+1,{{\bf{u}}}_{11}^{+k},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11}+1,{{\bf{u}}}_{11}^{+k},{n}_{01},{{\bf{u}}}_{01}| \rangle \\ &&\qquad -\displaystyle \frac{1}{2}\displaystyle \sum _{k\in {{\bf{u}}}_{11}\cup {{\bf{u}}}_{10}}\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle \\ &&\qquad -\left.\displaystyle \frac{1}{2}\displaystyle \sum _{k\in {{\bf{u}}}_{11}\cup {{\bf{u}}}_{01}}\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle \right],\end{eqnarray} \tag{ 2.25 }$$

where ${{\bf{u}}}_{11}^{+k}=\{{i}_{1},...,{i}_{{n}_{11}},k\}$ is the set ${{\bf{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₀₀, ${n}_{11}+{n}_{10}$, and ${n}_{11}+{n}_{01}$ summands regardless of the specific sets of unique two-level system labels ${{\bf{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 $\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle$ with the same numbers n₁₁, n₁₀, n₀₁ but different sets ${{\bf{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

$$\begin{eqnarray}&&\langle | {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \rangle \equiv \rho [{n}_{11},{n}_{10},{n}_{01}],\end{eqnarray} \tag{ 2.26 }$$

which contains all information about the system. The action of all Liouvillians on the density matrix elements defined by (2.26) (n_xy fixed, ${{\bf{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 $\rho [{n}_{11},{n}_{10},{n}_{01}]$, i.e. the total number degrees of freedom of the system is

$$\begin{eqnarray}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{n}_{11}+{n}_{10}}{+{n}_{01}+{n}_{00}=N}}1 & = & \displaystyle \frac{1}{6}(N+1)(N+2)(N+3)\\ & = & \displaystyle \frac{1}{6}({N}^{3}+6{N}^{2}+11\;N+6).\end{eqnarray} \tag{ 2.27 }$$

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).

Inserting (2.26) into (2.25) yields

$$\begin{eqnarray}&&{\partial }_{t}\rho [{n}_{11},{n}_{10},{n}_{01}]{| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\displaystyle \frac{\gamma }{2}[2(N-{n}_{11}-{n}_{10}-{n}_{01})\rho [{n}_{11}+1,{n}_{10},{n}_{01}]\\ &&\qquad -(2{n}_{11}+{n}_{10}+{n}_{01})\rho [{n}_{11},{n}_{10},{n}_{01}]]\phantom{\rule[-0000]{0ex}{2.20ex}}.\end{eqnarray} \tag{ 2.28 }$$

The action of this Liouvillian on the density matrix element $\rho [{n}_{11},{n}_{10},{n}_{01}]$, containing in- and out-scattering 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}_{11}+1$. The out-scattering contribution (second term) is proportional to the same density matrix element $\rho [{n}_{11},{n}_{10},{n}_{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}_{10}+{n}_{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}_{10}={n}_{01}=0$) and the loss of coherence in the off-diagonals (${n}_{10}\ne 0\ne {n}_{01}$).

A sketch in figure 3 illustrates this mechanism. The arrow pointing from the n₁₁ to the n₀₀ circle depicts the process in which states with higher numbers n₁₁ of excited two-level systems (and thus lower n₀₀) decay to the states with higher n₀₀ (and thus lower n₁₁), which results in a reduction of the excited state population. The arrows pointing from the n₁₀ and n₀₁ circles to the outside depict the destruction of quantum coherence.

Zoom In Zoom Out Reset image size

In principle we could stop here and implement a simulation based on the quantity $\rho [{n}_{11},{n}_{10},{n}_{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.

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 $\left\langle {\sum }_{i}{\sigma }_{11}^{i}\right\rangle$, using (2.22)

$$\begin{eqnarray}\langle \displaystyle \sum _{i}{\sigma }_{11}^{i}\rangle & = & \langle \displaystyle \sum _{i}{\sigma }_{11}^{i}\displaystyle \sum _{n=0}^{N}\displaystyle \sum _{\genfrac{}{}{0em}{}{{\rm{all}}}{{\rm{comb.}}}}| n,{\bf{u}},0,\varnothing \rangle \langle n,{\bf{u}},0,\varnothing | \rangle \\ & = & \displaystyle \sum _{n=0}^{N}\displaystyle \left(\genfrac{}{}{0em}{}{N}{n}\right)\langle \displaystyle \sum _{i}{\sigma }_{11}^{i}| n,{\bf{u}},0,\varnothing \rangle \langle n,{\bf{u}},0,\varnothing | \rangle \\ & = & \displaystyle \sum _{n=0}^{N}\displaystyle \left(\genfrac{}{}{0em}{}{N}{n}\right)\quad n\quad \rho [n,0,0].\end{eqnarray} \tag{ 2.29 }$$

Even though the information about the sets is redundant there is still a multinomial number of basis elements for each $\rho [i,j,k]$, which is a binomial in this case, since $j=k=0$.

From a numerical point of view this expression is problematic: when evaluating the binomial for e.g. 100 two-level systems we need to calculate and store numbers of dramatically different magnitudes i.e. $\left(\genfrac{}{}{0em}{}{100}{50}\right)\sim {10}^{29}$ versus $\left(\genfrac{}{}{0em}{}{100}{100}\right)=1$. 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.

$$\begin{eqnarray}&&\mathrm{tr}[\rho ]=\displaystyle \sum _{n=0}^{N}\displaystyle \left(\genfrac{}{}{0em}{}{N}{n}\right)\rho [n,0,0]=1,\end{eqnarray} \tag{ 2.30 }$$

it is clear that the magnitudes of those elements $\rho [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 $\rho [...]$ at the level of the equations of motion. Since all density matrix entries of the form $\rho [{n}_{11},{n}_{10},{n}_{01}]$ are equal, we can simply multiply the respective equations (like (2.28)) by the constant ${ \mathcal C }({n}_{11},{n}_{10},{n}_{01})$ (see (2.21)). By defining

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01}]={ \mathcal C }({n}_{11},{n}_{10},{n}_{01})\rho [{n}_{11},{n}_{10},{n}_{01}]\end{eqnarray} \tag{ 2.31 }$$

we can rewrite for instance (2.28)

$$\begin{eqnarray}{\partial }_{t}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01}] & = & \displaystyle \frac{\gamma }{2}[2({n}_{11}+1){ \mathcal P }[{n}_{11}+1,{n}_{10},{n}_{01}]\\ & & -(2{n}_{11}+{n}_{10}+{n}_{01}){ \mathcal P }[{n}_{11},{n}_{10},{n}_{01}]]\phantom{\rule[-0000]{0ex}{2.20ex}}\end{eqnarray} \tag{ 2.32 }$$

and observe that only some prefactors change in (2.28). Upon rewriting the expression for the expectation value (2.29) we observe that this indeed fixes the problem of numerical instability, since

$$\begin{eqnarray}&&\langle \displaystyle \sum _{i}{\sigma }_{11}^{i}\rangle =\displaystyle \sum _{n=0}^{N}n\quad { \mathcal P }[n,0,0].\end{eqnarray} \tag{ 2.33 }$$

This expression requires arguments of comparable magnitude only, which greatly improves numerical convergence and accuracy. The physical interpretation of ${ \mathcal P }[{n}_{11},{n}_{10},{n}_{01}]$ is even more accessible than $\rho [{n}_{11},{n}_{10},{n}_{01}]$: ${ \mathcal P }[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}_{10},{n}_{01}\ne 0$ ${ \mathcal P }[{n}_{11},{n}_{10},{n}_{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).

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\rangle$ and the basis for the joint system is constructed via

$$\begin{eqnarray}&&| k\rangle | n,{{\bf{u}}}_{n}\rangle \to { \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p],\end{eqnarray} \tag{ 2.34 }$$

where k and p are the photon degrees of freedom, i.e. $\langle ...| k\rangle \langle p| \rangle ={ \mathcal P }[...;\;k,p]$.

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].

The contributions of the Lindblad dissipators to the equation of motion are

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{ \mathcal D }[{\sigma }_{10}^{i}]}\\ &&\quad =\displaystyle \frac{P}{2}[2({n}_{00}+1){ \mathcal P }[{n}_{11}-1,{n}_{10},{n}_{01};k,p]\\ &&\qquad -(2{n}_{00}+{n}_{10}+{n}_{01}){ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]]\phantom{\rule[-0000]{0ex}{2.20ex}},\end{eqnarray} \tag{ 2.35 }$$

for the incoherent pumping

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{ \mathcal D }[{\sigma }_{z}^{i}]}=-{\gamma }^{\prime }({n}_{10}+{n}_{01}){ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p],\end{eqnarray} \tag{ 2.36 }$$

for the pure dephasing, and

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{ \mathcal D }[b]}\\ &&\quad =\displaystyle \frac{\kappa }{2}[2\sqrt{(k+1)(p+1)}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k+1,p+1]\\ &&\qquad -(k+p){ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]]\phantom{\rule[-0000]{0ex}{2.20ex}}\end{eqnarray} \tag{ 2.37 }$$

for the cavity photon decay. The action of the three dissipators acting on the two-level systems is sketched in figure 4 (a).

Zoom In Zoom Out Reset image size

The contribution of the interaction Hamiltonian in the RWA ${H}_{10}^{{\rm{RWA}}}$ (see (2.3)) to the quantum master equation is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{H}_{10}^{{\rm{RWA}}}}\\ &&\quad ={\rm{i}}g[({n}_{01}+1)\sqrt{k+1}\quad { \mathcal P }[{n}_{11}-1,{n}_{10},{n}_{01}+1;k+1,p]\\ &&\qquad +({n}_{00}+1)\sqrt{k+1}\quad { \mathcal P }[{n}_{11},{n}_{10}-1,{n}_{01};k+1,p]\\ &&\qquad +({n}_{10}+1)\quad \sqrt{k}\quad \quad \quad { \mathcal P }[{n}_{11},{n}_{10}+1,{n}_{01};k-1,p]\\ &&\qquad +({n}_{11}+1)\quad \sqrt{k}\quad \quad \quad { \mathcal P }[{n}_{11}+1,{n}_{10},{n}_{01}-1;k-1,p]\\ &&\qquad -({n}_{10}+1)\sqrt{p+1}\quad { \mathcal P }[{n}_{11}-1,{n}_{10}+1,{n}_{01};k,p+1]\\ &&\qquad -({n}_{00}+1)\sqrt{p+1}\quad { \mathcal P }[{n}_{11},{n}_{10},{n}_{01}-1;k,p+1]\\ &&\qquad -({n}_{01}+1)\quad \sqrt{p}\quad \quad \quad { \mathcal P }[{n}_{11},{n}_{10},{n}_{01}+1;k,p-1]\\ &&\qquad -({n}_{11}+1)\quad \sqrt{p}\quad \quad \quad { \mathcal P }[{n}_{11}+1,{n}_{10}-1,{n}_{01};k,p-1]]\phantom{\rule[-0000]{0ex}{2.20ex}},\end{eqnarray} \tag{ 2.38 }$$

where we have used ${n}_{00}=N-{n}_{11}-{n}_{10}-{n}_{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).

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 entanglement of Dicke states is reduced under open system conditions.

Looking back at the definition of ${ \mathcal P }[...]$ equation (2.31) we see using equations (2.21) and (2.26)

$$\begin{eqnarray}{ \mathcal P }[{n}_{11},{n}_{10},{n}_{01}] & = & { \mathcal C }({n}_{11},{n}_{10},{n}_{01})\rho [{n}_{11},{n}_{10},{n}_{01}]\\ & = & \mathrm{tr}\left[\rho \displaystyle \sum _{\genfrac{}{}{0em}{}{{\rm{all}}}{{\rm{comb.}}}}| {n}_{11},{{\bf{u}}}_{11},{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{11},{{\bf{u}}}_{11},{n}_{01},{{\bf{u}}}_{01}| \right],\end{eqnarray} \tag{ 2.39 }$$

where the sum in the last line again runs over all possible sets ${{\bf{u}}}_{11},...$. This expression produces a totally symmetric superposition of all the operators defined by fixed n₁₁, n₁₀, n₀₁.

We now consider the superradiant projector

$$\begin{eqnarray}&&\left|\displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right.\rangle \left.\langle \displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right|,\end{eqnarray} \tag{ 2.40 }$$

see appendix B (equations (B.5)–(B.7)) for a definition of the Dicke basis. Here n denotes the number of excited two-level systems, $\tfrac{N}{2}$ is the total pseudo spin, and $n-\tfrac{N}{2}$ is the population inversion. The population inversion is zero if half the two-level systems are excited, negative if less and positive if more than half of the two-level systems are excited. In the case of maximal pseudo spin $l=\tfrac{N}{2}$ the Dicke state is a symmetric superposition of all states of matching population inversion. Upon taking the expectation value of this projector, inserting (B.8), rearranging the indices and inserting (2.39) we observe

$$\begin{eqnarray}&&\mathrm{tr}\left[\left|\displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right.\rangle \left.\langle \displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right|\rho \right]\\ &&\quad =\displaystyle \frac{(N-n)!n!}{N!}\displaystyle \sum _{p=0}^{m}\langle \displaystyle \sum _{\genfrac{}{}{0em}{}{{\rm{all}}}{{\rm{comb.}}}}| n-p,{{\bf{u}}}_{11},p,{{\bf{u}}}_{10}\rangle \langle n-p,{{\bf{u}}}_{11},p,{{\bf{u}}}_{01}| \rangle \\ &&\quad =\displaystyle \frac{(N-n)!n!}{N!}\displaystyle \sum _{p=0}^{m}{ \mathcal P }[n-p,p,p],\qquad m={\rm{min}}(n,N-n).\end{eqnarray} \tag{ 2.41 }$$

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

$$\begin{eqnarray}&&{\partial }_{t}\left\langle \left|\displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right.\rangle \left.\langle \displaystyle \frac{N}{2},n-\displaystyle \frac{N}{2}\right|\right\rangle {| }_{{ \mathcal D }[{\sigma }_{z}^{k}]}\propto -2{\gamma }^{\prime }\displaystyle \sum _{p=0}^{m}p\quad { \mathcal P }[n-p,p,p].\end{eqnarray} \tag{ 2.42 }$$

This constitutes clear view on how the entanglement of a Dicke state is suppressed in the presence of environmental interactions. Different contributions to the entanglement decay on different time scales, just depending on how far away from the main diagonal the corresponding density matrix entry is. The decay rate is $2p{\gamma }^{\prime }$, where $2p$ 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].

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 ${| 0\rangle }_{i}$, ${| 1\rangle }_{i}$ and ${| 2\rangle }_{i}$

$$\begin{eqnarray}{\sigma }_{22}^{i} & = & {| 2\rangle }_{i}{\langle 2| }_{i},\quad {\sigma }_{21}^{i}={| 2\rangle }_{i}{\langle 1| }_{i},\quad {\sigma }_{20}^{i}={| 2\rangle }_{i}{\langle 0| }_{i},\\ {\sigma }_{12}^{i} & = & {| 1\rangle }_{i}{\langle 2| }_{i},\quad {\sigma }_{11}^{i}={| 1\rangle }_{i}{\langle 1| }_{i},\quad {\sigma }_{10}^{i}={| 1\rangle }_{i}{\langle 0| }_{i},\\ {\sigma }_{02}^{i} & = & {| 0\rangle }_{i}{\langle 2| }_{i},\quad {\sigma }_{01}^{i}={| 0\rangle }_{i}{\langle 1| }_{i},\quad {\sigma }_{00}^{i}={| 0\rangle }_{i}{\langle 0| }_{i},\end{eqnarray} \tag{ 3.1 }$$

where the ${\sigma }_{{kk}}^{i}$ are projection and the ${\sigma }_{{kl}}^{i}$ general flip operators. Again the basis states for the collective three-level system can be defined as direct product states

$$\begin{eqnarray}&&| {n}_{2},{{\bf{u}}}_{{n}_{2}},{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle =\displaystyle \prod _{k\in {{\bf{u}}}_{{n}_{2}}}{\sigma }_{20}^{k}\quad \displaystyle \prod _{l\in {{\bf{u}}}_{{n}_{1}}}{\sigma }_{10}^{l}{| 0\rangle }_{N},\end{eqnarray} \tag{ 3.2 }$$

where the three-level systems in ${{\bf{u}}}_{{n}_{2}}$ and ${{\bf{u}}}_{{n}_{1}}$ are in state $| 2\rangle$ and $| 1\rangle$ respectively. The two sets ${{\bf{u}}}_{{n}_{2}}$ and ${{\bf{u}}}_{{n}_{1}}$ are disjoint. The corresponding Hamiltonians are

$$\begin{eqnarray}&&{H}_{0}={\rm{\hslash }}\displaystyle \sum _{l}{\omega }_{l}{b}_{l}^{\dagger }{b}_{l}+{\rm{\hslash }}\displaystyle \sum _{i=0}^{2}{\omega }_{i}\displaystyle \sum _{j}{\sigma }_{{ii}}^{j},\end{eqnarray} \tag{ 3.3a }$$

$$\begin{eqnarray}&&{H}_{{\rm{I}}}={\rm{\hslash }}\displaystyle \sum _{i\ne j}\underset{{H}_{{ij}}}{\underbrace{\displaystyle \sum _{l}\;{g}_{{ij}}^{l}\displaystyle \sum _{k}({\sigma }_{{ij}}^{k}+{\sigma }_{{ji}}^{k})({a}_{l}^{\dagger }+{a}_{l})}},\end{eqnarray} \tag{ 3.3b }$$

where we have included multiple electromagnetic (cavity) modes. Possible Lindblad dissipators are again e.g. spontaneous radiative decay between individual levels ${ \mathcal D }[{\sigma }_{{ij}}^{k}]$ with $i\lt j$ if ${\rm{\hslash }}{\omega }_{i}\lt {\rm{\hslash }}{\omega }_{j}$, where again the parameter needs to be identical for all individual three-level systems. For instance let us consider a Ξ system, where $2\to 1$ and $1\to 0$ are dipole allowed and $2\to 0$ is not dipole allowed and the energy eigenvalues for the states are ascending in the order $| 0\rangle$, $| 1\rangle$, and $| 2\rangle$, see ${\gamma }_{12}$ and ${\gamma }_{01}$ in figures 6(a) and 7(b). Then we would have the spontaneous emission Liouvillians ${ \mathcal D }[{\sigma }_{12}^{k}]$ and ${ \mathcal D }[{\sigma }_{01}^{k}]$. For more information please refer to appendix A.

The associated Liouville space is spanned by the basis elements

$$\begin{eqnarray}&&| {n}_{2},{{\bf{u}}}_{{n}_{2}},{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle \langle {m}_{2},{{\bf{u}}}_{{m}_{2}},{m}_{1},{{\bf{u}}}_{{m}_{1}}| .\end{eqnarray} \tag{ 3.4 }$$

Here n₂ gives the number of three-level systems that are in the Ket state ${| 2\rangle }_{k}$ in above expression, the associated labels are contained in set ${{\bf{u}}}_{{n}_{2}}$. The associated Bra state can either be ${\langle 2| }_{k}$, ${\langle 1| }_{k}$, or ${\langle 0| }_{k}$. Analogously to the previous section we distinguish these cases and write

$$\begin{eqnarray}&&{n}_{2}={n}_{22}+{n}_{21}+{n}_{20},\quad {{\bf{u}}}_{{n}_{2}}={{\bf{u}}}_{{n}_{22}}\cup {{\bf{u}}}_{{n}_{21}}\cup {{\bf{u}}}_{{n}_{20}}.\end{eqnarray} \tag{ 3.5 }$$

For the other contributions we get

$$\begin{eqnarray*}{n}_{1} & = & {n}_{12}+{n}_{11}+{n}_{10},\quad {{\bf{u}}}_{{n}_{1}}={{\bf{u}}}_{{n}_{12}}\cup {{\bf{u}}}_{{n}_{11}}\cup {{\bf{u}}}_{{n}_{10}}\\ {m}_{2} & = & {n}_{22}+{n}_{12}+{n}_{02},\quad {{\bf{u}}}_{{m}_{2}}={{\bf{u}}}_{{n}_{22}}\cup {{\bf{u}}}_{{n}_{12}}\cup {{\bf{u}}}_{{n}_{02}}\\ {m}_{1} & = & {n}_{21}+{n}_{11}+{n}_{01},\quad {{\bf{u}}}_{{m}_{1}}={{\bf{u}}}_{{n}_{21}}\cup {{\bf{u}}}_{{n}_{11}}\cup {{\bf{u}}}_{{n}_{01}}.\end{eqnarray*}$$

For completeness we also define the contribution for the ground state population ${\sigma }_{00}^{i}$, even though we can omit it in the equations of motion as we did in the two-level system case. n₀₀ and ${{\bf{u}}}_{{n}_{00}}$ can be defined as

$$\begin{eqnarray}&&{n}_{00}=N-\displaystyle \sum _{i,j\ne 0}{n}_{{ij}},\quad {{\bf{u}}}_{{n}_{00}}={{\bf{u}}}_{N}\setminus \bigcup _{i,j\ne 0}{{\bf{u}}}_{{ij}},\end{eqnarray} \tag{ 3.6 }$$

where, of course, the indices of summation and union exclude only the case $i=j=0$. The fact that the last contribution n₀₀ can be expressed through the other 8 contributions allows the omission of the index. As long as the number of three- (or general multi-) level systems is fixed we can always eliminate the number connected to one basis element (3.1) and the selection of ${\sigma }_{00}$ is merely convention.

Analogous to (2.18) and (2.19) the following relations are fulfilled

$$\begin{eqnarray}{{\bf{u}}}_{N} & = & {{\bf{u}}}_{22}\cup {{\bf{u}}}_{21}\cup \cdots \cup \;{{\bf{u}}}_{00},\\ N & = & {n}_{22}+{n}_{21}\;+\cdots +\;{n}_{00},\end{eqnarray} \tag{ 3.7 }$$

which just states that the basis elements for the collection of N three-level systems (see (3.4)) consists of N operators of the individual three-level systems (see (3.1))—one for each three-level system.

If all coupling parameters in the quantum master equation do not discriminate between the different systems, i.e. if they are independent of the label of the individual systems, the information about the sets is redundant again. In complete analogy to the previous section closed equations can be obtained by introducing the quantities

$$\begin{eqnarray}&&\quad \mathrm{tr}[| {n}_{2},{{\bf{u}}}_{{n}_{2}},{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle \langle {m}_{2},{{\bf{u}}}_{{m}_{2}},{m}_{1},{{\bf{u}}}_{{m}_{1}}| \rho ]\\ &&\quad \equiv \quad \mathrm{tr}[| {n}_{22},{{\bf{u}}}_{22},{n}_{21},{{\bf{u}}}_{21}...{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{22},...{n}_{01},{{\bf{u}}}_{01}| \rho ]\\ &&\quad \equiv \quad \rho [{n}_{22},{n}_{21},...{n}_{01}]\end{eqnarray} \tag{ 3.8 }$$

and

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{22},{n}_{21},...{n}_{01}]={ \mathcal C }({n}_{22},{n}_{21},...{n}_{01})\rho [{n}_{22},{n}_{21},...{n}_{01}]\end{eqnarray} \tag{ 3.9 }$$

with the multinomial coefficient

$$\begin{eqnarray}&&{ \mathcal C }({n}_{22},{n}_{21},...{n}_{01})=\displaystyle \frac{N!}{{n}_{22}!{n}_{21}!...{n}_{01}!\left(N-{\displaystyle \sum }_{i,j\ne 0}{n}_{{ij}}\right)!},\end{eqnarray} \tag{ 3.10 }$$

where the sum ${\sum }_{i,j\ne 0}$ sums over all contributions except n₀₀. The total number of degrees of freedom of this theory is

$$\begin{eqnarray}&&\displaystyle \sum _{{n}_{22}+...{n}_{00}=N}1=\displaystyle \left(\genfrac{}{}{0em}{}{N+8}{N}\right)\quad \propto {N}^{8}.\end{eqnarray} \tag{ 3.11 }$$

The interaction of levels 0 and 1 with a single cavity mode is described by the H₁₀ Hamiltonian (see (3.3b))

$$\begin{eqnarray}&&{H}_{10}={\rm{\hslash }}{g}_{10}({a}^{\dagger }+a)\displaystyle \sum _{k}({\sigma }_{10}^{k}+{\sigma }_{01}^{k}).\end{eqnarray} \tag{ 3.12 }$$

In figure 5 the action the of the H₁₀ 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.

Zoom In Zoom Out Reset image size

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].

For a general $(l+1)$-level system we can introduce the basis matrices

$$\begin{eqnarray}&&{\sigma }_{{jk}}^{i}={| j\rangle }_{i}{\langle k| }_{i},\quad j,k\in \{0,1,...l\},\end{eqnarray} \tag{ 3.13 }$$

with the projection ${\sigma }_{{kk}}^{i}$ and the general flip ${\sigma }_{{kl}}^{i}$ operators. For a $(l+1)$-level system there are ${(l+1)}^{2}$ different operators of the from (3.13). The Hilbert space of the collection of multi-level systems is spanned by

$$\begin{eqnarray}&&| {n}_{l},{{\bf{u}}}_{{n}_{l}},{n}_{(l-1)},{{\bf{u}}}_{{n}_{(l-1)}}...,{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle =\displaystyle \prod _{k\in {{\bf{u}}}_{{n}_{l}}}{\sigma }_{l0}^{k}...\quad \displaystyle \prod _{m\in {{\bf{u}}}_{{n}_{1}}}{\sigma }_{10}^{m}{| 0\rangle }_{N}.\end{eqnarray} \tag{ 3.14 }$$

The collective basis elements for Liouville space are given by

$$\begin{eqnarray}&&| {n}_{l},{{\bf{u}}}_{{n}_{l}},...,{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle \langle {m}_{l},{{\bf{u}}}_{{m}_{l}},...,{m}_{1},{{\bf{u}}}_{{m}_{1}}| .\end{eqnarray} \tag{ 3.15 }$$

Analogous to two- and three-level systems we divide the numbers ${n}_{l},...$, ${m}_{l},...$ and the associated sets ${{\bf{u}}}_{{n}_{l}},...$, ${{\bf{u}}}_{{m}_{l}},...$ into ${(l+1)}^{2}$ different numbers n_jk and ${(l+1)}^{2}$ different sets ${{\bf{u}}}_{{jk}}$—one number and one set per basis element (see (3.13)):

$$\begin{eqnarray}{n}_{l} & = & {n}_{{ll}}+{n}_{l(l-1)}\cdots {n}_{l0},\quad {{\bf{u}}}_{{n}_{l}}={{\bf{u}}}_{{n}_{{ll}}}\cup {{\bf{u}}}_{{n}_{l(l-1)}}\cdots {{\bf{u}}}_{{n}_{l0}}\\ & & \vdots \\ {n}_{1} & = & {n}_{1l}+{n}_{1(l-1)}\cdots {n}_{10},\quad {{\bf{u}}}_{{n}_{1}}={{\bf{u}}}_{{n}_{1l}}\cup {{\bf{u}}}_{{n}_{1(l-1)}}\cdots {{\bf{u}}}_{{n}_{10}}\\ {m}_{l} & = & {n}_{{ll}}+{n}_{(l-1)l}\cdots {n}_{0l},\quad {{\bf{u}}}_{{m}_{l}}={{\bf{u}}}_{{n}_{{ll}}}\cup {{\bf{u}}}_{{n}_{(l-1)l}}\cdots {{\bf{u}}}_{{n}_{0l}}\\ & & \vdots \\ {m}_{1} & = & {n}_{l1}+{n}_{(l-1)1}\cdots {n}_{01},\quad {{\bf{u}}}_{{m}_{1}}={{\bf{u}}}_{{n}_{l1}}\cup {{\bf{u}}}_{{n}_{(l-1)1}}\cdots {{\bf{u}}}_{{n}_{01}}\end{eqnarray} \tag{ 3.16 }$$

and

$$\begin{eqnarray}&&{n}_{00}=N-\displaystyle \sum _{i,j\ne 0}{n}_{{ij}},\quad {{\bf{u}}}_{{n}_{00}}={{\bf{u}}}_{N}\setminus \bigcup _{i,j\ne 0}{{\bf{u}}}_{{ij}},\end{eqnarray} \tag{ 3.17 }$$

where the summation and union indices again exclude the case $i=j=0$, only. Equivalently

$$\begin{eqnarray}{{\bf{u}}}_{N} & = & {{\bf{u}}}_{{ll}}\cup {{\bf{u}}}_{l(l-1)}\cup \cdots {{\bf{u}}}_{00},\\ N & = & {n}_{{ll}}+{n}_{l(l-1)}+\cdots {n}_{00},\end{eqnarray} \tag{ 3.18 }$$

are also fulfilled. Closed equations can be achieved when expanding the Lindblad equation in

$$\begin{eqnarray}&&\quad \mathrm{tr}[| {n}_{l},{{\bf{u}}}_{{n}_{l}},...{n}_{1},{{\bf{u}}}_{{n}_{1}}\rangle \langle {m}_{l},{{\bf{u}}}_{{m}_{l}},...{m}_{1},{{\bf{u}}}_{{m}_{1}}| \rho ]\\ &&\quad \equiv \quad \mathrm{tr}[| {n}_{{ll}},{{\bf{u}}}_{{ll}},...{n}_{10},{{\bf{u}}}_{10}\rangle \langle {n}_{{ll}},{{\bf{u}}}_{{ll}}...{n}_{01},{{\bf{u}}}_{01}| \rho ]\\ &&\quad \equiv \quad \rho [{n}_{{ll}},{n}_{l(l-1)},...{n}_{01}],\end{eqnarray} \tag{ 3.19 }$$

upon requiring that the parameters of the theory are independent of the individual system labels. Equivalently to the numerically stable ${ \mathcal P }$ for the two-level systems we write

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{{ll}},...{n}_{01}]={ \mathcal C }({n}_{{ll}},...{n}_{01})\rho [{n}_{{ll}},...{n}_{01}],\end{eqnarray} \tag{ 3.20 }$$

where ${ \mathcal C }({n}_{{ll}},...{n}_{01})$ is again a multinomial coefficient

$$\begin{eqnarray}&&{ \mathcal C }({n}_{{ll}},...{n}_{01})=\displaystyle \frac{N!}{{n}_{{ll}}!...{n}_{01}!\left(N-{\displaystyle \sum }_{i,j\ne 0}{n}_{{ij}}\right)!},\end{eqnarray} \tag{ 3.21 }$$

where the sum ${\sum }_{i,j\ne 0}$ sums over all contributions except n₀₀. The number of degrees of freedom of the many $(l+1)$-level system solution is fixed by the number of indices counting the basis elements, i.e. the number of different n_ij. It amounts to the number of possibilities of satisfying the relation

$$\begin{eqnarray}&&N=\underset{m\quad {\rm{summands}}}{\underbrace{\displaystyle \sum _{i,j}{n}_{{ij}}}},\end{eqnarray} \tag{ 3.22 }$$

which is

$$\begin{eqnarray}&&\displaystyle \left(\genfrac{}{}{0em}{}{N+m-1}{N}\right)\quad \propto {N}^{m-1}.\end{eqnarray} \tag{ 3.23 }$$

Hence the full solution of the many $(l+1)$-level system scales as

$$\begin{eqnarray}&&\displaystyle \sum _{0\leqslant {n}_{{nn}}\cdots +{n}_{01}\leqslant N}1=\displaystyle \left(\genfrac{}{}{0em}{}{N+{(l+1)}^{2}-1}{N}\right)\end{eqnarray} \tag{ 3.24 }$$

since ${(l+1)}^{2}$ is the number of basis matrices for the individual $(l+1)$-level system.

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.

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 three-level laser scheme (system I) would look like

$$\begin{eqnarray}&&{\dot{\rho }}_{{\rm{I}}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{{\rm{I}}},{H}_{10}+{H}_{20}]+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{{\rm{I}}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{{\rm{I}}}.\end{eqnarray} \tag{ 3.25 }$$

Here H₁₀ is the interaction with the cavity mode (see (3.12)), H₂₀ is the optical pumping of the system (see (3.12) with subscript 1 replaced by 2), ${ \mathcal D }[{\sigma }_{12}]{\rho }_{{\rm{I}}}$ the incoherent relaxation process that brings occupation in the upper lasing level, and ${ \mathcal D }[{\sigma }_{01}]{\rho }_{{\rm{I}}}$ is the spontaneous emission into non-lasing modes. The level scheme of this setup is shown in figure 6 (a). The associated diagrams are shown in figure 7. The dissipator contributions—figure 7 (b)—connect density degrees of freedom ($\sim {n}_{{ii}}$) with densities, whereas polarization degrees of freedom ($\sim {n}_{{ij}},i\ne j$) are only dephased. At this level the polarizations are decoupled from the densities and one could omit the polarizations entirely, as will become clear below. However the Hamiltonians H₁₀ and H₂₀ connect all polarization degrees of freedom to densities (see figures 7 (a) and (c)), therefore we have to treat the system in its full complexity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This behavior is changed when we replace the pumping Hamiltonian H₂₀ by an incoherent pumping Lindblad dissipator ${ \mathcal D }[{\sigma }_{20}]\rho$, omitting intermediate quantum coherences. The corresponding master equation (system II) is

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{II}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{II}},{H}_{10}]+{ \mathcal D }[{\sigma }_{20}^{k}]{\rho }_{\mathrm{II}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{II}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{II}}.\end{eqnarray} \tag{ 3.26 }$$

The incoherent setup corresponds to replacing the arrows belonging to the pumping Hamiltonian g₂₀ 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₂₀ 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}_{20},{n}_{21},{n}_{02},{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}_{20},{n}_{21},{n}_{02},{n}_{21}\ne 0$ remain zero throughout the whole 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

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01}],\end{eqnarray} \tag{ 3.27 }$$

where now the new relation

$$\begin{eqnarray}&&N={n}_{22}+{n}_{11}+{n}_{10}+{n}_{01}+{n}_{00}\end{eqnarray} \tag{ 3.28 }$$

holds. Using (3.23) we see that this solution scales with

$$\begin{eqnarray}&&\displaystyle \left(\genfrac{}{}{0em}{}{N+4}{N}\right)\sim \displaystyle \frac{1}{4!}{N}^{4}.\end{eqnarray} \tag{ 3.29 }$$

This scaling is considerably lower than that of the coherently driven three level systems (see (3.25), $\sim {N}^{8}$) and of course also the naive brute force solution ($\sim {3}^{2N}$). 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{III}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{III}},{H}_{30}+{H}_{21}]+{ \mathcal D }[{\sigma }_{23}^{k}]{\rho }_{\mathrm{III}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{III}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{III}},\end{eqnarray} \tag{ 3.30 }$$

for a coherent pump (system III) and with

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{IV}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{IV}},{H}_{21}]+{ \mathcal D }[{\sigma }_{30}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{23}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{IV}},\end{eqnarray} \tag{ 3.31 }$$

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₃₀ 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ${ \mathcal P }({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⁷ 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 $\sim {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.

We start by giving the equations of motion associated to (3.25)

$$\begin{eqnarray}&&{\dot{\rho }}_{{\rm{I}}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{{\rm{I}}},{H}_{10}+{H}_{20}]+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{{\rm{I}}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{{\rm{I}}},\end{eqnarray} \tag{ C.1 }$$

which are given in terms the elements

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{22},{n}_{21},{n}_{20},{n}_{12},{n}_{11},{n}_{10},{n}_{02},{n}_{01};k,p],\end{eqnarray} \tag{ C.2 }$$

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. The contribution of the Hamiltonian

$$\begin{eqnarray}&&{H}_{10}^{{\rm{RWA}}}={\rm{\hslash }}{g}_{10}\displaystyle \sum _{k}({a}^{\dagger }{\sigma }_{01}^{k}+a{\sigma }_{10}^{k})\end{eqnarray} \tag{ C.3 }$$

is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{21},{n}_{20},{n}_{12},{n}_{11},{n}_{10},{n}_{02},{n}_{01};\;k,p]{| }_{{H}_{10}^{{\rm{RWA}}}}\\ &&\quad ={\rm{i}}{g}_{10}[({n}_{10}+1)\quad \sqrt{k}\quad \quad \quad \quad { \mathcal P }[...{n}_{10}+1...;\;k-1,p]\\ &&\qquad +({n}_{11}+1)\quad \sqrt{k}\quad \quad \quad \quad { \mathcal P }[...{n}_{11}+1...{n}_{01}-1;\;k-1,p]\\ &&\qquad +({n}_{12}+1)\quad \sqrt{k}\quad \quad \quad \quad { \mathcal P }[...{n}_{12}+1...{n}_{02}-1...;\;k-1,p]\\ &&\qquad +({n}_{00}+1)\quad \sqrt{k+1}\quad { \mathcal P }[...{n}_{10}-1...;\;k+1,p]\\ &&\qquad +({n}_{01}+1)\quad \sqrt{k+1}\quad { \mathcal P }[...{n}_{11}-1...{n}_{01}+1;\;k+1,p]\\ &&\qquad +({n}_{02}+1)\quad \sqrt{k+1}\quad { \mathcal P }[...{n}_{12}-1...{n}_{02}+1...;\;k+1,p]\\ &&\qquad -({n}_{00}+1)\quad \sqrt{p+1}\quad { \mathcal P }[...{n}_{01}-1;\;k,p+1]\\ &&\qquad -({n}_{10}+1)\quad \sqrt{p+1}\quad { \mathcal P }[...{n}_{11}-1,{n}_{10}+1...;\;k,p+1]\\ &&\qquad -({n}_{20}+1)\quad \sqrt{p+1}\quad { \mathcal P }[...{n}_{21}-1,{n}_{20}+1...;\;k,p+1]\\ &&\qquad -({n}_{01}+1)\quad \sqrt{p}\quad \quad \quad \quad { \mathcal P }[...{n}_{01}+1;\;k,p-1]\\ &&\qquad -({n}_{11}+1)\quad \sqrt{p}\quad \quad \quad \quad { \mathcal P }[...{n}_{11}+1,{n}_{10}-1...;\;k,p-1]\\ &&\qquad -({n}_{21}+1)\quad \sqrt{p}\quad \quad \quad \quad { \mathcal P }[...{n}_{21}+1,{n}_{20}-1...;\;k,p-1]].\end{eqnarray} \tag{ C.4 }$$

Here the pumping Hamiltonian will be assumed to be a semiclassical coherent pumping Hamiltonian, again using RWA

$$\begin{eqnarray}&&{H}_{20}^{{\rm{RWA}}}={\rm{\hslash }}\mu \displaystyle \sum _{k}({E}_{t}^{*}{\sigma }_{02}^{k}+{E}_{t}{\sigma }_{20}^{k}),\end{eqnarray} \tag{ C.5 }$$

with

$$\begin{eqnarray}&&{E}_{t}={E}_{0}{{\rm{e}}}^{-{\rm{i}}{\omega }_{\mathrm{ext}}t},\end{eqnarray} \tag{ C.6 }$$

where ${\omega }_{\mathrm{ext}}$ is the frequency and E₀ the amplitude of the semiclassical driving field. This contribution is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{21},{n}_{20},{n}_{12},{n}_{11},{n}_{10},{n}_{02},{n}_{01};k,p]{| }_{{H}_{20}}^{{\rm{RWA}}}\\ &&\quad ={\rm{i}}\mu {E}_{0}[({n}_{20}+1)\quad { \mathcal P }[...{n}_{20}+1...;\;k,p]\\ &&\qquad +({n}_{21}+1)\quad { \mathcal P }[...{n}_{21}+1...{n}_{01}-1;\;k,p]\\ &&\qquad +({n}_{22}+1)\quad { \mathcal P }[...{n}_{22}+1...{n}_{02}-1...;\;k,p]\\ &&\qquad +({n}_{00}+1)\quad { \mathcal P }[...{n}_{20}-1...;\;k,p]\\ &&\qquad +({n}_{01}+1)\quad { \mathcal P }[...{n}_{21}-1...{n}_{01}+1;\;k,p]\\ &&\qquad +({n}_{02}+1)\quad { \mathcal P }[...{n}_{22}-1...{n}_{02}+1...;\;k,p]\\ &&\qquad -({n}_{00}+1)\quad { \mathcal P }[...{n}_{02}-1...;\;k,p]\\ &&\qquad -({n}_{10}+1)\quad { \mathcal P }[...{n}_{12}-1...{n}_{10}+1...;\;k,p]\\ &&\qquad -({n}_{20}+1)\quad { \mathcal P }[...{n}_{22}-1...{n}_{20}+1...;\;k,p]\\ &&\qquad -({n}_{02}+1)\quad { \mathcal P }[...{n}_{02}+1...;\;k,p]\\ &&\qquad -({n}_{12}+1)\quad { \mathcal P }[...{n}_{12}+1...{n}_{10}-1...;\;k,p]\\ &&\qquad -({n}_{22}+1)\quad { \mathcal P }[{n}_{22}+1...{n}_{20}-1...;\;k,p]].\end{eqnarray} \tag{ C.7 }$$

The two spontaneous emission contributions are

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{21},{n}_{20},{n}_{12},{n}_{11},{n}_{10},{n}_{02},{n}_{01};k,p]{| }_{{ \mathcal D }[{\sigma }_{12}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{12}}{2}[2({n}_{22}+1){ \mathcal P }[{n}_{22}+1...{n}_{11}-1...;\;k,p]\\ &&\qquad -(2{n}_{22}+{n}_{12}+{n}_{02}+{n}_{21}+{n}_{20}){ \mathcal P }[...;\;k,p]]\end{eqnarray} \tag{ C.8 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{21},{n}_{20},{n}_{12},{n}_{11},{n}_{10},{n}_{02},{n}_{01};\;k,p]{| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{01}}{2}[2({n}_{11}+1){ \mathcal P }[...{n}_{11}+1...;\;k,p]\\ &&\qquad -(2{n}_{11}+{n}_{12}+{n}_{21}+{n}_{01}+{n}_{10}){ \mathcal P }[...;\;k,p]].\end{eqnarray} \tag{ C.9 }$$

The equations of motion for the second setup (3.26)

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{II}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{II}},{H}_{10}]+{ \mathcal D }[{\sigma }_{20}^{k}]{\rho }_{\mathrm{II}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{II}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{II}},\end{eqnarray} \tag{ C.10 }$$

are given in terms of the elements

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01};k,p],\end{eqnarray} \tag{ C.11 }$$

see figure 8. Again in the RWA, the contribution of the Hamiltonian ${H}_{10}^{{\rm{RWA}}}$ is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01};\;k,p]{| }_{{H}_{10}^{{\rm{RWA}}}}\\ &&\quad ={\rm{i}}{g}_{10}[({n}_{10}+1)\quad \sqrt{k}\quad \quad \quad \quad { \mathcal P }[...{n}_{10}+1...;\;k-1,p]\\ &&\qquad +({n}_{11}+1)\quad \sqrt{k}\quad \quad \quad \quad { \mathcal P }[...{n}_{11}+1...{n}_{01}-1;\;k-1,p]\\ &&\qquad +({n}_{00}+1)\quad \sqrt{k+1}\quad { \mathcal P }[...{n}_{10}-1...;\;k+1,p]\\ &&\qquad +({n}_{01}+1)\quad \sqrt{k+1}\quad { \mathcal P }[...{n}_{11}-1...{n}_{01}+1;\;k+1,p]\\ &&\qquad -({n}_{00}+1)\quad \sqrt{p+1}\quad { \mathcal P }[...{n}_{01}-1;\;k,p+1]\\ &&\qquad -({n}_{10}+1)\quad \sqrt{p+1}\quad { \mathcal P }[...{n}_{11}-1,{n}_{10}+1...;\;k,p+1]\\ &&\qquad -({n}_{01}+1)\quad \sqrt{p}\quad \quad \quad \quad { \mathcal P }[...{n}_{01}+1;\;k,p-1]\\ &&\qquad -({n}_{11}+1)\quad \sqrt{p}\quad \quad \quad \quad { \mathcal P }[...{n}_{11}+1,{n}_{10}-1...;\;k,p-1]].\end{eqnarray} \tag{ C.12 }$$

The three dissipator contributions are

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{ \mathcal D }[{\sigma }_{20}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{20}}{2}[2({n}_{00}+1){ \mathcal P }[{n}_{22}-1...;\;k,p]\\ &&\qquad -(2{n}_{00}+{n}_{01}+{n}_{10}){ \mathcal P }[...;\;k,p]],\end{eqnarray} \tag{ C.13 }$$

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01};\;k,p]{| }_{{ \mathcal D }[{\sigma }_{12}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{12}}{2}[2({n}_{22}+1){ \mathcal P }[{n}_{22}+1...{n}_{11}-1...;\;k,p]\\ &&\qquad -2{n}_{22}{ \mathcal P }[...;\;k,p]]\end{eqnarray} \tag{ C.14 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{22},{n}_{11},{n}_{10},{n}_{01};k,p]{| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{01}}{2}[2({n}_{11}+1){ \mathcal P }[...{n}_{11}+1...;\;k,p]\\ &&\qquad -(2{n}_{11}+{n}_{01}+{n}_{10}){ \mathcal P }[...;\;k,p]].\end{eqnarray} \tag{ C.15 }$$

The equations of motion for (3.30)

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{III}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{III}},{H}_{30}+{H}_{21}]+{ \mathcal D }[{\sigma }_{23}^{k}]{\rho }_{\mathrm{III}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{III}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{III}},\end{eqnarray} \tag{ C.16 }$$

are given in terms of the elements

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};k,p],\end{eqnarray} \tag{ C.17 }$$

see figure 9. The contribution of the Hamiltonian

$$\begin{eqnarray}&&{H}_{21}^{{\rm{RWA}}}={\rm{\hslash }}{g}_{21}\displaystyle \sum _{k}({a}^{\dagger }{\sigma }_{12}^{k}+a{\sigma }_{21}^{k})\end{eqnarray} \tag{ C.18 }$$

is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};k,p]{| }_{{H}_{21}^{{\rm{RWA}}}}\\ &&\quad ={\rm{i}}{g}_{21}[({n}_{12}+1)\sqrt{k+1}\quad { \mathcal P }[...{n}_{22}-1...{n}_{12}+1...;\;k+1,p]\\ &&\qquad +({n}_{11}+1)\sqrt{k+1}\quad { \mathcal P }[...{n}_{21}-1...{n}_{11}+1...;\;k+1,p]\\ &&\qquad +({n}_{21}+1)\quad \sqrt{k}\quad \quad \quad { \mathcal P }[...{n}_{21}+1...{n}_{11}-1...;\;k-1,p]\\ &&\qquad +({n}_{22}+1)\quad \sqrt{k}\quad \quad \quad { \mathcal P }[...{n}_{22}+1...{n}_{12}-1...;\;k-1,p]\\ &&\qquad -({n}_{21}+1)\sqrt{p+1}\quad { \mathcal P }[...{n}_{22}-1,{n}_{21}+1...;\;k,p+1]\\ &&\qquad -({n}_{11}+1)\sqrt{p+1}\quad { \mathcal P }[...{n}_{12}-1,{n}_{11}+1...;\;k,p+1]\\ &&\qquad -({n}_{12}+1)\quad \sqrt{p}\quad \quad \quad { \mathcal P }[...{n}_{12}+1,{n}_{11}-1...;\;k,p-1]\\ &&\qquad -({n}_{22}+1)\quad \sqrt{p}\quad \quad \quad { \mathcal P }[...{n}_{22}+1,{n}_{21}-1...;\;k,p-1]].\end{eqnarray} \tag{ C.19 }$$

The contribution of the Hamiltonian

$$\begin{eqnarray}&&{H}_{30}^{{\rm{RWA}}}={\rm{\hslash }}\mu \displaystyle \sum _{k}({E}_{t}^{*}{\sigma }_{03}^{k}+{E}_{t}{\sigma }_{30}^{k}),\end{eqnarray} \tag{ C.20 }$$

is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};k,p]{| }_{{H}_{30}^{{\rm{RWA}}}}\\ &&\quad ={\rm{i}}\mu {E}_{0}[({n}_{03}+1)\quad { \mathcal P }[...{n}_{33}-1...{n}_{03}+1...;\;k,p]\\ &&\qquad +({n}_{00}+1)\quad { \mathcal P }[...{n}_{30}-1...;\;k,p]\\ &&\qquad +({n}_{30}+1)\quad { \mathcal P }[...{n}_{30}+1...;\;k,p]\\ &&\qquad +({n}_{33}+1)\quad { \mathcal P }[...{n}_{33}+1...{n}_{03}-1...;\;k,p]\\ &&\qquad -({n}_{30}+1)\quad { \mathcal P }[...{n}_{33}-1,{n}_{30}+1...;\;k,p]\\ &&\qquad -({n}_{00}+1)\quad { \mathcal P }[...{n}_{03}-1...;\;k,p]\\ &&\qquad -({n}_{03}+1)\quad { \mathcal P }[...{n}_{03}+1...;\;k,p]\\ &&\qquad -({n}_{33}+1)\quad { \mathcal P }[...{n}_{33}+1,{n}_{30}-1...;\;k,p]].\end{eqnarray} \tag{ C.21 }$$

see figure 9. The three dissipator contributions are

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};k,p]{| }_{{ \mathcal D }[{\sigma }_{23}^{k}]}\\ &&\qquad =\displaystyle \frac{{\gamma }_{23}}{2}[2({n}_{33}+1){ \mathcal P }[{n}_{33}+1...{n}_{22}-1...;\;k,p]\\ &&\qquad \quad -(2{n}_{33}+{n}_{03}+{n}_{30}){ \mathcal P }[...;k,p]],\end{eqnarray} \tag{ C.22 }$$

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};\;k,p]{| }_{{ \mathcal D }[{\sigma }_{12}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{12}}{2}[2({n}_{22}+1){ \mathcal P }[...{n}_{22}+1...{n}_{11}-1...;\;k,p]\\ &&\qquad -(2{n}_{22}+{n}_{12}+{n}_{21}){ \mathcal P }[...;\;k,p]],\end{eqnarray} \tag{ C.23 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{30},{n}_{22},{n}_{21},{n}_{12},{n}_{11},{n}_{03};\;k,p]{| }_{{ \mathcal D }[{\sigma }_{01}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{01}}{2}[2({n}_{11}+1){ \mathcal P }[...{n}_{11}+1...;\;k,p]\\ &&\qquad -(2{n}_{11}+{n}_{12}+{n}_{21}){ \mathcal P }[...;\;k,p]].\end{eqnarray} \tag{ C.24 }$$

The equations of motion for (3.31)

$$\begin{eqnarray}&&{\dot{\rho }}_{\mathrm{IV}}=\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\rho }_{\mathrm{IV}},{H}_{21}]+{ \mathcal D }[{\sigma }_{30}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{23}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{01}^{k}]{\rho }_{\mathrm{IV}}+{ \mathcal D }[{\sigma }_{12}^{k}]{\rho }_{\mathrm{IV}},\end{eqnarray} \tag{ C.25 }$$

are given in terms of the elements

$$\begin{eqnarray}&&{ \mathcal P }[{n}_{33},{n}_{22},{n}_{21},{n}_{12},{n}_{11};\;k,p],\end{eqnarray} \tag{ C.26 }$$

see figure 10. The ${H}_{21}^{{\rm{RWA}}}$, ${ \mathcal D }[{\sigma }_{12}^{k}]$, and ${ \mathcal D }[{\sigma }_{01}^{k}]$ do not change. The pumping dissipator contribution is given by

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{22},{n}_{21},{n}_{12},{n}_{11};\;k,p]{| }_{{ \mathcal D }[{\sigma }_{30}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{30}}{2}[2({n}_{00}+1){ \mathcal P }[{n}_{33}-1...;\;k,p]\\ &&\qquad -2{n}_{00}{ \mathcal P }[...;\;k,p]]\end{eqnarray} \tag{ C.27 }$$

and the ${ \mathcal D }[{\sigma }_{23}^{k}]$ contribution changes slightly since we can neglect the coherences n₀₃ and n₃₀

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal P }[{n}_{33},{n}_{22},{n}_{21},{n}_{12},{n}_{11};k,p]{| }_{{ \mathcal D }[{\sigma }_{23}^{k}]}\\ &&\quad =\displaystyle \frac{{\gamma }_{23}}{2}[2({n}_{33}+1){ \mathcal P }[{n}_{33}+1...{n}_{22}-1...;\;k,p]\\ &&\qquad -2{n}_{33}{ \mathcal P }[...;\;k,p]],\end{eqnarray} \tag{ C.28 }$$