1 Introduction↩︎

In the seminal paper dated back to 1958, John J. Hopfield proposed a model of photons interacting with excitons in the regime of strong coupling, where the system is quantized in terms of hybrid light-matter particles named polaritons [1]. Since then, experiments have confirmed this theory in diverse settings. The polariton family now includes, among others, photons strongly coupled to excitons, phonons, plasmons, polarons, magnons, Rydberg and ultracold atoms [2]–[12], existing in various materials and geometries including cavities, waveguides, nanorods, nanotubes, two-dimensional materials and moiré structures [13]–[20]. Depending on the physical platform, they exhibit unique physical properties such as extremely low effective mass, bosonic condensation, hyperbolic dispersion, ultrastrong coupling, or strong interactions [21]–[25]. Applications of polaritons include ultralow threshold lasing, topological lasers, optical logic, neurmorphic computing and quantum simulations [26]–[36].

Polariton systems are usually well described by the Hopfield Hamiltonian, which in the rotating wave approximation (RWA) reduces to the simple eigenmode problem \[\begin{align} \label{eq:Hopfield}\begin{pmatrix} E_C(\mathbf{k}) & \Omega/2 \\ \Omega/2 & E_X \end{pmatrix}\begin{pmatrix} C \\ X \end{pmatrix} = E(\mathbf{k}) \begin{pmatrix} C \\ X \end{pmatrix} \end{align}\tag{1}\] where \(C,X\) are quantum amplitudes of photons and excitons (or other matter modes) in the polariton eigenstate, \(\mathbf{k}\) is the momentum, \(E_C(\mathbf{k})\) is the photon energy dispersion, \(E_X\) is matter excitation energy (here constant), and \(\Omega\) is the Rabi energy quantifying the light-matter interaction. The resulting upper and lower polariton mode energies are \[\begin{align} \label{eq:Epm}E_\pm(\mathbf{k})&=\frac{E_C(\mathbf{k})+E_X}{2}\pm\frac{1}{2}\sqrt{(E_C(\mathbf{k})-E_X)^2+\Omega^2} \end{align}\tag{2}\] This model is also valid in systems that are partly confined, such as microcavities, where the momentum \(\mathbf{k}\) is replaced by the momentum parallel to the cavity plane \(\mathbf{k}_\parallel\). While the predictions of this model are often very accurate, the parameters are almost always obtained by fitting to the experimental data. On the other hand, one can ask when this simple model is a correct description of physical systems. The most obvious case which goes beyond its applicability is when more than two modes of light or matter interact, but in this case a simple multimode extension of 1 usually suffices. When translational variance is absent, \(\mathbf{k}\) is not a good quantum number, but as long as spatial modes can be well defined using another quantum number, Eq. 1 can be still used. It is also possible to include polarization degrees of freedom [37], effects of dissipation and pumping [38] and ultrastrong light-matter coupling (beyond RWA) [1], [39] by suitable extensions of the Hopfield model.

Nevertheless, there are situations where descriptions based on Eq. 1 and its generalizations fail. In the multimode case, when light-matter interaction energy is stronger than the bare energy level spacing, it modifies mode shapes so much that polariton modes can no longer be approximated as superpositions of a few bare modes. In the extreme case, an infinite number of bare modes may be needed to reconstruct the polariton mode - for example when polariton modes are localized while bare modes are plane waves. An example of such a system is a planar dielectric microcavity that contains a finite volume of an active medium where excitons can form (see Fig. 1). Other examples include non-orthogonal modes in nanophotonic structures and excitons bound by the interaction with light [40], [41].

A question arises whether it is possible to determine an accurate model of polaritons valid in an arbitrary physical structure that would be applicable in both classical and quantum regimes. Ideally, it would provide the form of a Hamiltonian (in the conservative case) or a quantum master equation (in the dissipative case), with coefficients resulting from material properties, without any fitting parameters. Moreover, the model should be efficient for analytical and numerical computations, meaning that the basis of modes should be optimal for the problem at hand. This issue is particularly important in the quantum regime where the sub-optimal choice of basis could lead to a massive increase of computational complexity.

Here, we propose a recipe for constructing a quantum mechanical model of polaritons in an arbitrary dielectric environment. We determine eigenmodes of a light-matter system starting from the Lagrangian density including light-matter interaction. In the conservative case, we employ Bogoliubov transformation to obtain the diagonal form of the Hamiltonian, with eigenmodes that are in correspondence to the solutions of classical Maxwell equations coupled to matter polarization. Our method supplements the treatments based on the approach of Huttner and Barnett [42] by simplifying calculations and providing a straightforward quantum-classical correspondence. Next, we consider dissipative systems in the linear regime. While many theoretical approaches have been proposed in the past to describe photons interacting with the environment [42]–[54], none of them demonstrated a straightforward and general method to obtain a quantum master equation with no fitting parameters, valid in the strong coupling regime. In this work, by means of the technique of third quantization [55], we find the optimal basis that corresponds to normal superoperator modes of the dissipative system with linear losses. In analogy to the conservative case, we find that the solutions of Maxwell equations including optical and material losses, as well as losses resulting from boundary conditions, correspond to the normal modes of the quantum Liouvillian. We show how this allows one to determine the exact form of the master equation in an optimal, diagonal basis.

As an application of our theory, we consider a nonlinear system of interacting polaritons. We show how our theory can be used to engineer polariton many-body Hamiltonians. For instance, constricting the volume of emitters leads to the increase of two-body interactions between polaritons even in the absence of increased optical confinement. This provides a new strategy in the quest for strongly interacting exciton-polaritons, which can boost existing methods that range from optical confinement [56]–[58], dipolar interactions [15], [59], [60] to polaron-mediated interactions [9], [61]. We also consider interaction between two spatially separated polariton modes and predict strong nonlocal interactions, which can be used to induce nonlocal correlations and mode entanglement of emitted light.

2 Conservative system↩︎

The problem of interaction of light and matter at the quantum level has been considered in great detail in the past [42]–[54], including the seminal work of Huttner and Barnett [42], which also treated interaction with matter reservoir modes. In this section, we propose an alternative approach to solving the diagonalization problem, which simplifies the analysis and emphasizes a physically important aspect, the quantum-classical correspondence [62]. Additionally, it provides a basis for the consideration of the dissipative case in the next section, ultimately leading to the formulation of the quantum master equation.

2.1 Classical regime↩︎

We consider coupling of the electromagnetic field and matter polarization under the assumption of the small size of emitters compared to the wavelength of light. Lagrangian density of electromagnetic field in a dielectric medium coupled to localized matter modes can be written in the Coulomb gauge as \(\mathcal{L}=\mathcal{L}_{E M}+\mathcal{L}_X+\mathcal{L}_{EM-X}\), with [1], [42], [63] \[\begin{align} & \mathcal{L}_{E M}=\frac{\epsilon_b({\boldsymbol{r}})}{2} \dot{\mathbf{A}}^2-\frac{1}{2 \mu_0}(\nabla \times \mathbf{A})^2 \\ & \mathcal{L}_X=\frac{1}{2} \rho \dot{\mathbf{X}}^2-\frac{1}{2} \rho \omega_0^2 \mathbf{X}^2 \\ & \mathcal{L}_{E M-X}=-\alpha({\boldsymbol{r}}) \mathbf{A} \cdot \dot{\mathbf{X}} \end{align}\] where \(\epsilon_b({\boldsymbol{r}})\) is the space-dependent background dielectric constant, \(\rho\) is the mass density, \(\omega_0\) is matter oscillator frequency, \(\alpha({\boldsymbol{r}})\) is the light-matter coupling coefficient, \(\mathbf{E}=-\dot{\mathbf{A}}\), and \(\mathbf{B}=\nabla \times \mathbf{A}\). Polarization density can be identified as \(\mathbf{P}=\alpha({\boldsymbol{r}})\mathbf{X}\). The electromagnetic part of the Lagrangian corresponds to Maxwell equations in a dielectric. The matter part is described as harmonic oscillators. For simplicity, we assume that frequency and mass density of the matter component are space-independent, but these assumptions are not crucial for the following considerations. The background refractive index \(\epsilon_b({\boldsymbol{r}})\) is here frequency independent, but later we will show how to incorporate dispersion of the background material. The above model is in some sense minimal, but can already precisely describe a broad range of inhomogeneous polariton systems. We will also see that generalizations to incorporate e.g. anisotropy, magnetic response or electric currents are easy to implement.

While the above is written as a Lagrangian function in terms of continuous fields, we will actually have in mind a version of the model where a finite number of electromagnetic field modes and matter modes are considered. Physically, this corresponds to considering finite system volume \(V\), with an appropriate UV cutoff for the electromagnetic field. For matter modes, the cutoff corresponds to the size of the emitter, such as an exciton or an atom. On the other hand, we assume that the spatial variation of the electromagnetic field component occurs on much larger spatial scales than the size of emitters, which justifies the use of the continuous field to describe matter polarization. Despite such discretization of modes, we will use continuous integrals over volume \(V\) for convenience, keeping in mind that they should be replaced by corresponding summations in the case of matter fields. We will not use continuity of these fields anywhere in the derivation. The total number of modes of the model will be denoted by \(2N=2N_{\rm el}+2N_{\rm mat}\), where \(2N_{\rm el}\) and \(2N_{\rm mat}\) are the number of bare electromagnetic and matter modes, respectively, and the coefficient 2 stems from time-reversal symmetry. This discretization avoids certain computational issues occurring in the infinite-mode case.

