Quantum dynamics of open many-qubit systems strongly coupled to a quantized electromagnetic field in dissipative cavities

Solid-state cavity quantum electrodynamics (QED) attracted much interest as a promising platform for quantum information and quantum sensing systems; see, e.g., thorma2015 ; lodahl2015 ; degen2017 ; dovzhenko2018 ; bitton2019 for recent reviews. A typical scenario involves an ensemble of quantum emitters (ideally, two-level systems) such as quantum dots, defects in crystals, or molecules, strongly coupled to a quantized electromagnetic (EM) field in a dielectric or plasmonic cavity. We will call these quantum emitters qubits for brevity, to be distinguished from the logical qubits which may involve many two-level systems as well as photonic or mixed degrees of freedom. Several or many qubits are required for most applications. Although direct near-field coupling of qubits with each other is possible and desired for some gating protocols, in the nanophotonics context such a coupling would require deterministic placement of qubits with sub-nm accuracy, which is challenging. A simpler scenario which still permits various ways of the quantum state manipulation is the one in which the qubits are coupled only through the common cavity mode(s). This is precisely the situation considered in this paper. We derive step by step analytic solutions for the quantum dynamics of an ensemble of qubits strongly coupled to a quantized cavity mode, in the presence of decoherence and coupling to dissipative reservoirs, starting from the simplest case of $N$ identical qubits, then including a spread of transition frequencies in two or more two-level qubits, and finally the multilevel fermionic systems with arbitrary energy level distribution in the bands. Compact analytic solutions for the quantum states and observables in strongly coupled open quantum systems become possible within the formalism of the stochastic equation for the state vector which we recently developed tokman2020 ; chen2021 .

As we show in the paper, this approach brings significant advantages as compared to the existing methods of treating open quantum systems. Indeed, the Heisenberg-Langevin formalism faces challenges when dealing with nonperturbative dynamics of strongly coupled systems as it leads to nonlinear operator-valued equations. The master equation formalism quickly becomes intractable analytically and allows only numerical treatment as the number of the degrees of freedom increases; see Appendix A for a detailed comparison. Furthermore, determining the measure of entanglement within the master equation formalism is a separate problem since convenient universal figures of merit exist only for simplest systems. In contrast, our approach gives an analytic solution for a nonperturbative evolution of a strongly coupled system with an arbitrary number of the degrees of freedom. Furthermore, the explicit expression for the state vector makes the determination of entanglement or lack thereof much more transparent; it allows one to compare the resulting quantum state with Bell states and other quantum states popular in quantum information applications; it shows how the observables characterizing a given state gets affected by relaxation and eventually drown in an increasing noise, etc. Unlike the Wigner-Weisskopf method, our approach treats the effects of noise correctly and yields a correct equilibrium state at long times.

Of course the possibility of analytic treatment comes with its limitations. We use the rotating wave approximation (RWA) throughout the paper, which means that all frequency scales such as Rabi frequencies and detunings are much smaller than the cavity mode frequency and the optical transition frequencies in the qubits. Beyond RWA, general analytic solutions become impossible in the fully quantized and nonlinear (nonperturbative) regime of light-matter interaction between fermionic and bosonic fields. Note that some recent experiments with low-frequency transitions (e.g. microwave or terahertz) in strongly coupled multi-electron systems went beyond the RWA into the ultra strong coupling regime; e.g., todorov2010 ; forndiaz2017 ; kono2019 .

The simplest situation is when the qubits are identical, for example if the qubit is a particular transition between two electron or vibrational states in a given molecule and the system includes several or many identical molecules. We solve it in Sec. II in the dissipationless case and then in Secs. IV and V in the presence of dissipation and decoherence. It is common knowledge thorma2015 ; dovzhenko2018 ; haroche that a system of $N$ identical qubits responds to the EM field as one qubit with a “giant” collective dipole moment which scales as $\sqrt{N}$. At the same time, the qubits themselves demonstrate a much more complex dynamics. Depending on the initial excited quantum state of the qubits, they may develop either constructive or destructive interference during their evolution, which leads to either a complete energy transfer to the cavity mode or a complete suppression of such transfer. For example, a system in which only one qubit is initially excited may evolve into an entangled state of all qubits which is decoupled from the cavity field. Moreover, in the typical case relevant for the plasmonic nanocavity chikkaraddy2016 ; benz2016 ; park2016 ; pelton2018 ; gross2018 ; park2019 , when the dissipation of a cavity mode is much faster than the decoherence in a qubit, the field dissipation inevitably drives the system into such a decoupled entangled state which remains preserved until the dissipation in an ensemble of qubits kicks in. Note that the existence of entangled dark states has been recently predicted in gray2015 ; gray2016 by calculating the concurrence with the master equation in the Lindblad approximation.