The Lagrangian \(L=\int_V\mathcal{L}\,d{\boldsymbol{r}}\) corresponds to Maxwell equations coupled to polarization field. Their \(2N\) solutions \((\mathbf{E}_j({\boldsymbol{r}})\), \(\mathbf{X}_j({\boldsymbol{r}}))\) with \(j=1\dots 2N\) obey \[\begin{align} \nabla \times \nabla \times \mathbf{E}_j=&\frac{\epsilon_b({\boldsymbol{r}})\omega_j^2}{\epsilon_0 c^2} \mathbf{E}_j+ \frac{\alpha({\boldsymbol{r}})\omega_j^2}{\epsilon_0 c^2} \mathbf{X}_j \tag{3}\\ (\omega_0^2 - \omega_j^2) \tag{4} \mathbf{X}_j=&\frac{\alpha({\boldsymbol{r}})}{\rho} \mathbf{E}_j \end{align}\]

where \(\omega_j\) is the oscillation frequency. We can move to the complex representation of the electromagnetic field where negative frequency eigenmodes are complex conjugates of their positive frequency partners. For each mode with frequency \(\omega_j\), there exists a mode with frequency \(-\omega_j\) and conjugated \(E_j^*\) profile (i.e. the magnetic field profile has a flipped sign).

The above equations can be recast as an eigenvalue problem for a Hermitian operator \(\hat{O}\) on the composite six-dimensional vector \(\mathbf{F} =\left(\mathbf{G}({\boldsymbol{r}}),\mathbf{Y}({\boldsymbol{r}})\right)^T\) with \(\mathbf{G}({\boldsymbol{r}})=\sqrt{\epsilon_b({\boldsymbol{r}})} \mathbf{E}({\boldsymbol{r}})\) and \(\mathbf{Y}({\boldsymbol{r}})=\rho^{1/2} \omega_0\mathbf{X}({\boldsymbol{r}})\). Its eigenmodes form an orthogonal set \[\begin{align} \label{eq:normalization} \int_V {\boldsymbol{F}}^*_i({\boldsymbol{r}}) {\boldsymbol{F}}_j({\boldsymbol{r}})d {\boldsymbol{r}}&= \int_V \left({\boldsymbol{G}}^*_i({\boldsymbol{r}}) {\boldsymbol{G}}_j({\boldsymbol{r}}) + {\boldsymbol{Y}}^*_i({\boldsymbol{r}}) {\boldsymbol{Y}}_j({\boldsymbol{r}}) \right)d {\boldsymbol{r}}=\nonumber\\&= \mathcal{N}_j\delta_{ij} \end{align}\tag{5}\] of solutions with real eigenvalues \(\omega_j^2\), complete in a certain subspace. In particular, in the case \(\epsilon_b({\boldsymbol{r}})=1\), \(\alpha({\boldsymbol{r}})=0\) the subspace contains linear combinations of all matter modes and transversely polarized electromagnetic plane waves [63]. Here \(\int_V d {\boldsymbol{r}}|{\boldsymbol{G}}_j({\boldsymbol{r}})|^2\) and \(\int_V d {\boldsymbol{r}}|{\boldsymbol{Y}}_j({\boldsymbol{r}})|^2\) quantify the photonic and matter components of the polariton, in analogy to the squared Hopfield coeficients \(|X|^2\) and \(|C|^2\) in Eq. 1 . Indeed, it can be verified that Eqs. 3 and 4 transform to the Hopfield problem Eq. 1 in the homogeneous case and under RWA \(|\omega_j-\omega_0|\ll\omega_0\), provided that \(\mathbf{G}=C\hat{e}\,e^{i\mathbf{k}{\boldsymbol{r}}}\), \(\mathbf{Y}=X\hat{e}\,e^{i\mathbf{k}{\boldsymbol{r}}}\), \(\Omega=-\hbar\alpha/\sqrt{\rho\varepsilon}\), \(E_X=\hbar\omega_0\) and \(E_C(k)=\hbar |\mathbf{k}|c\sqrt{\varepsilon_0/\varepsilon}\), where \(\hat{e}\) is the polarization unit vector perpendicular to momentum \(\mathbf{k}\).

The normalization coefficients \(\mathcal{N}_j\) can be chosen to be unity, leading to an orthonormal set of eigenmodes. However, for reasons that will become clear later, we choose \(\mathcal{N}_j\) in such a way that the energy of each mode (value of the classical Hamiltonian functional \(H\)) is equal to \(\hbar |\omega_j|\). Hermiticity of operator \(\hat{O}\) represented by Eqs. 3 and 4 in the general case \(\langle {\boldsymbol{F}}_1|\hat{O} {\boldsymbol{F}}_1\rangle=\langle \hat{O} {\boldsymbol{F}}_1|{\boldsymbol{F}}_2\rangle\) can be checked directly by integration by parts. We will assume that all of these classical solutions are stable, i.e. \(\omega_j^2>0\) for all \(j\), which is a physically justified requirement of the stability of vacuum, as energy cannot appear from nowhere. Note that this requirement is not fulfilled for all parameters, for example it is not the case when matter oscillators have a negative mass density \(\rho\).

In the absence of coupling, \(\alpha=0\), the first equation Eq. 3 describes electromagnetic field in a dielectric medium. In the presence of coupling, the solutions of the coupled equations 3 and 4 are polariton solutions in the classical regime. The corresponding Hamiltonian density \(\mathcal{H}=\mathcal{H}_{E M}+\mathcal{H}_X+\mathcal{H}_{EM-X}\) and Hamiltonian \(H=\int_V \mathcal{H}\,d {\boldsymbol{r}}\) are obtained by introducing canonical momenta \(\Pi_{A,i}=\frac{\partial \mathcal{L}}{\partial \dot{A}_i}=\epsilon_b(\mathbf{r}) \dot{A}_i\) (i.e. \(\mathbf{\Pi}_A=-\mathbf{D}\)) and \({\Pi}_{X,i}=\frac{\partial\mathcal{L}}{\partial \dot{{X_i}}}=\alpha {A_i}+\rho \dot{{X_i}}\,\)

\[\begin{align} \mathcal{H}_{EM} =&\frac{\mathbf{\Pi}_A^2}{2\epsilon_b(\mathbf{r})}+\frac{1}{2\mu_0}(\mathbf{\nabla}\times \mathbf{A})^2\,,\\ \mathcal{H}_{X}=&\frac{1}{2\rho}\mathbf{\Pi}_X^2+\frac{1}{2} \rho \omega_0^2 \mathbf{X}^2\,,\tag{6}\\ \mathcal{H}_{EM-X}=&-\frac{\alpha}{\rho}\mathbf{\Pi}_X \mathbf{A} + \frac{\alpha^2}{2\rho}\mathbf{A}^2\,,\tag{7} \end{align}\] where \(\mathcal{H}_{EM-X}\) includes the term \(\frac{\alpha^2}{2\rho}\mathbf{A}^2\), which, in the quantum regime, corresponds to interaction of light with vacuum fluctuations of the matter field.

2.2 Quantization↩︎

Canonical quantization of light and matter fields is performed by demanding appropriate commutation relations between variables and their conjugates. First, we determine the diagonal form of the above Hamiltonian. Instead of doing this explicitly, we will make use of known results. In the absence of light-matter coupling, \(\alpha({\boldsymbol{r}})=0\), it is well known that both the electromagnetic field described by \(\mathcal{H}_{EM}\) and the matter field described by \(\mathcal{H}_{X}\) can be quantized as independent sets of harmonic oscillators. Assuming space of finite volume, the resulting diagonal form is \(\hat{H}_{EM}=\sum_{+i} \hbar \omega_i \left(\hat{a}_i^\dagger \hat{a}_i+\frac{1}{2}\right)\) [63], where the "\(+i\)" summation is over positive frequencies only, and in the Heisenberg picture the corresponding electromagnetic field and vector potential are \[\begin{align} \hat{\boldsymbol{E}}({\boldsymbol{r}},t)&=\sum_{+i} \hat{a}_i {\boldsymbol{E}}_i^{(\alpha=0)}({\boldsymbol{r}}) + {\rm h.c.}\\ \nonumber \hat{\boldsymbol{A}}({\boldsymbol{r}},t)&=-i\sum_{+i}\frac{1}{\omega_i} \hat{a}_i {\boldsymbol{E}}_i^{(\alpha=0)}({\boldsymbol{r}}) + {\rm h.c.} \end{align}\] where \(E_i^{(\alpha=0)}({\boldsymbol{r}})\) are eigensolutions to uncoupled Eq. (3 ). The difference with respect to [63] results from different choice of normalization coefficients \(\mathcal{N}_i\). On the other hand, \(\hat{H}_{X}=\sum_\nu\int_V dr \hbar \omega_0 \left(\hat{b}_\nu^\dagger(r) \hat{b}_\nu(r)+\frac{1}{2}\right)\) where \(\nu=x,y,z\), while \(\hat{a}_i\), \(\hat{b}_\nu({\boldsymbol{r}})\) are annihilation operators obeying bosonic commutation relations. Matter operator \[\hat{\boldsymbol{X}}({\boldsymbol{r}},t)=\left(\frac{\hbar}{2\rho\omega_0}\right)^{1/2} \sum_\nu \hat{b}_\nu({\boldsymbol{r}}) \vec{e}_\nu + {\rm h.c.}\] is simply coordinate operator of a 3D harmonic oscillator at each (discrete) point in space where emitter mode exists. Since the total Hamiltonian is quadratic, conjugate momenta \({\boldsymbol{\Pi}}\) and \({\boldsymbol{\Pi}}_X\) can be written as linear combinations of operators \(\hat{a}_i\), \(\hat{a}_i^\dagger\), and \(\hat{b}_\nu(r)\), \(\hat{b}^\dagger_\nu(r)\), respectively.

2.3 Diagonalization in polariton modes↩︎