In Secs. III and IV the quantum evolution of two qubits with arbitrary vacuum Rabi frequencies and frequency detunings from the cavity mode are described. The formation of dark states decoupled from the cavity mode due to destructive quantum interference is predicted. In Sec. V we calculate the emission spectra for an ensemble of $N$ qubits in a leaky cavity.

Finally, in Sec. VI, we solve for the evolution of the electron systems in which transition frequencies may differ by a large value, much larger than the homogeneous linewidth of an individual emitter and the linewidth of a cavity mode. This introduces an additional effective decoherence for an ensemble of qubits, in addition to the coupling of electrons and photons to their dissipative reservoirs. At the same time, different transition frequencies introduce extra degrees of freedom and add more functionality to the system. Examples of such systems include semiconductor quantum wells, quantum dots, or defects in crystals.

Since all results in the paper are in the analytic form and the plots are normalized by the single-qubit vacuum Rabi frequency, here we list typical values of the parameters in experimental solid-state nanophotonic systems that determine the strength of light-matter coupling and relaxation rates. To determine if the strong coupling regime and quantum entanglement in the electron-photon system can be achieved, one has to compare the relaxation rates with the characteristic Rabi frequency which determines the oscillations of energy between the cavity field and the qubits. The strong-coupling conditions for an ensemble of $N$ qubits are easier to implement since the effective Rabi frequency scales as $\sqrt{N}$. At the same time, placing a large number of qubits in a nanocavity may present a challenge by itself as individual quantum emitters will experience nonuniform field distribution of a cavity mode and the coupling between them is difficult to control and may lead to additional decoherence.

The parameters mentioned below were obtained from recent experimental reports and reviews such as thorma2015 ; lodahl2015 ; degen2017 ; dovzhenko2018 ; bitton2019 ; todorov2010 ; forndiaz2017 ; kono2019 ; chikkaraddy2016 ; benz2016 ; park2016 ; pelton2018 ; gross2018 ; park2019 ; lukin2016 ; deppe ; reithmaier ; sercel2019 ; fieramosca2019 . This paper does not attempt to provide a comprehensive overview of the rapidly growing literature on the subject.

In electron-based quantum emitters the largest oscillator strengths in the visible/near-infrared range have been observed for excitons in organic molecules, followed by perovskites and more conventional inorganic semiconductor quantum dots. The typical variation of the dipole matrix element of the optical transition which enters the Rabi frequency is from tens of nm to a few Angstrom (in units of the electron charge). Within the same kind of a system, the dipole moment grows with increasing wavelength. The relaxation times are strongly temperature and material quality-dependent, varying from tens or hundreds of ps for electric-dipole allowed interband transitions in quantum dots at 4 K to the tens and hundreds of $\mu$s for spin qubits based on defects in semiconductors and diamond at mK temperatures. Dipole-forbidden optical transitions, e.g. dark excitons in quantum dots, can have the lifetime up to milliseconds, but their coupling to light could be too weak to reach the strong coupling regime. At room temperature the decoherence rates for the optical transitions are in the tens of meV range.

The photon decay times are longest for dielectric micro- and nano-cavities: photonic crystal cavities, nanopillars, distributed Bragg reflector mirrors, microdisk whispering gallery mode cavities, etc. Their quality factors are typically between $10^{3}-10^{7}$, corresponding to photon lifetimes from sub-ns to $\mu$s range. However, the field localization in the dielectric cavities is diffraction-limited, which limits the attainable single-qubit vacuum Rabi frequency values to hundreds of $\mu$eV deppe ; reithmaier . The effective decay rate of the eigenstates (e.g. exciton-polaritons) in dielectric cavity QED systems is typically limited by the relaxation in the fermion quantum emitter subsystem.

In plasmonic cavities, field localization on a nm and even sub-nm scale has been achieved, but the photon decay time is in the ps or even fs range and therefore, the photon losses dominate the overall decoherence rate. Still, when it comes to strong coupling at room temperature to a single quantum emitter such as a single molecule or a quantum dot, the approach utilizing plasmonic nanocavities has seen more success so far. In these systems the Rabi splitting of the order of 100-200 meV has been observed chikkaraddy2016 ; benz2016 ; park2016 ; pelton2018 ; gross2018 ; park2019 .

Consider first a textbook problem of $N$ identical two-level systems with states $\left|0_{j}\right\rangle$ and $\left|1_{j}\right\rangle$, where $j=1,...N$, with energy levels $0$ and $W$. We introduce fermionic operators of annihilation and creation of an excited state $\left|1_{j}\right\rangle$,

$$\hat{\sigma}_{j}=\left|0_{j}\right\rangle\left\langle 1_{j}\right|,\ \ \ \hat{\sigma}_{j}^{\dagger}=\left|1_{j}\right\rangle\left\langle 0_{j}\right|,$$

(1)

the dipole moment operator,

$$\mathbf{\hat{d}}=\sum_{j=1}^{N}\left(\mathbf{d}\hat{\sigma}_{j}^{\dagger}+\mathbf{d^{*}}\hat{\sigma}_{j}\right),$$

(2)

and the Hamiltonian for all qubits,

$$\hat{H}_{a}=W\sum_{j=1}^{N}\hat{\sigma}_{j}^{\dagger}\hat{\sigma}_{j}.$$

(3)

Here $\mathbf{d=}\left\langle 1_{j}\right|\mathbf{\hat{d}}\left|0_{j}\right\rangle$ is the same for all $j$. Our $N$-qubit system interacts with a single-mode field

$$\mathbf{\hat{E}}=\mathbf{E}\left(\mathbf{r}\right)\hat{c}+\mathbf{E}^{\ast}\left(\mathbf{r}\right)\hat{c}^{\dagger},$$

(4)

where $\hat{c}$ and $\hat{c}^{\dagger}$ are standard annihilation and creation operators for bosonic Fock states.

The function $\mathbf{E}\left(\mathbf{r}\right)$ determines the spatial structure of the electric field in a cavity. It is normalized as Tokman2016

$$\int_{V}\frac{\partial\left[\omega^{2}\varepsilon\left(\omega,\mathbf{r}\right)\right]}{\omega\partial\omega}\mathbf{E}^{\ast}\left(\mathbf{r}\right)\mathbf{E}\left(\mathbf{r}\right)d^{3}r=4\pi\hbar\omega$$

(5)

to preserve the standard form of the field Hamiltonian,

$$\hat{H}_{em}=\hbar\omega\left(\hat{c}^{\dagger}\hat{c}+\frac{1}{2}\right).$$

(6)

The relation between the modal frequency $\omega$ and the function $\mathbf{E}\left(\mathbf{r}\right)$ can be found by solving the boundary-value problem of the classical electrodynamics. Here $V$ is a quantization volume and $\varepsilon\left(\omega,\mathbf{r}\right)$ is the dielectric function of a dispersive medium that fills the cavity.

The total Hamiltonian after adding interaction with the field is

$$\hat{H}_{a+em}=\hat{H}_{a}+\hat{H}_{em}-\mathbf{\hat{d}}\cdot\left(\mathbf{E}\hat{c}+\mathbf{E}^{\ast}\hat{c}^{\dagger}\right),$$

(7)

$$\hat{H}=\hbar\omega\left(\hat{c}^{\dagger}\hat{c}+\frac{1}{2}\right)+\sum_{j=1}^{N}\left[W\hat{\sigma}_{j}^{\dagger}\hat{\sigma}_{j}-\left(\mathbf{d}\cdot\mathbf{E}\hat{\sigma}_{j}^{\dagger}\hat{c}+\mathbf{d}^{\ast}\cdot\mathbf{E}\hat{\sigma}_{j}\hat{c}+\mathbf{d}\cdot\mathbf{E}^{\ast}\hat{\sigma}_{j}^{\dagger}\hat{c}^{\dagger}+\mathbf{d}^{\ast}\cdot\mathbf{E}^{\ast}\hat{\sigma}_{j}\hat{c}^{\dagger}\right)\right],$$