Considering now the coupled case, \(\alpha({\boldsymbol{r}})\neq0\), it is clear that the coupling Hamiltonian density \(\mathcal{H}_{EM-X}\), Eq. (7 ), after substitution of the above formulas, becomes quadratic in creation and annihilation modes of the uncoupled system. Note that adding the coupling not only introduces a coupling term in the Hamiltonian but also alters the conjugate momentum corresponding to a particular coordinate, which is taken into account by the \({\boldsymbol{A}}^2\) term in Eq. (7 ). This term, however, is also quadratic. The complete Hamiltonian has the general form \(\hat{H}=\frac{1}{2}\sum_{ij}\hat{c}^\dagger_i \mathcal{A}_{ij} \hat{c}_j + \hat{c}^\dagger_i \mathcal{B}_{ij}\hat{c}^\dagger_j + {\rm h.c.}\), with \(\mathcal{A},\mathcal{B}\in \mathbb{C}^{N\times N}\) where the set of operators \(\hat{c}_i\) is composed of both \(\hat{a}_i\) and \(\hat{b}_\nu(r)\). Recall that \(\hat{b}_\nu(r)\) includes a discrete set of emitter positions so this set is also finite. This Hamiltonian can be diagonalized using Bogoliubov transformation \[\hat{c}_i = \sum_{+j} \left( u_{ij} \hat{P}_j + v_{ij}^* \hat{P}_j^\dagger \right)\] into \[\label{eq:Hpol} \hat{H}=\sum_{+j} \hbar \omega_j \left(\hat{P}_j^\dagger \hat{P}_j+\frac{1}{2}\right)\tag{8}\] where \(\hat{P}_j\) are operators of polariton modes, which obey bosonic commutation relations due to symplecticity of the transformation. While the existence of such transformation is not guaranteed for bosons, it always exists when the eigenmodes of the associated Bogoliubov-de Gennes matrix \(\hat{H}_{\text{BdG}}\) defined by \[\begin{align} \hat{H}+\mathrm{const}&= \frac{1}{2}\left(c^\dagger\,c\right) \begin{pmatrix} \mathcal{A} & \mathcal{B} \\ \mathcal{B}^\dagger & \mathcal{A}^T \end{pmatrix} \begin{pmatrix} c \\ c^\dagger \end{pmatrix} = \nonumber\\&= \frac{1}{2}\left(c^\dagger\, c\right)\hat{H}_{\text{BdG}} \begin{pmatrix} c \\ c^\dagger \end{pmatrix} \end{align}\] are stable, which is equivalent to the positivity of eigenvalues or stability of the corresponding classical system [64]. This condition is fulfilled whenever there is a well defined ground state of the system, i.e. \(\hat{H}_{\text{BdG}}\) is bounded from below. Due to the conservation of the number of modes under the transformation, the number of polariton modes in Eq. 8 is equal to the sum of bare electromagnetic and matter bosonic modes (equal to N since each bosonic mode corresponds to a pair of positive and negative frequency classical modes).

2.4 Quantum-classical correspondence↩︎

While the existence of Bogoliubov transformation does not give information about the explicit form of \(\hat{P}_j\) operators, or coefficients \(u_{ij},v^*_{ij}\), physical observable operators can be determined by consideration of the classical limit, i.e. a coherent state corresponding to \(j\)-th polariton eigenmode \(\hat{p}_j\), which we denote by \(|\alpha_j\rangle\), where \(\alpha_j\) is a complex number. Expectation values of electric field and matter coordinates have to obey Eqs. (3 ) and (4 ). Since only quadratic terms are present in the Hamiltonian, this suggests that they are linear combinations of polariton operators of the form \[\begin{align} \hat{\boldsymbol{E}}({\boldsymbol{r}},t)&=\sum_{+j} \hat{P}_j {\boldsymbol{E}}_j({\boldsymbol{r}}) + {\rm h.c.} \tag{9}\\ \hat{\boldsymbol{X}}({\boldsymbol{r}},t)&=\sum_{+j} \hat{P}_j {\boldsymbol{X}}_j({\boldsymbol{r}}) + {\rm h.c.} \tag{10} \end{align}\] Since the Hamiltonian is diagonal in the basis \(\hat{P}_j\), the first terms on the right hand sides are positive frequency parts \(\hat{\boldsymbol{E}}^{(+)}\) and \(\hat{\boldsymbol{X}}^{(+)}\), and their Hermitian conjugates are negative frequency parts \(\hat{\boldsymbol{E}}^{(-)}\) and \(\hat{\boldsymbol{X}}^{(-)}\). Normally ordered terms, where positive frequency operators are on the right of negative frequency operators, are to be used to determine physical observables [63], [65]. It is straightforward to check that formulas of Eqs. (9 ) and (10 ) give the correct values of correlation functions in the classical limit. Consider the two-point first-order spatial correlation function \[\begin{align} \label{eq:E1E2} G^{(1)}_{\mu,\nu}&({\boldsymbol{r}}_1,{\boldsymbol{r}}_2;t,t) = \langle \alpha_j | { E}_\nu^{(-)}({\boldsymbol{r}}_1,t) { E}_\mu^{(+)}({\boldsymbol{r}}_2,t)| \alpha_j\rangle =\\ \nonumber&= |\alpha_j|^2 { E}_{j,\nu}^*({\boldsymbol{r}}_1){ E}_{j,\mu}({\boldsymbol{r}}_2) = \langle n_j \rangle { E}_{j,\nu}^*({\boldsymbol{r}}_1){ E}_{j,\mu}({\boldsymbol{r}}_2) \end{align}\tag{11}\] where \(\mu,\nu=x,y,z\), and we used that for a coherent state \(\hat{P}_j|\alpha_j\rangle = \alpha_j |\alpha_j\rangle\). This is exactly the correlation function that one would expect for a coherent state in the classical limit, in the state given by \({\boldsymbol{E}}_j\) and \({\boldsymbol{X}}_j\), considering that \({\boldsymbol{E}}_j\) is the electric field corresponding to a coherent state with one polariton on average, \(\langle n_j \rangle=1\). Recall that \({\boldsymbol{E}}_j\) and \({\boldsymbol{X}}_j\) are normalized in Eq. 5 such that the expectation value of the Hamiltonian is equal to \(\hbar |\omega_j|\), and any real-valued classical electromagnetic field necessarily appears as an equal contribution of positive- and corresponding negative-frequency modes.

Analogously, correlations between \({\boldsymbol{E}}\) and \({\boldsymbol{X}}\) fields can be obtained. It is also clear that no linear combinations of polariton operators other those given by Eqs. (9 ) and (10 ) can result in correct first-order correlation functions of all kinds, while nonlinear terms are not present in the Bogoliubov transformation.

The above relations can be inverted by multiplication of both sides (9 ) by \(\epsilon_b({\boldsymbol{r}}) {\boldsymbol{E}}_j^*({\boldsymbol{r}})\) and both sides of (10 ) by \(\rho \omega_0^2 {\boldsymbol{X}}_j^*({\boldsymbol{r}})\), adding to each other and integrating over space to yield the explicit form of the polariton operator \[\hat{P}_j = \mathcal{N}_j^{-1}\int_V d{\boldsymbol{r}} \left( \hat{\boldsymbol{E}} ({\boldsymbol{r}}) \epsilon_b({\boldsymbol{r}}){\boldsymbol{E}}_j^*({\boldsymbol{r}}) + \hat{\boldsymbol{X}} ({\boldsymbol{r}}) \rho \omega_0^2 {\boldsymbol{X}}_j^*({\boldsymbol{r}})\right)\] where we have used orthogonality of \((\sqrt{\epsilon_b({\boldsymbol{r}})}{\boldsymbol{E}}({\boldsymbol{r}}),\rho^{1/2} \omega_0 {\boldsymbol{X}}({\boldsymbol{r}}))\). The polariton annihilation operator annihilates both photons and excitons according to the amplitudes \(\epsilon_b({\boldsymbol{r}}){\boldsymbol{E}}_j^*({\boldsymbol{r}})\) and \(\rho \omega_0^2 {\boldsymbol{X}}_j^*({\boldsymbol{r}})\), respectively.

In summary, an arbitrary linear light-matter system in the conservative case can be diagonalized using Bogoliubov transformation. Polariton eigenmodes are given explicitly in the electromagnetic field and matter polarization operators, Eqs. 9 and 10 , where the spatial shapes of the modes correspond to classical solutions of Maxwell equations coupled to the matter field. We note that this result has been formulated in previous works, in particular in [42] in the homogeneous case, and in [54] in the inhomogeneous case, which also included the effects of frquency-dependent dielectric constant. In our approach, explicit determination of the transformation and the corresponding coefficients is not required, since all the physically relevant operators can be obtained by consideration of the classical limit. Therefore, the calcuations become much more concise. Importantly, as we demonstrate in the next section, this approach can be generalized to the dissipative case, which allows one to determine the form of the quantum master equation.

3 Dissipative system↩︎

In the dissipative case, we consider a quadratic system as described above, but coupled to reservoir modes, which results in loss/gain channels to exciton, photon, polariton, or other modes. In the Born-Markov approximation, reduced density matrix of the system follows Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) master equation \[\label{eq:Lindblad} \dot{\hat{\rho}} = \widehat{\mathcal{L}}\hat{\rho} = \frac{1}{i\hbar} [\hat{H},\hat{\rho}] +\sum_\mu \gamma_\mu \left(\hat{L}_\mu \hat{\rho} \hat{L}_\mu^\dagger -\frac{1}{2}\left\{ \hat{L}_\mu^\dagger \hat{L}_\mu,\hat{\rho}\right\}\right),\tag{12}\] where \(\widehat{\mathcal{L}}\) is the Liouvillian superoperator (we will indicate superoperators using "wide hats"), \(\hat{H}\) is the Hamiltonian part operator given by Eq. 8 , and \(\mu\) enumerates different gain/loss channels with rates \(\gamma_\mu\). In the polariton basis introduced in the previous section, these channels are defined as \[\hat{L}_\mu = \sum_j \left( l_{\mu,j} \hat{P}_j+g_{\mu,j} \hat{P}^\dagger_j\right)\] with arbitrary complex coefficients \(l_{\mu,j}\) and \(g_{\mu,j}\). We consider only linear Lindbladian terms to preserve the quadratic form of the Liouvillian, and assume that the system has a stable nonequilibrium steady state. In particular, in the absence of gain \(g_{\mu,j}=0\), the steady state is the polariton vacuum which fulfills \(P_j |0\rangle_P=0\). However, we will consider the general case, without assuming that the steady state is a pure state.