(8)

where we assumed that the field amplitude is the same at the position of all qubits. If the field frequency is close to resonance,

$$\frac{W}{\hbar}\approx\omega$$

(9)

we can use the RWA Hamiltonian in which counter-rotating terms are dropped:

$$\hat{H}=\hbar\omega\left(\hat{c}^{\dagger}\hat{c}+\frac{1}{2}\right)+\sum_{j=1}^{N}\left[W\hat{\sigma}_{j}^{\dagger}\hat{\sigma}_{j}-\left(\mathbf{d}\cdot\mathbf{E}\hat{\sigma}_{j}^{\dagger}\hat{c}+\mathbf{d}^{\ast}\cdot\mathbf{E}^{\ast}\hat{\sigma}_{j}\hat{c}^{\dagger}\right)\right],$$

(10)

The Schrödinger equation with the Hamiltonian (10) is separable into low-dimensional blocks and solvable analytically for any excitation energy and for an arbitrary multiphoton state; see the Appendix A for a general treatment. However, the generation of nonclassical multiphoton states is still an experimental challenge and there is no compelling reason in presenting more general but also more cumbersome formulas for multiphoton excitations. That is why we choose the initial conditions for which the transitions to states beyond single-photon excitations are forbidden by energy conservation. These are most popular states widely used in quantum information applications, including for example Bell states and their generalizations to many-qubit systems. Then the states that can be reached as a result of evolution of the system include the ground state $\left|0\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle$ and the states with energies close to the single-photon energy:

$$\Psi=C_{00}\left|0\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle+C_{10}\left|1\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle+\sum_{j=1}^{N}C_{0j}\left|0\right\rangle\left|1_{j}\right\rangle\Pi_{m\neq j}^{N}\left|0_{m}\right\rangle.$$

(11)

Substituting the state vector (11) into Eq. (10) leads to the equations for coefficients,

$$\dot{C}_{00}+i\frac{\omega}{2}C_{00}=0,$$

(12)

$$\dot{C}_{10}+i\frac{3\omega}{2}C_{10}-i\Omega_{R}^{\ast}\sum_{j=1}^{N}C_{0j}=0,$$

(13)

$$\dot{C}_{0j}+i\left(\frac{\omega}{2}+\frac{W}{\hbar}\right)C_{0j}-i\Omega_{R}C_{10}=0,$$

(14)

where the Rabi frequency $\Omega_{R}=\frac{\mathbf{d\cdot E}}{\hbar}$. Let $\sum_{j=1}^{N}C_{0j}=F$. Then the equation for $C_{00}(t)$ is decoupled,

$$C_{00}(t)=C_{00}(0)e^{-i\frac{\omega}{2}t};$$

(15)

and the other two equations,

$$\dot{C}_{10}+i\frac{3\omega}{2}C_{10}-i\Omega_{R}^{\ast}F=0,$$

(16)

$$\dot{F}+i\left(\frac{\omega}{2}+\frac{W}{\hbar}\right)F-iN\Omega_{R}C_{10}=0,$$

(17)

can be easily solved. Indeed, after making the substitution

$$\left(\begin{array}[]{c}C_{10}\\ F\end{array}\right)=\left(\begin{array}[]{c}G_{10}\\ G_{F}\end{array}\right)e^{-i\frac{3\omega}{2}t},$$

we obtain

$$\dot{G}_{10}-i\Omega_{R}^{\ast}G_{F}=0,$$

(18)

$$\dot{G}_{F}+i\Delta G_{F}-i\Omega_{R}NG_{10}=0,$$

(19)

where $\Delta=\frac{W}{\hbar}-\omega.$

This leads to the following equation for the variable $G_{10}$ which describes the quantum state of the field,

$$\ddot{G}_{10}+i\Delta\dot{G}_{10}+N|\Omega_{R}|^{2}G_{10}=0,$$

(20)

which looks exactly like the one for a single qubit, but with an $N$-qubit Rabi frequency obtained by $\left|\Omega_{R}\right|^{2}\Longrightarrow N\left|\Omega_{R}\right|^{2}$. It has the solution

$$G_{10}=A_{1}e^{\Gamma_{1}t}+A_{2}e^{\Gamma_{2}t},$$