3.1 Third quantization treatment↩︎

In general, conservative and dissipative parts of 12 do not have the same diagonal basis of eigenmodes. To find a common basis and diagonalize the complete Liouvillian, one should perform a transformation that is an analog of Bogoliubov transformation, but in the dissipative case. Such a transformation for quadratic bosonic systems with linear Lindbladians has been proposed under the name of third quantization [55]. It turns out that while it is generally not possible to find a common basis in the operator space, it is possible to find an adequate superoperator basis. Note that third quantization can be implemented in several different bases [66], [67], each of them having certain advantages, but here we follow the original approach of [55], which will allow to use very similar arguments as in Sec. 2.

In third quantization, the central objects are superoperators, such as \(\widehat{\mathcal{L}}\), acting on operators such as the density matrix operator denoted as \(|\hat{\rho}\rangle\!\rangle\). One can think of \(|\hat{\rho}\rangle\!\rangle\) as a vectorized representation of the density matrix and superoperators as operators acting on these vectors. We are interested in expectation values of observables denoted by \((\!(\hat{A}|\hat{\rho}\rangle\!\rangle:=\mathrm{Tr}\,\hat{A}\hat{\rho}\), where A is an operator representing the observable, a member of dual space with respect to the space of density operators, in the sense that \(\hat{A}\hat{\rho}\) products are trace-class [55]. Basic superoperators are simply annihilation \(\hat{a}_j\) or creation \(\hat{a}^\dagger_j\) operators acting from the left or from the right on the density matrix. We define superoperators \(\widehat{a}^{\dagger L}|\hat{\rho}\rangle\!\rangle = |\hat{a}^{\dagger}\hat{\rho}\rangle\!\rangle\) and \(\widehat{a}^{\dagger R}|\hat{\rho}\rangle\!\rangle = |\hat{\rho} \hat{a}^{\dagger}\rangle\!\rangle\), while \((\!(\hat{A}|\widehat{a}^{\dagger L}= (\!(\hat{A}\hat{a}^{\dagger}|\) and \((\!(\hat{A}|\widehat{a}^{\dagger R}= (\!(\hat{a}^{\dagger}\hat{A}|\) etc. For each bosonic mode, one may define four "basic" superoperators \[\begin{align} \label{eq:a01} \widehat{a}_{0,j}&=\widehat{a}_j^L, & \widehat{a}'_{0,j}&=\widehat{a}^{\dagger L}_j-\widehat{a}^{\dagger R}_j, \\ \nonumber \widehat{a}_{1,j}&=\widehat{a}^{\dagger R}_j, & \widehat{a}'_{1,j}&=\widehat{a}^{R}_j-\widehat{a}^{L}_j, \end{align}\tag{13}\] where \(j=1..N\) is the bosonic mode index. It is straightforward to check that these superoperators satisfy commutation relations \([\widehat{a}_{\mu,j},\widehat{a}'_{\nu,k}]=\delta_{\mu,\nu}\delta_{j,k}\), \([\widehat{a}_{\mu,j},\widehat{a}_{\nu,k}]=[\widehat{a}'_{\mu,j},\widehat{a}'_{\nu,k}]=0\) with \(\mu,\nu\in\{0,1\}\). One also has \(\widehat{a}_{\nu,j} |0\rangle\!\rangle=0\) where \(|0\rangle\!\rangle\) is the density matrix of vacuum, and, using the cyclic property of trace for trace-class products, \((\!(1|\widehat{a}'_{\nu,j}=0\), where \((\!(1|\) is the identity operator. This allows to define dual-Fock bases of density operators and observable operators \[\label{eq:Fock} |\underline{m}\rangle\!\rangle = \prod_{\nu,j} \frac{(\widehat{a}'_{\nu,j})^{m_{\nu,j}}}{\sqrt{m_{\nu,j}!}}\,|0\rangle\!\rangle, \qquad (\!( \underline{m}| = (\!( 1|\;\prod_{\nu,j} \frac{(\widehat{a}_{\nu,j})^{m_{\nu,j}}}{\sqrt{m_{\nu,j}!}}.\tag{14}\] where the \(2N\)-component index \(\underline{m}=(m_{\nu,j};\nu\in\{0,1\},j\in\{1..N\})^T\) and Fock states fulfill the biorthonormality condition \((\!(\underline{m}'|\underline{m}\rangle\!\rangle=\delta_{\underline{m}', \underline{m}}\). Note that elements of the dual-Fock basis are generally not density matrices nor physical observables, as they are not necessarily Hermitian or unit-trace. However, carefully chosen superpositions of dual-Fock basis vectors can be used to construct physical operators.

In [55], it was shown that quadratic Liouvillians of the form 12 generating stable and non-degenerate dynamics can be diagonalized as \[\label{eq:Lzeta} \widehat{\mathcal{L}}= -\sum_{r=1}^{2N} i \omega_r \widehat{\zeta}'_r \widehat{\zeta}_r\tag{15}\] where \(\omega_r\) can be interpreted as the complex mode frequency, while \(\widehat{\zeta}_r\) and \(\widehat{\zeta}'_r\) are normal modes superoperators given by a symplectic transformation on the extended Nambu vector of superoperators \(\widehat{b}=(\widehat{a}_{0,j},\widehat{a}_{1,j},\widehat{a}'_{0,j},\widehat{a}'_{1,j})^T\) \[\label{eq:Vb}(\widehat{\zeta}_r,\widehat{\zeta}'_r)^T =\mathbf{V} \widehat{b}\tag{16}\] where \(\mathbf{V}\) is a complex \(4N\times4N\) matrix. In general, \(\widehat{\zeta}_r\) and \(\widehat{\zeta}'_r\) are not Hermitian conjugate to each other, but due to symplecticity of the transformation, they satisfy commutation relations \([\widehat{\zeta}_r,\widehat{\zeta}'_s]=\delta_{r,s}\), \([\widehat{\zeta}_r,\widehat{\zeta}_s]=[\widehat{\zeta}'_r,\widehat{\zeta}'_s]=0\) ¹. Moreover, \((\!(1|\widehat{\zeta}'_r=0\), which ensures the conservation of trace in evolution under Eq. 15 . The form of 15 guarantees that in the Fock space defined in Eq. 14 , a stable system has a unique nonequilibrium steady state, the vacuum of \(\widehat{\zeta}_r\) superoperators, which we denote by \(|0_\zeta\rangle\!\rangle:=|{\rm NESS}\rangle\!\rangle\). Superoperators \(\widehat{\zeta}_r\) and \(\widehat{\zeta}'_r\) can be used, in analogy to \(\widehat{a}_{\mu,\nu}\) and \(\widehat{a}'_{\mu,\nu}\) in Eq. 14 , to construct dual-Fock bases of density operators and observable operators by acting on \(|{0_\zeta}\rangle\!\rangle\) and \((\!(1|\), respectively \[\label{eq:FockZeta} |\underline{m}\rangle\!\rangle_\zeta =\prod_{r}\frac{\bigl(\widehat{\zeta}'_{r}\bigr)^{m_{r}}}{\sqrt{m_{r}!}}\, |0_\zeta\rangle\!\rangle, \quad (\!( \underline{m}| _\zeta = (\!( 1|\;\prod_{r} \frac{(\widehat{\zeta}_{r})^{m_{r}}}{\sqrt{m_{r}!}},\tag{17}\] where \(\underline{m}=(m_{r};r\in\{1..2N\})^T\).

3.2 Coherent states of dissipative modes↩︎

Our approach is based on the quantum-classical correspondence. To make use of it, we define analogs of coherent states in the space of superoperators. Let us define \(\widehat{\zeta}^{+}_r\), \(\widehat{\zeta}'^{+}_r\) as superoperators that result in Hermitian conjugates with respect to \(\widehat{\zeta}_r\), \(\widehat{\zeta}'_r\), i.e. \(|\widehat{\zeta}_r^{+}\hat{\rho}^\dagger\rangle\!\rangle=|(\widehat{\zeta}_r \hat{\rho})^\dagger\rangle\!\rangle\) and \(|\widehat{\zeta}_r'^{+}\hat{\rho}^\dagger\rangle\!\rangle=|(\widehat{\zeta}_r' \hat{\rho})^\dagger\rangle\!\rangle\) (for all operators \(\hat{\rho}\) from the Fock space defined via Eq. 17 , even if \(\hat{\rho}\neq\hat{\rho}^\dagger\)) ². It is easy to see from 13 that these operators can be constructed by replacing \(\widehat{a}_{\nu,j}\) with \(\widehat{a}_{1-\nu,j}\). Taking Hermitian conjugate of Eq. 15 , it follows from the conservation of Hermiticity during evolution that for every normal mode superoperator \(\widehat{\zeta}_r\) in 15 with complex frequency \(\omega_r\), \(\widehat{\zeta}^{+}_r\) is a normal mode superoperator with frequency \(-\omega^*_r\). Therefore, it is either another element of the set of normal modes in Eq. 15 , or a linear combination of degenerate normal modes with the same frequency \(-\omega^*_r\) (our assumptions exclude \(\omega_r=-\omega^*_r\)). In the latter case, \(\widehat{\zeta}^{+}_r\) can be incorporated into the normal mode basis through symplectic Gram-Schmidt procedure. Hence, operators \(\widehat{\zeta}^{+}_r\), \(\widehat{\zeta}'^{+}_r\) fulfill the commutation relations just as the other members of the set of normal modes, in particular \([\widehat{\zeta}'^{+}_r,\widehat{\zeta}^{+}_s]=\delta_{r,s}\), \([\widehat{\zeta}'^{+}_r,\widehat{\zeta}_s]=[\widehat{\zeta}^{+}_r,\widehat{\zeta}_s]=[\widehat{\zeta}'^{+}_r,\widehat{\zeta}'_s]=0\). The Liouvillian can be rewritten as \[\label{eq:Lzeta95pp} \widehat{\mathcal{L}}= \sum_{+r}\left( -i \omega_r \widehat{\zeta}'_r \widehat{\zeta}_r + i \omega_r^* \widehat{\zeta}'^+_r \widehat{\zeta}^+_r \right)\tag{18}\] where the summation is over \(N\) normal modes with \(\Re(\omega_r)>0\).

We define the coherent state corresponding to \(r\)-th normal mode as \[\label{eq:alpha} |\alpha\rangle\!\rangle_r = {\rm exp}\left(\alpha\widehat{\zeta}'_r+\alpha^*\widehat{\zeta}'^{+}_r\right)|0_\zeta\rangle\!\rangle,\tag{19}\] where \(\alpha\) is a complex number. These operators are manifestly Hermitian, unlike, for instance, the operator \({\rm e}^{\alpha\widehat{\zeta}'_r -\alpha^*\widehat{\zeta}_r}|0_\zeta\rangle\!\rangle\). Since \((\!(1|\widehat{\zeta}'_r=(\!(1|\widehat{\zeta}'^+_r=0\), we have \(\mathrm{Tr}|\alpha\rangle\!\rangle_r = (\!(1|\alpha\rangle\!\rangle_r=1\). Note that these coherent states differ from these defined in [66], [67], which correspond to the original basis of \(\widehat{a}_{\mu,j}\) operators and not to normal mode superoperators \(\widehat{\zeta}_r\). Using commutation relations, it is straightforward to show that \[\label{eq:zeta95alpha} \widehat{\zeta}_s |\alpha\rangle\!\rangle_r = \alpha |\alpha\rangle\!\rangle_r \delta_{rs}, \qquad \widehat{\zeta}^{+}_s |\alpha\rangle\!\rangle_r = \alpha^* |\alpha\rangle\!\rangle_r \delta_{rs}.\tag{20}\] Defining the "Heisenberg" time-dependent superoperators as \(\widehat{\zeta}_r(t):={\rm e}^{-\widehat{\mathcal{L}}t}\widehat{\zeta}_r \,{\rm e}^{\widehat{\mathcal{L}}t}\), we have \(d\widehat{\zeta}_r(t)/dt=[\widehat{\zeta}_r(t),\widehat{\mathcal{L}}]\). Using commutation relations again, one can verify that \([\widehat{\zeta}_r(t),\widehat{\zeta}_s']=\delta_{r,s}{\rm e}^{ -i\omega_rt}\) and \(\widehat{\zeta}_r(t)=\widehat{\zeta}_r{\rm e}^{-i\omega_r t}\). With this result at hand, we find that the evolution of coherent state operators fulfills \[\begin{align} \label{eq:alphamotion} \widehat{\zeta}_r|\alpha_r(t)\rangle\!\rangle &= \widehat{\zeta}_r {\rm e}^{\widehat{\mathcal{L}}t} |\alpha\rangle\!\rangle_r = {\rm e}^{\widehat{\mathcal{L}}t} \widehat{\zeta}_r(t) |\alpha\rangle\!\rangle_r =\\\nonumber &= \alpha {\rm e}^{-i\omega_r t} |\alpha_r(t)\rangle\!\rangle, \end{align}\tag{21}\] and similarly \(\widehat{\zeta}_r^+|\alpha_r(t)\rangle\!\rangle=\alpha^* {\rm e}^{i\omega_r^* t}|\alpha_r(t)\rangle\!\rangle\), i.e. coherent state operator remains coherent during evolution while the parameter \(\alpha\) evolves with complex frequency \(\omega_r\).

Since \((\!(1|\widehat{\zeta}'_r=0\), primed superoperators are linear superpositions of primed basic superoperators \(\widehat{a}'_{\mu,j}\) defined in Eq. 13 , i.e. \(\widehat{\zeta}'_r=\sum_j(c_{r,j}\widehat{a}_{0,j}'+d_{r,j}\widehat{a}_{1,j}')\). In result, \(|\alpha\rangle\!\rangle_r\) given by Eq. 19 can be, in the case of a pure steady state \(|0_\zeta\rangle\!\rangle=|0_\zeta\rangle\langle 0_\zeta|\), expressed as

where "h.c.\(\!^R\)" is a Hermitian conjugate with operators acting from the right, \(\hat{D}_{r,j}( \alpha) = \exp\left(\alpha \hat{a}^{\dagger}_{r,j} - \alpha^*\hat{a}_{r,j}\right)\) with \(\hat{a}^\dagger_{r,j}=c_{r,j} \hat{a}_j^{\dagger}- d_{r,j} \hat{a}_j\), and \(\hat{D}_r(\alpha)=\prod \hat{D}_{r,j}(\alpha)\) is the displacement operator corresponding to normal mode \(r\), defined as \[\label{eq:ar} \hat{a}_r^\dagger=\sum_j (c_{r,j}\hat{a}_j^\dagger+d_{r,j}\hat{a}_j),\tag{22}\] In the case when the steady state is a mixed state, \(|0_\zeta\rangle\!\rangle=\sum_np_n|n_\zeta\rangle\langle n_\zeta|\), we have \(|\alpha\rangle\!\rangle_r = \sum_n p_n \hat{D}_r(\alpha)|n_\zeta\rangle\langle n_\zeta| \hat{D}^\dagger_r(\alpha)\). We conclude that the dissipative coherent state operator \(|\alpha\rangle\!\rangle_r\) is positive semi-definite, and it corresponds to a physical state.

Clearly, the dissipative coherent state is a displaced nonequilibrium steady state, and it remains so during evolution. This generalizes the known result that pure coherent states remain pure and coherent during single-mode evolution with linear decay [68]. The operators \(\hat{a}_r\) corresponding to different normal modes do not commute with each other in general, \([\hat{a}^\dagger_r,\hat{a}_s]\neq\delta_{r,s}\), since the classical solutions are not necessarily orthogonal in a non-Hermitian system. We emphasize that these coherent states are not necessarily pure states since the nonequilibrium steady state \(|0_\zeta\rangle\!\rangle\) can be a mixed state. Nevertheless, the measured correlation functions exhibit perfect first-order coherence in the large particle number limit, as will be demonstrated below.

3.3 Quantum-classical correspondence↩︎

In analogy to the conservative case, we will demonstrate the correspondence of coherent states to classical states obeying Maxwell equations, but this time including the effects of dissipation. As shown above, \(N\) bosonic modes give rise to \(2N\) normal modes of the Liouvillian. On the other hand, a single bosonic mode corresponds in the classical limit to a pair of Maxwell equations eigenmodes with positive and negative frequency. Hence, the \(2N\) eigenmodes of Maxwell equations have \(2N\) corresponding normal modes of the Liouvillian. To demonstrate direct quantum-classical correspondence, we must find the form of the electric field and matter operators in normal modes and compare them with solutions in the classical limit.

As described in Sec. 3.1, in third quantization operators of observables are constructed by acting from the right side on the identity operator \((\!(1|\), see Eq. 17 . Since electric field intensity and matter particle density are linear in the total particle number, their superoperators must include linear terms in \(\widehat{\zeta}_r\) only. In analogy to the Hermitian case, Eqs. 9 and 10 , we propose that \[\begin{align} (\!(\hat{\boldsymbol{E}}({\boldsymbol{r}})|&=(\!(1|\widehat{\boldsymbol{E}}({\boldsymbol{r}})=(\!(1|\sum_r \left(\widehat{\zeta}_r {\boldsymbol{E}}_r({\boldsymbol{r}}) + \widehat{\zeta}^{+}_r {\boldsymbol{E}}^*_r({\boldsymbol{r}})\right) \tag{23}\\ (\!(\hat{\boldsymbol{X}}({\boldsymbol{r}})|&=(\!(1|\widehat{\boldsymbol{X}}({\boldsymbol{r}})=(\!(1|\sum_r \left(\widehat{\zeta}_r {\boldsymbol{X}}_r({\boldsymbol{r}}) + \widehat{\zeta}^{+}_r {\boldsymbol{X}}^*_r({\boldsymbol{r}})\right). \tag{24} \end{align}\] The first and the second terms on the right hand sides correspond to positive and negative frequency parts, respectively. To check if this guess is correct, let us assume that the state of the system corresponds to the single-mode coherent state \(|\alpha\rangle\!\rangle_r\), and calculate arbitrary first-order correlation functions, in analogy to Eq. 11 \[\begin{align} \label{eq:E1E2zeta} \nonumber G_{\mu,\nu}^{(1)}&({\boldsymbol{r}}_1,{\boldsymbol{r}}_2;t,t) = (\!({ E}_\nu^{L(-)}({\boldsymbol{r}}_1,t) { E}_\mu^{L(+)}({\boldsymbol{r}}_2,t)| \alpha\rangle\!\rangle_r =\\ & = (\!( { E}_{r,\nu}^*({\boldsymbol{r}}_1){ E}_{r,\mu}({\boldsymbol{r}}_2) \widehat{\zeta}^{+}_r\widehat{\zeta}_r| \alpha\rangle\!\rangle_r + G^{(1)}_{NESS} = \\ \nonumber & =|\alpha|^2 { E}_{r,\nu}^*({\boldsymbol{r}}_1){ E}_{r,\mu}({\boldsymbol{r}}_2) + G^{(1)}_{NESS} \\\nonumber &= \langle n \rangle { E}_{r,\nu}^*({\boldsymbol{r}}_1){ E}_{r,\mu}({\boldsymbol{r}}_2) + G^{(1)}_{NESS}, \end{align}\tag{25}\] where \(G^{(1)}_{NESS}\) is the correlation function \(G_{\mu,\nu}^{(1)}({\boldsymbol{r}}_1,{\boldsymbol{r}}_2;t,t)\) in the nonequilibrium steady state \(| 0_\zeta\rangle\!\rangle\). This term appears since the electric field superoperator \(\widehat{\boldsymbol{E}}\) may contain also primed \(\widehat{\zeta}'_r\) terms, which does not contradict Eq. 23 since \((\!(1|\widehat{\zeta}'_r=0\). Nevertheless, in the classical limit we obtain first-order correlation functions that correspond exactly to expected correlations in solutions of Maxwell equations. Clearly, no other choice of space-dependent coefficients in Eqs. 23 and 24 can lead to the correct result.

Consequently, coherent states of normal modes \(|\alpha\rangle\!\rangle_r\) of the above system can be assigned to solutions of Maxwell equations coupled to matter modes, Eqs. (3 ) and (4 ), but with the addition of decay/gain channels. Possible channels include matter decay \(\gamma_x\) and photon decay described with the imaginary part of the background refractive index \(\epsilon_{bI}({\boldsymbol{r}})\). The corresponding system of eigenmode equations reads \[\begin{align} \nabla \times \nabla \times \mathbf{E}_r&=\frac{\left(\epsilon_b({\boldsymbol{r}})+i\epsilon_{bI}({\boldsymbol{r}})\right)\omega_r^2}{\epsilon_0 c^2} \mathbf{E}_r+ \frac{\alpha({\boldsymbol{r}})\omega_r^2}{\epsilon_0 c^2} \mathbf{X}_r, \nonumber\\ \label{eq:MaxwellD} (\omega_0^2 - 2 i \gamma_x \omega_r &-\omega_r^2) \mathbf{X}_r=\frac{\alpha({\boldsymbol{r}})}{\rho} \mathbf{E}_r. \end{align}\tag{26}\]

Moreover, we can include decay resulting from outgoing-wave boundary conditions. In this case, eigenmodes of wave equations such as Eqs. 26 are known as quasinormal modes (QNMs), and have been successfully used to describe systems ranging from black holes to subwavelength plasmonic resonators [69]–[72]. The difficulty with using QNMs lies in the fact that they are divergent at \(|{\boldsymbol{r}}| \to \infty\), which is a direct consequence of natural temporal divergence at \(t \to -\infty\) [73]. However, we can make use of the regularization developed in classical systems to cope with this divergence. Assuming that our system of interest is localized in space, we can surround it with an absorbing material, for example an equivalent of a uniaxial perfectly matched layer (PML), which is a layer of anisotropic material designed to absorb any outgoing waves [74]. This leads to regularized QNMs of Maxwell equations which are exactly equivalent to QNMs of the original problem within the region of interest, do not suffer from the divergence problem, and can be properly normalized [69]. Note that a simple isotropic lossy material would also suffice to absorb outgoing waves, as long as it is "switched on" sufficiently slowly as one moves to \(|{\boldsymbol{r}}| \to \infty\).

The second equation of the system 26 can be incorporated in the first equation as a frequency-dependent contribution to the refractive index. This equation can be then solved analytically or numerically using a Maxwell equation solver that can find eigenmodes in frequency domain with frequency-dependent complex dielectric constant and outgoing wave boundary conditions. Calculated eigenmodes and complex eigenvalues correspond to normal modes \(r\) of the Liouvillian and their complex frequencies via Eqs. 23 and 24 .

In fact, this equivalence of frequency-dependent refractive index and coupling to additional modes is the basis of "auxiliary modes" technique for formulating dispersive Maxwell equations as a non-dispersive problem [75], [76]. Likewise, in our case we can model dispersive background material by auxiliary quantum matter modes which are weakly coupled to light. This way one can introduce frequency-dependent refractive index of the dielectric, in a way similar to that presented in [42], [54].

3.4 Lindbladian form of master equation↩︎

We can transform the Liouvillian 18 into a more familiar Lindbladian-like form that does not involve superoperators on the right hand side.

Let us consider an arbitrary set of \(N\) bosonic operators satisfying canonical commutation relations \([\hat{b}_r,\hat{b}^\dagger_s]=\delta_{rs}\), \([\hat{b}^\dagger_r,\hat{b}^\dagger_s]=[\hat{b}_r,\hat{b}_s]=0,\) with a common vacuum \(\left| 0_b\right\rangle\). We define a new set of superoperators \[\begin{align} \label{eq:eta} &\widehat{\eta}_r=\widehat{b}_r^L,& &\widehat{\eta}'_r=\widehat{b}^{\dagger L}_r-\widehat{b}^{\dagger R}_r, \\\nonumber &\widehat{\eta}^+_r=\widehat{b}_r^{\dagger R},& &\widehat{\eta}'^+_r=\widehat{b}^R_r-\widehat{b}^L_r, \end{align}\tag{27}\] cf. Eq. 13 . It is easy to see that \([\widehat{\eta}_r,\widehat{\eta}_s'] = \delta_{rs}, [\widehat{\eta}'_r,\widehat{\eta}'_s] = [\widehat{\eta}_r,\widehat{\eta}_s] = 0\). The two Fock bases \[\begin{align} \label{eq:twinFock} |\underline{m}\rangle\!\rangle_\zeta =&\prod_{r}\frac{\bigl(\widehat{\zeta}'_{r}\bigr)^{m_{0,r}}}{\sqrt{m_{0,r}!}} \prod_{r}\frac{\bigl(\widehat{\zeta}'^+_{r}\bigr)^{m^+_{1,r}}}{\sqrt{m^+_{1,r}!}}\, |0_\zeta\rangle\!\rangle, \\\nonumber |\underline{m}\rangle\!\rangle_{{\eta}} =&\prod_{r}\frac{\bigl(\widehat{\eta}'_{r}\bigr)^{m_{0,r}}}{\sqrt{m_{0,r}!}} \prod_{r}\frac{\bigl(\widehat{\eta}'^+_{r}\bigr)^{m^+_{1,r}}}{\sqrt{m^+_{1,r}!}}\, |0_b\rangle\!\rangle, \end{align}\tag{28}\] define two Fock spaces \(\mathcal{F}_\zeta\) and \(\mathcal{F}_\eta\), where the \(2N\)-component index \(\underline{m}=(m_{\nu,r};\nu\in\{0,1\},r\in\{1..N\})^T\). Note that \(\mathcal{F}_\zeta\) is identical to the one defined in Eq. 17 . The linear map \(f\) between \(\mathcal{F}_\zeta\) and \(\mathcal{F}_\eta\) such that \(f(| \underline{m}\rangle\!\rangle_\zeta)=| \underline{m}\rangle\!\rangle_{{\eta}}\) is an isomorphism with respect to the action of superoperators. Consider now the following "twin Liouvillian" acting on operators from the space spanned by \(|\underline{m}\rangle\!\rangle_{\eta}\) \[\label{eq:tildeL} \widehat{\mathcal{L}}^b = \sum_{+r}\Bigl(-i\omega_r{\widehat{\eta}}_r'\widehat{\eta}_r +i\omega_r^*{\widehat{\eta}}_r'^+\widehat{\eta}_r^+\Bigr),\tag{29}\] cf. Eq. 18 . Considering that \[\begin{align} &\widehat{\zeta}_r\, |0_\zeta\rangle\!\rangle = 0, \quad &\widehat{\zeta}_r^+\, |0_\zeta\rangle\!\rangle &= 0, \\\nonumber &\widehat{\eta}_r\, |0_b\rangle\!\rangle = \widehat{b}_r^L\,|0_b\rangle\!\rangle = 0, \quad &\widehat{\eta}^+_r\, |0_b\rangle\!\rangle &= \widehat{b}_r^{\dagger R}\,|0_b\rangle\!\rangle = 0 \end{align}\] and that both sets of superoperators fulfill the same commutation relations, it follows that if \(\hat{\rho}_\zeta\) and \(\hat{\rho}_\eta\) are operators from Fock spaces \(\mathcal{F}_\zeta\) and \(\mathcal{F}_\eta\), respectively, such that \(\hat{\rho}_\eta=f(\hat{\rho}_\zeta)\), then \(\widehat{\mathcal{L}}^b\hat{\rho}_\eta=f(\widehat{\mathcal{L}}\hat{\rho}_\zeta)\). Therefore, \(\widehat{\mathcal{L}}^b\) is equivalent conjugate (similar operator) to \(\widehat{\mathcal{L}}\) via isomorphism \(f\). In other words, matrix representation of \(\widehat{\mathcal{L}}\) in the \(| \underline{m}\rangle\!\rangle_{\zeta}\) basis is identical to matrix representation of \(\widehat{\mathcal{L}}^b\) in the \(| \underline{m}\rangle\!\rangle_{\eta}\) basis. Moreover, the isomorphic mapping between \(\hat{\rho}_\eta\) and \(\hat{\rho}_\zeta\) is conserved in time if \(\hat{\rho}_\eta(t)\) evolves according to 29 and \(\hat{\rho}_\zeta(t)\) evolves according to 18 .

Now, we can rewrite \(\widehat{\mathcal{L}}^b\) superoperator in a diagonal form using Eqs.@eq:eq:eta as \(\widehat{\mathcal{L}}^b= \sum_{+r}\widehat{\mathcal{L}}^b_r\), where \(\widehat{\mathcal{L}}^b_r\) contains \(\hat{b}_r\) and \(\hat{b}^\dagger_r\) operators only \[\begin{align} \label{eq:master95eta} \widehat{\mathcal{L}}^b_{r}\hat{\rho}_\eta =& \left(-i\omega_r{\widehat{\eta}}_r'\widehat{\eta}_r +i\omega_r^*{\widehat{\eta}}_r'{}^+\widehat{\eta}_r^+\right)\hat{\rho}_\eta=\\ \nonumber =&\left(-i\omega_r(\widehat{b}_r^{\dagger L }-\widehat{b}_r^{\dagger R })\widehat{b}_r^L + i\omega_r^*(\widehat{b}_r^R - \widehat{b}_r^L)\widehat{b}_r^{\dagger R}\right)\hat{\rho}_\eta,\\ \nonumber =& -i\omega_r(\hat{b}_r^\dagger \hat{b}_r\hat{\rho}_\eta - \hat{b}_r\hat{\rho}_\eta \hat{b}_r^\dagger) + i\omega_r^*(\hat{\rho}_\eta \hat{b}_r^\dagger \hat{b}_r - \hat{b}_r\hat{\rho}_\eta \hat{b}_r^\dagger), \end{align}\tag{30}\] which leads to

Where we omitted the index \(\eta\). In result, evolution governed by a quadratic Liouvillian is isomorphic to evolution under a set of independent single-mode quadratic Liouvillians of the GKLS form. The description of evolution with the Lindbladian equation in the non-Hermitian normal basis is thus just as simple as Lindbladian evolution in the Hermitian basis, cf. Eq. 12 , but has the important advantage of having a diagonal form. This greatly reduces the size of the Hilbert space necessary for analyzing the evolution. That said, one must take care of the nontrivial nature of the vacuum state \(| 0_{{\zeta}}\rangle\!\rangle\) in Eq. 28 and the isomorphism between superoperators when considering specific initial conditions, coupling the system to external modes, or adding non-quadratic terms.

We emphasize that our approach goes beyond the previous attempts of quantization in the basis of quasinormal modes [70], [71], where the resulting Liovillian contained mode-mixing terms, and was therefore non-separable. In our case, both the Hamiltonian and the Lindbladian part of the master equation are obtained in the same diagonal basis, which is the crucial feature allowing for mode separation. It comes at the price of considering an isomorphic representation rather than the original one, but this has no influence on the form of the resulting master equation.

3.5 Summary↩︎

We showed that in analogy to the conservative case, solutions of Maxwell equations coupled to matter form a basis of normal modes of the Liovillian superoperator. This results provides a simple recipe for determining the quantum master equation for a linear system, which consists of: (1) Determining classical solutions of Maxwell equations coupled to matter with appropriate boundary conditions (QNMs), typically by incorporating the response of resonant matter modes into the complex refractive index; (2) Extracting the frequencies with positive real parts and the corresponding eigenmodes and using them as a basis for master equation of the form 31 .

4 Applications↩︎

Transformations introduced in previous sections focused on the determination of a diagonal basis in a quadratic system. While they resulted in optimal basis in quadratic systems, they can still result in a convenient basis in a non-quadratic system. One may draw analogy to the convenient plane-wave basis which is commonly used to analyze systems in which plane waves are not eigenstates of the Hamilonian. Below, we will demonstrate in several practical examples how our approach can be useful in non-quadratic, interacting systems.

4.1 Interactions↩︎

For the interaction Hamiltonian, we consider \[\label{eq:X4} \mathcal{H}_{\rm int}=\frac{g}{6} {\boldsymbol{\hat{X}}}^4\tag{32}\] Substituting Eq. (10 ), we obtain under the RWA \[\label{eq:Hint} H_{\rm int}=g \sum_{jl} \hat{P}^\dagger_l\hat{P}^\dagger_j\hat{P}_j\hat{P}_l \int_V d{\boldsymbol{r}}|{\boldsymbol{X}}_j({\boldsymbol{r}})|^2|{\boldsymbol{X}}_l({\boldsymbol{r}})|^2 + O(n)\tag{33}\] The coefficient \(g\) has now a clear physical interpretation and is simply the interaction constant of the matter component, which corresponds to interaction strength of polaritons when the matter Hopfield coefficient is close to unity. General interaction Hamiltonians with other forms than Eq. 32 may be treated in the same way. Note that we have not made any approximation regarding the weakness of the nonlinearity of the system, and our treatment is still exact. However, the use of the polariton basis \(\hat{P}_i\) will be useful in practice mainly in the case when one is able to identify a small subset of polariton eigenmodes that are important. This is mostly the case when interactions are weak compared to other energy scales in the system.

4.2 Quantum fluctuations↩︎

We consider a single polariton mode, and take into account the lower-order terms that were neglected in the previous subsection (we stay within the RWA) \[H_{\rm int}=g \int_V |{\boldsymbol{X}}({\boldsymbol{r}})|^4 d{\boldsymbol{r}}\left(\hat{P}^\dagger\hat{P}^\dagger\hat{P}\hat{P} + 2 \hat{P}^\dagger\hat{P} + \frac{1}{2}\right).\] The two last terms in brackets can be interpreted as the interaction between non-vacuum polaritons and vacuum fluctuations, and the interactions within the polariton vacuum, respectively. These terms can give rise to physically measurable effects. For instance, the second term leads to interaction-induced shift of the energy of polariton excitations. This can be observed in spectroscopic measurements as a blueshift of the polariton energy in the linear regime, \(\hbar\omega_i^{\rm int}=\hbar\omega_i+2g\int|{\boldsymbol{X}}({\boldsymbol{r}})|^4 d{\boldsymbol{r}}\), where \(\omega_i^{\rm int}\) is the low-density polariton excitation frequency including the effect of interactions. The magnitude of this energy shift is comparable to the single-polariton nonlinearity. It is density-independent, but proportional to the square of the exciton fraction, in contrast to other density-independent effects which are simply proportional to the exciton fraction.

4.3 Exciton confinement↩︎

One of the outstanding challenges in exciton-polariton physics is the achievement of polariton blockade, a regime where strong interactions result in highly non-classical polariton statistics. This regime occurs when polariton-polariton interactions are much stronger than decoherence. Typically, this is the case when the condition \(U/\gamma>1\) is fulfilled with interaction coefficient \(U\) and polariton decay rate \(\gamma\). However, in typical systems one has \(U/\gamma\ll1\). Existing methods and proposals to enhance the interaction strength include strong confinement of the electromagnetic field [56]–[58], [77], inducing electric dipole moment for excitons [15], [59], [60], [78], and exploiting interactions with electrons [9], [61], [79]. In any case, strong confinement of light mode volume seems to be required, but it is difficult to squeeze light much below the wavelength of light. This may be achieved using plasmonic excitations, but at the cost of high loss rates.

Here we consider using weak exciton confinement rather than photon confinement to increase the strength of effective interactions between polaritons. In contrast to quantum dots, by weak exciton confinement we understand confinement on spatial scales much larger than the exciton size, such that excitons can be described as harmonic oscillators, or bosons in the quantum case. Using our theoretical method, we analyze in detail the possibility of obtaining a high ratio of \(U/\gamma\) by engineering the structure geometry in a vertical microcavity structure depicted in Fig. 2 (a). The structure (not shown in the Figure in its full extent) includes a pair of dielectric Bragg mirrors, each composed of 35 \(\rm Al_{0.95}Ga_{0.05}As/Al_{0.2}Ga_{0.8}As\) layer pairs. Between the mirrors, a \(3\lambda/2\) \(\rm Al_{0.95}Ga_{0.05}As\) cavity is placed, with a 10 nm InGaAs QW layer at the antinode of the photonic cavity mode. The transverse size of the QW volume is 0.5\(\mu\)m and the extent in the direction perpendicular to the Figure plane is 1\(\mu\)m. Refractive indices of these materials (including attenuation) are available in public databases, so it is possible to model this structure without fitting parameters. This setting is similar as in [57], but we assume no specific lateral optical confinement. Instead, there is a small volume of "active" material in which the energy of the excitonic excitations is tuned to the energy of the photonic mode, while in the rest of the volume the exciton energy is out of resonance. Such a weak spatial confinement in quantum wells may be achieved e.g. by selective interdiffusion method [80] or additional electrical gates [81].

We determine polariton eigenmodes by solving Maxwell equations coupled to the material response 26 including a complex, space dependent dielectric constant. The calculations of eigenmodes are performed with the plane wave admittance transfer method [82] implemented in the PLaSK software [83] including absorbing boundary conditions with perfectly matching layers. Light-matter coupling strength of the active material was obtained by modeling the system of [57] and reproducing the polariton energy splitting. The interaction and exciton decay rates were also taken from [57]. The solutions of Maxwell equations allow us to estimate the low-intensity mode energy, decay rate, light-matter interaction and polariton-polariton interaction constant using Eq. 33 . To keep computational time in reasonable limits, we determined two-dimensional solutions of Maxwell equations assuming a fixed dimension of \(1\mu\)m in the direction perpendicular to the Figure plane. While this assumption may not be very accurate taking into account the wavelength of light in dielectric \(\lambda/n\approx0.3\mu\)m, the physical principle that we aim to demonstrate should hold also in the fully three-dimensional case.

Figure 2 (b) shows the resulting photonic intensity distribution in the lowest energy polariton mode. The photonic part is confined, but at the spatial scale much larger than the size of the active QW volume, notice the difference in the scales in Fig. 2 (a) and (b). On the other hand, the exciton part is fully confined in the active volume, leading to strong polariton-polariton interactions. According to the results of our theory, a single-mode system can be modeled with quantum master equation of the form \[\begin{align} \dot{\hat{\rho}} =& \bigg(i \Re (\omega)\left[\hat{\rho}, a^\dagger a\right] - \Im(\omega) \left(2a \hat{\rho} a^\dagger -\left\{ a^\dagger a,\hat{\rho}\right\}\right) \bigg) + \nonumber \\&+ U \left(\hat{a}^\dagger\hat{a}^\dagger\hat{a}\hat{a} + 2 \hat{a}^\dagger\hat{a} + \frac{1}{2}\right), \end{align}\] where \(\omega\) is the complex frequency of the polariton mode resulting from Maxwell equations and \(U=g \int_V |{\boldsymbol{X}}({\boldsymbol{r}})|^4 d{\boldsymbol{r}}\). However, in this particular structure, the excitons have a two-dimensional character, so we use a corresponding two-dimensional interaction integral \(U=g_{\rm 2D} \int_S |{\boldsymbol{X}}_\parallel({\boldsymbol{r}})|^4 d{\boldsymbol{r}}_\parallel\) with \({\boldsymbol{X}}_\parallel({\boldsymbol{r}})\) being the two-dimensional polarization density. Indeed, as shown in Fig. 2 (c), the reduction of the active layer lateral size leads to the strong increase of the polariton interaction strength \(U\). At the same time, the mismatch between photonic and excitonic distributions results in the reduction of the light-matter interaction strength \(\Omega\). The loss rate is also slightly increased for smaller layer sizes. However, in a certain range of active layer sizes, the strong interaction regime \(U/\gamma>1\) holds together with the strong coupling condition \(\Omega/\gamma>1\).

In Fig. 2 (d) we present another proposed structure composed of an atomic layer flake of a two-dimensional material \(\rm MoS_2\) sandwiched between \(\rm SiO_2/Si_3N_4\) dielectric Bragg mirrors with 25 layer pairs each. Two-dimensional materials are characterized by a very strong light-matter interaction and the silicon-based Bragg mirrors have very low absorption rates, which may be used to fabricate high-Q cavities [5], [60], [84]–[86]. The calculated light intensity is shown in Fig. 2 (e). In this case, for the \(1\mu\)m\(\times1\mu\)m flake size and non-radiative losses rate \(\gamma_{\rm NR}=10\mu eV\) we estimate the interaction coefficient to be \(U=2.59\mu eV\) and the ratio \(U/\gamma=12.95\).

Finally, Fig. 2 (f) shows a possible waveguide polariton structure with a standing-wave polariton mode in the longitudinal direction, based on \(\rm AlGaAs\) incorporating a 20 nm thick GaAs QW volume. The localized active volume acts as a defect, pinning the polariton mode. Due to the total internal reflection, the total decay rate of polaritons in this configuration is only \(\gamma\approx 7 \mu\)eV assuming nonradiative exciton decay rate of \(10 \mu\)eV. The interaction coefficient depends on the particular configuration used. Particularly interesting are systems with electric field induced dipolar excitons, where the interaction strength may be of the order of 1 meV μm\(^2\) [15]. In this case, the QW layer of size 1 μm\(^2\) results in an interaction coefficient of the order of 1 meV and \(U/\gamma\approx 100\). However, even in the case of non-dipolar excitons, ratios of \(U/\gamma\approx1\) could be achieved.

4.4 Engineering nonlocal interactions↩︎

We extend the structure of the previous section to incorporate two confined, but interacting spatial modes, see Fig. 3 (a). To this end, we incorporated in the cavity two rectangular regions, filled with \(\rm Al_{0.2}Ga_{0.8}As\) material. The refractive index of this material is higher than that of \(\rm Al_{0.95}Ga_{0.05}As\), which results in an effective double well potential for photons in the horizontal direction, as depicted schematically in the inset of Fig. 3 (a). The position and composition of these regions is not crucial, and in practice the potential could be implemented in other ways, e.g. using cavity mesas. The inset also depicts the spatial cross-sections of light amplitude in the symmetric and antisymmetric lower polariton modes obtained numerically. The full light intensity distribution of the symmetric mode is shown in Fig. 3 (b). To enhance non-local interactions, the active material was positioned between the two potential wells. The horizontal width of the active volume marked in blue is equal to 400 nm.

Symmetric mode (S) and antisymmetric mode (AS) form a lower polariton doublet, shifted from other modes (in particular the upper polariton doublet) by 535 \(\mu\)eV, which is much larger than the 142 \(\mu\)eV mode splitting in the doublet. We perform a unitary transformation to the basis of polariton modes localized in the left an right potential well via \(\hat{P}_{\rm L,R}=(1/\sqrt{2})(\hat{P}_{\rm S}\mp \hat{P}_{\rm AS})\). After this transformation, master equation takes the form \[\begin{align} \dot{\hat{\rho}} =\frac{1}{i\hbar} [\hat{H},\hat{\rho}] +\sum_{j=L,R} \gamma_j \left(a_j \hat{\rho} a_j^\dagger -\frac{1}{2}\left\{ a_j^\dagger a_j,\hat{\rho}\right\}\right), \end{align}\] where \[\begin{align} \label{eq:twomode} \hat{H} &= \hbar \omega_{LR} \left(a_{\rm L}^\dagger a_{\rm L} + a_{\rm R}^\dagger a_{\rm R}\right) + \\ &+ J \left(a_{\rm L}^\dagger a_{\rm R} + a_{\rm R}^\dagger a_{\rm L}\right)+ \frac{1}{2} \Big(U_{LL} a_{\rm L}^\dagger a_{\rm L}^\dagger a_{\rm L} a_{\rm L} + \\ & + U_{RR} a_{\rm R}^\dagger a_{\rm R}^\dagger a_{\rm R} a_{\rm R} + 2 U_{LR} a_{\rm L}^\dagger a_{\rm R}^\dagger a_{\rm R} a_{\rm L}\Big).\nonumber \end{align}\tag{34}\] Here, localized mode energy is \(\omega_{\rm LR}=\Re (\omega_{\rm S}+\omega_{\rm AS})/2\), coupling coefficient is \(J=\Re (\omega_{\rm S}-\omega_{\rm AS})/2\) and loss coeffcients are \(\gamma_{\rm L,R}=-\Im (\omega_{\rm S}+\omega_{\rm AS})\) while the interaction coefficients are \(U_{ij}=g \int_V |{\boldsymbol{X}_i}({\boldsymbol{r}})|^2 |{\boldsymbol{X}_j}({\boldsymbol{r}})|^2 d{\boldsymbol{r}}\) with \(U_{LL}=U_{RR}=:U_{11}\) due to mirror symmetry and \(U_{LR}=:U_{12}\). We neglected interactions with quantum fluctuations leading to constant frequency shifts. Due to the weakness of dissipative effects compared to mode energies, in this case we can also neglect higher order effects discussed in Sec. 3.4.

We calculate the second-order equal-time polariton correlation functions in a steady state under coherent pumping by adding the term \(\sum_{i=L,R}\Omega_p(e^{i\omega_p t}\hat{P}_i+e^{-i\omega_p t}\hat{P}_i^\dagger)\) to the Hamiltonian. Since polariton energy is much higher than the light-matter interaction, we work under RWA and transform the Hamiltonian to the frame rotating with pump frequency \(\omega_p\), introducing the detuning \(\delta=\hbar(\omega_{\rm LR}-\omega_p)\). Steady state intra-mode and inter-mode correlations are shown in Fig. 3 (c) and Fig. 3 (d), respectively. The system exhibits standard bunching/antibunching behavior when crossing the resonance [87]. Interestingly, cross-mode correlations between left and right polariton modes \(g_{12}^{(2)}(0)\) exhibit also a strong dependence on detuning, showing strong preference for simultaneous population in both modes (positive cross-correlation) at large detunings. The range of parameters for which the Cauchy-Schwartz inequality is violated, marking the emergence of quantum entanglement, is marked with shaded area. This phenomenon is significantly boosted by non-local interactions, as evident from the comparison to the results of simulations where the cross-interaction term was artificially removed, see dashed line in Fig. 3 (d). This illustrates that engineering nonlocal interactions can have a substantial effect on polariton statistics in nanostructures, and could be used to construct sources of quantum light.

5 Conclusions↩︎

In conclusion, we introduced a method for obtaining a quantum model corresponding to a confined coupled light-matter system, both in the conservative and dissipative regime, with no assumptions on the system geometry. The method is based on quantization in normal modes basis, and allows to treat without approximations systems that are not well described by the Hopfield model. It requires solving the classical limit equations to determine the complete form of the quantum model in an optimal basis, with no fitting parameters. Even though the method is based on diagonalization of a quadratic system, we showed how it can be used in practice to describe non-quadratic, interacting systems. We proposed practical applications resulting in the increase of effective polariton interactions and engineering non-local interactions in semiconductor exciton-polariton structures.

In comparison to previous approaches [42]–[54], the main advantage of our treatment is the simplicity of calculations and of the resulting model. We also emphasize that the applicability of the presented method is much broader than the context presented here. For example, one can go beyond the assumption of small emitter size, and the assumption of infinite emitter mass, as long as the classical description of the system is known. It is also possible to add to the model interactions with other relevant excitations, such as phonons. As long as the system is well described as a collection of bosonic modes, the quantum-to classical correspondence can be used to determine the form of the quantum model. Moreover, even when one needs to describe modes which are not bosonsic (for example, spin systems), analogous transformations may be found [88].

The presented theoretical method can be used to engineer light-matter systems that include semiconductor nanostructures and cold atom systems. Possible practical applications include engineering systems for quantum simulations [36], where reproducing a given, precisely defined quantum model is of great importance, and designing sources of quantum light. Finally, our method can provide better understanding of nonconventional polariton systems that are not well described within simplified theoretical models such as the Hopfield model.

This project received funding from the European Union’s Horizon Europe research and innovation programme under grant agreements No. ID 101115575 (Q-ONE) and ID 101130384 (QUONDENSATE). AO acknowledges support from the National Science Center, Poland, project No. 2024/52/C/ST3/00324.

References↩︎