where the constants $A_{1,2}$ are determined by initial conditions and

$$\Gamma_{1,2}=-i\frac{\Delta}{2}\pm i\sqrt{\frac{\Delta^{2}}{4}+N\left|\Omega_{R}\right|^{2}}.$$

(21)

Then the solution for coefficients $C_{0j}$ can be obtained by substituting the expression for $C_{10}$ into Eq. (14),

$$C_{0j}(t)=e^{-i\left(\frac{3\omega}{2}+\Delta\right)t}\left(C_{0j}(0)+i\Omega_{R}\int_{0}^{t}e^{i\Delta\tau}G_{10}(\tau)\,d\tau\right).$$

(22)

Therefore, a system of $N$ identical qubits responds to the EM field as one qubit with a “giant” dipole moment $\propto\sqrt{N}$. At the same time, the qubits themselves demonstrate a much more complex dynamics. Their interaction with a quantum field leads to their entanglement not only with the field but also with each other; see, e.g., TMD2015 . These quantum correlations within the ensemble of qubits can profoundly alter the evolution of the whole system. We illustrate it with a few examples.

First, consider our system in the factorized initial state $\Psi(0)=\left|1\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle$, in which the photon mode is excited, all qubits are in their ground states, and there is no entanglement between any degrees of freedom. Considering for simplicity the case $\Delta=0$, the subsequent evolution of the state is given by

$$\Psi=e^{-i\frac{3\omega}{2}t}\left[\cos\left(\sqrt{N}|\Omega_{R}|t\right)\left|1\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle+\frac{ie^{i\theta}}{\sqrt{N}}\sin\left(\sqrt{N}|\Omega_{R}|t\right)\sum_{j=1}^{N}\left|0\right\rangle\left|1_{j}\right\rangle\Pi_{m\neq j}^{N}\left|0_{m}\right\rangle\right],$$

(23)

where $\theta=\mathrm{Arg}[\Omega_{R}]$. At the moments of time $\sqrt{N}|\Omega_{R}|t=\frac{\pi}{2}+\pi M$ the qubits are entangled with each other in a Bell-type state but not entangled with the field, whereas at an arbitrary moment of time, there are entangled both with each other and with the field. The excitation energy periodically oscillates between the cavity mode and an entangled state of the qubits. If the system is open, the excitation will dissipate through relaxation in both photon and qubit subsystems.

Now consider a different factorized initial state $\Psi(0)=\left|0\right\rangle\left|1_{1}\right\rangle\Pi_{j=2}^{N}\left|0_{j}\right\rangle$, in which only one qubit is excited, for definiteness the one with $j=1$. This implies that qubits can be individually addressed, for example by a beam of light or a near-field nanotip with excitation radius much smaller than the distance between the qubits.

Taking again $\Delta=0$ for simplicity, we obtain the solution

$$\begin{array}[]{c}\Psi=e^{-i\frac{3\omega}{2}t}\left[\frac{ie^{-i\theta}}{\sqrt{N}}\sin\left(\sqrt{N}|\Omega_{R}|t\right)\left|1\right\rangle\Pi_{j=1}^{N}\left|0_{j}\right\rangle+\left(1-\frac{1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)}{N}\right)\left|0\right\rangle\left|1_{1}\right\rangle\Pi_{j=2}^{N}\left|0_{j}\right\rangle\right.\\ \left.-\frac{1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)}{N}\sum_{j=2}^{N}\left|0\right\rangle\left|1_{j}\right\rangle\Pi_{m\neq j}^{N}\left|0_{m}\right\rangle\right],\end{array}$$

(24)

It follows from the time-dependent coefficients in Eq. (24) that for $N\gg 1$ the occupation of the initial state given by the second term on the right-hand side of Eq. (24) oscillates around 1 with a small amplitude $\frac{1}{N}$, whereas the occupations of all other states oscillate around zero with equally small amplitudes. This means that the probability of the energy transfer from the excited qubit to the field or to an ensemble of ground-state qubits is about $N$ times smaller than the probability that the initial state is preserved, which scales as $\left|1-\frac{1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)}{N}\right|^{2}\geq\left(1-\frac{2}{N}\right)^{2}$.

This result holds true whether or not we know which particular qubit is excited, as long as the excitation radius remains much smaller than the distance between the qubits. Contrary to famous double-slit experiments, for the qubit excitation the “which qubit” information constitutes a classical uncertainty. Indeed, imagine that we illuminate the qubits with a tightly confined beam from the far field or with a near field of the tip. The answer to the question ”which qubit got excited” depends on the spatial mode structure which is determined by classical Maxwell’s equations. Therefore, the created state will have the structure $\Psi(0)=\left|0\right\rangle\left|1_{j}\right\rangle\Pi_{q\neq j}^{N}\left|0_{q}\right\rangle$, where the classical uncertainty in the value of the qubit number $j$ does not give rise to the entangled quantum state of the type $\Psi(0)=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}\left|0\right\rangle\left|1_{j}\right\rangle\Pi_{q\neq j}^{N}\left|0_{q}\right\rangle$ (as quantum uncertainty would).

This reasoning can be generalized to the situation in which a group of $J$ qubits fits inside the excitation radius, where $J\ll N$. We show in Sec. IV below that if a group of $J\ll N$ qubits is initially excited, their deexcitation will be suppressed too and the initial state will be preserved with probability $\sim 1-\frac{J}{N}$. Here again, the uncertainty of the type “which group is illuminated” is classical. At the same time, for a beam carrying single-photon energy there is a quantum uncertainty “which qubit inside the group will get excited”.

The presence of $N>1$ qubits allows the existence of entangled qubit states that are completely decoupled from the cavity field. Indeed, for the initial condition $C_{10}=0$ we can find any number of entangled qubit states satisfying the equation $\left|C_{00}\right|^{2}+\sum_{j=1}^{N}\left|C_{0j}\right|^{2}=1$, $\sum_{j=1}^{N}C_{0j}=0$, i.e. the states that are decoupled from the field. This solution is of the type

$$\left(\begin{array}[]{c}C_{10}\\ \vdots\\ C_{0j}\\ \vdots\end{array}\right)\propto e^{-i\frac{3\omega}{2}t}\left(\begin{array}[]{c}0\\ \vdots\\ A_{0j}e^{-i\Delta t}\\ \vdots\end{array}\right),\ \ \ \sum_{j=1}^{N}A_{0j}=0.$$

The partial and complete decoupling of qubits from the field can be illustrated by considering initial Bell-type states of the two qubits:

$$\Psi_{\pm}(0)=\left|0\right\rangle\frac{\left|1_{1}\right\rangle\left|0_{2}\right\rangle\pm\left|0_{1}\right\rangle\left|1_{2}\right\rangle}{\sqrt{2}}\Pi_{j=3}^{N}\left|0_{j}\right\rangle.$$

(25)

Consider again the exact resonance for simplicity, $\Delta=0$. If the field is not excited initially, i.e., $C_{10}(0)=0$, the solution for the complex probability amplitudes takes the form

$$C_{10}(t)=ie^{-i\theta-i\frac{3\omega}{2}t}\frac{F(0)}{\sqrt{N}}\sin\left(\sqrt{N}|\Omega_{R}|t\right),$$

(26)

$$C_{0j}(t)=e^{-i\frac{3\omega}{2}t}\left(C_{0j}(0)-\frac{F(0)}{N}\left[1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)\right]\right),$$

(27)

where $F(0)=\sum_{j=1}^{N}C_{0j}(0)$. In our case we have $C_{01}(0)=\frac{1}{\sqrt{2}}$ and $C_{02}(0)=\pm\frac{1}{\sqrt{2}}$ for $\Psi_{\pm}(0)$ respectively; $C_{0;j>2}=0$. Therefore, for $\Psi_{(-)}(0)$ we obtain $F(0)=0$, i.e., $\Psi_{(-)}(t)=e^{-i\frac{3\omega}{2}t}\Psi_{(-)}(0)$ and the state of two qubits is completely uncoupled from the field.

For $\Psi_{(+)}(0)$ we obtain $F(0)=\sqrt{2}$ and

$$C_{10}(t)=ie^{-i\theta-i\frac{3\omega}{2}t}\frac{\sqrt{2})}{\sqrt{N}}\sin\left(\sqrt{N}|\Omega_{R}|t\right),$$

(28)

$$C_{0;1,2}(t)=e^{-i\frac{3\omega}{2}t}C_{0;1,2}(0)\left(1-\frac{2}{N}\left[1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)\right]\right),$$

(29)

$$C_{0;j>2}(t)=C_{0;1,2}(0)\frac{2}{N}\left[1-\cos\left(\sqrt{N}|\Omega_{R}|t\right)\right].$$

(30)

Obviously, for $\Psi_{(+)}(0)$ the probability amplitudes are changing but the change is small as long as $\frac{2}{N}\ll 1$.

To summarize, if no qubits are excited in the initial state, their ensemble interacts with a photon as one giant qubit. If the system is in the initial state in which one qubit or a small number of qubits are initially excited, the energy transfer of its excitation to other qubits and the photon mode will be slowed down by a factor of $N$. This will significantly slow down the relaxation of the initial excitation if the relaxation in the qubit subsystem is slower than that for the cavity mode. Moreover, there are entangled qubit states that are completely decoupled from the cavity field. For a solid state system with large $N$, the excitation and control of such or any other specified entangled multi-qubit states can be tricky. However, as we show below, the EM field dissipation may drive the multi-qubit system towards such states.

The monolithic solid-state systems where a very large $N$ can be reached include, for example, intersubband transitions in doped quantum wells (QWs), excitons in semiconductor QWs or 2D semiconductors, and a 2D electron gas in a strong magnetic field. In these systems, strong and even ultra-strong coupling to a cavity field has been demonstrated, although so far with a classical EM field; e.g., todorov2010 ; kono2019 . Once the quantized field regime is reached, one can imagine coupling several N-qubit systems to a common cavity field, e.g. a multi-quantum well system or excitons in different valleys of 2D semiconductors TMD2015 . The control of specific single-qubit excitations in the monolithic systems is a challenge. If such a control is desirable, it could be easier to deal with a system of spatially isolated qubits in a cavity, e.g. quantum dots or molecular qubits, as proposed in gray2015 .

When the qubits are not identical, this opens more possibilities for quantum state manipulation but also makes the dynamics more complex. The conceptually simplest system includes two qubits which have different transition frequencies. Here we describe their dissipationless dynamics whereas the effects of dissipation are included in Sec. IV.

Assume that both transition frequencies are close enough to the field frequency, so that we can still use the RWA:

$$\hat{H}=\hbar\omega\left(\hat{c}^{\dagger}\hat{c}+\frac{1}{2}\right)+\sum_{j=1,2}\left[W_{j}\hat{\sigma}_{j}^{\dagger}\hat{\sigma}_{j}-\hbar\Omega_{Rj}\left(\hat{\sigma}_{j}^{\dagger}\hat{c}+\hat{\sigma}_{j}\hat{c}^{\dagger}\right)\right],$$

(31)

where $\Omega_{R1,2}=\frac{\mathbf{d}_{1,2}\mathbf{\cdot E}}{\hbar}$. However, the frequency detuning values can be arbitrary with respect to Rabi frequencies. Consider only the ground state and lowest excited state with excitation energy close to single photon energy (i.e. $n=1$). The wave function including only those states is

$$\Psi=C_{000}\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{100}\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{010}\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle+C_{001}\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle.$$

(32)

Here the order of the terms is always $\left|photon\right\rangle\left|atomN1\right\rangle\left|atomN2\right\rangle$.

From the Schrödinger equation we obtain

$$\dot{C}_{000}+i\frac{\omega}{2}C_{000}=0;$$

(33)

$$\dot{C}_{100}+i\omega\frac{3}{2}C_{100}-i\left(\Omega_{R1}^{*}C_{010}+\Omega_{R2}^{*}C_{001}\right)=0,$$

(34)

$$\dot{C}_{010}+i\left(\omega\frac{1}{2}+\frac{W_{1}}{\hbar}\right)C_{010}-i\Omega_{R1}C_{100}=0,$$

(35)

$$\dot{C}_{001}+i\left(\omega\frac{1}{2}+\frac{W_{2}}{\hbar}\right)C_{001}-i\Omega_{R2}C_{100}=0,$$

(36)

One can eliminate high frequencies with $\left(\begin{array}[]{c}C_{100}\\ C_{010}\\ C_{001}\end{array}\right)=\left(\begin{array}[]{c}G_{100}\\ G_{010}\\ G_{001}\end{array}\right)e^{-i\omega\frac{3}{2}t},$ which gives

$$\dot{G}_{100}-i\left(\Omega_{R1}^{*}G_{010}+\Omega_{R2}^{*}G_{001}\right)=0,$$

(37)

$$\dot{G}_{010}+i\Delta_{1}G_{010}-i\Omega_{R1}G_{100}=0,$$

(38)

$$\dot{G}_{001}+i\Delta_{2}G_{001}-i\Omega_{R2}G_{100}=0,$$

(39)

where $\Delta_{1,2}=\frac{W_{1,2}}{\hbar}-\omega.$ The solution for coefficients can be found as $\left(\begin{array}[]{c}G_{100}\\ G_{010}\\ G_{001}\end{array}\right)\propto e^{\Gamma t};$

$$\left(\begin{array}[]{ccc}\Gamma&-i\Omega_{R1}^{*}&-i\Omega_{R2}^{*}\\ -i\Omega_{R1}&\Gamma+i\Delta_{1}&0\\ -i\Omega_{R2}&0&\Gamma+i\Delta_{2}\end{array}\right)\left(\begin{array}[]{c}G_{100}\\ G_{010}\\ G_{001}\end{array}\right)=0,$$

which gives an equation for eigenvalues,

$$\Gamma\left(\Gamma+i\Delta_{1}\right)\left(\Gamma+i\Delta_{2}\right)+|\Omega_{R1}|^{2}\left(\Gamma+i\Delta_{2}\right)+|\Omega_{R2}|^{2}\left(\Gamma+i\Delta_{1}\right)=0$$

(40)

This is a cubic equation which can be solved analytically or on a computer.

Figures 1 and 2 show one example of the solution for eigenenergies and eigenstates, respectively, when $|\Omega_{R1}|=\Omega_{R}$, $|\Omega_{R2}|=2\Omega_{R}$, and $\Delta_{2}=\Delta_{1}+5\Omega_{R}$. Two anticrossing points are visible in the spectra in Fig. 1. As is clear from Fig. 2, by controlling the initial conditions for a given eigenstate and tuning through resonances with a cavity mode one can arrive at an arbitrary entangled state or transfer the excitation from one qubit to another or from any of the qubits to the photon state.

If $\Delta_{1}=\Delta_{2}=\Delta$, we obtain

$$\left(\Gamma+i\Delta\right)\left(\Gamma^{2}+i\Delta\Gamma+|\Omega_{R1}|^{2}+|\Omega_{R2}|^{2}\right)=0.$$

(41)

In addition to the solutions of $\ \Gamma^{2}+i\Delta\Gamma+|\Omega_{R1}|^{2}+|\Omega_{R2}|^{2}=0$ denoted as states $a$ and $b$ in Figs. 3 and 4 below, we have a root $\Gamma=-i\Delta$ for any $\Omega_{R1,2}$. This solution is a dark state which does not couple to the field:

$$\left(\begin{array}[]{c}C_{100}\\ C_{010}\\ C_{001}\end{array}\right)\propto\left(\begin{array}[]{c}0\\ 1\\ -\frac{\Omega_{R1}^{*}}{\Omega_{R2}^{*}}\end{array}\right).$$

(42)

Figures 3 and 4 illustrate this example, showing eigenenergies and eigenstates for $\Delta_{1}=\Delta_{2}=\Delta$. Clearly, one state denoted as Dark state on the figures is decoupled from light at any detunings whereas the two bright states $a$ and $b$ interact strongly through coupling to a cavity mode. By changing the initial conditions one can store the excitation in a dark state and read it out on demand, of course within the limits imposed by decoherence.

As a final example, consider the case $\Delta_{1}=0$, $\Delta_{2}=\Delta$, which gives

$$|\Omega_{R2}|^{2}\Gamma+\left(|\Omega_{R1}|^{2}+\Gamma^{2}\right)\left(\Gamma+i\Delta\right)=0.$$

(43)

In particular, when $\left|\Delta\right|\gg|\Omega_{R1,2}|$ there is an approximate analytic solution $\Gamma_{1}\approx-i\Delta$, $\Gamma_{2,3}\approx\pm|\Omega_{R1}|$ which means that the field couples only to a resonant qubit.