text:
# Superradiance of low density Frenkel excitons in a crystal slab of three-level atoms: Quantum interference effect

## I Introduction

It is well known that the fluorescence of an exciton exhibits super-radiant character due to the appearance of macroscopic transition dipole moment of the exciton . However this collective feature of exciton radiance depends upon dimensionality of crystal. In infinite bulk crystals, the exciton does not radiate because this exciton is dressed with a photon which has the same wave vector due to the transitional symmetry of the total system. As a result, a stable polariton is formed . In the case of lower dimensional systems, one can show that the exciton decays superradiantly due to breakdown of the translational symmetry. For example the decay rates of exciton are of the order $`(\lambda /a)\gamma `$ for one-dimensional (1D) crystals and $`(\lambda /a)^2\gamma `$ for 2D crystals , with $`\gamma `$ being the radiative decay rate for an isolated atom, $`a`$ the lattice constant, and $`\lambda `$ the light wavelength.

The enhanced factor $`(\lambda /a)^D`$ in the radiative width of exciton is now regarded commonly as one of evidences of “superradiance” , which has been demonstrated in experiments for Frenkel excitons  and Wannier excitons . Since the enhanced radiant effect can also appear in the usual superradiance of an ensemble of atoms , how can we distinguish the different features between the fluorescence of the excitons and that of the atomic ensemble? It is known that, in the fully population-inverted atom systems, the atoms radiate independently with each other in the initial stage. The back-action of the emitted photon to the atoms results in the correlation among atoms. Consequently, the atoms become cooperative and thus the fluorescence from the atomic ensemble will show different statistical properties in the initial and final stages. On the other hand, the exciton fluorescence exhibits identical statistical character during the whole process . Physically, this is because the initial dipole moments of the atoms are spatially random in an atomic ensemble, but in a semiconductor crystal, the dipole moment of exciton presents a macroscopic effect even at initial moment. So the optical properties of multidimensional quantum-confined semiconductor structures (MQCS), such as quantum wells, quantum wires and quantum dots, have their own specific features in physical processes.

The optical properties of the MQCS in a semiconductor microcavity (SMC) have attracted more and more attention in the past decades . The SMC with high-reflectivity dielectric mirrors leads to the realization of the strong coupling between radiation and matter. Moreover, since the optical mode structure of the SMC may alter around the MQCS, many new phenomena, such as tailoring the spontaneous-radiation rate and pattern , the coupled exciton-photon mode splitting in a SMC , have been demonstrated. The resonant interaction between excitons and a single-mode cavity field and the corresponding detuning effect were further investigated .

Most of these former works mentioned above dealed only with the two-level lattice atom case, however the three-level atom case may be very useful to implement quantum information encoding and processing . Over the last few years, the cavity QED with the collective excitations of ensembles of three-level atoms has attracted much attention for quantum computing implementations. In this case, many atoms are entangled through their interaction with the common cavity field. To maintain quantum coherence in this quantum information processing , it is important to reach the so-called strong coupling regime where the single-photon coherent coupling $`g_0\gamma \text{}\gamma _{\text{cav}}`$, the atomic and cavity dissipation (decoherence) rates respectively. It is the symmetric collective excitation that can reach the strong coupling regime without requiring a high finesse cavity as $`g_0\sqrt{N}`$, with the total number of atoms $`N`$ and here we need to consider the quantum decoherence induced by the spatial inhomogeneousness of couplings . In fact, the phenomena of superfluorescence or superradiance  constitute another example of collective state dynamics. Recent experimental success clearly demonstrates the power of such an atomic ensemble based on the system for entangling macroscopic objects .

In the present paper, we study the fluorescence of low density excitons in a crystal slab containing V-type three-level atoms. The purpose of this paper is to investigate the quantum interference effect  in the time evolution of light intensity. In sec. II, considering that the existence of the two-level atomic exciton is mathematically associated with the infinite dimensional reducible representation of SU(2) Lie algebra , we can rationally define the exciton operators for the three-level case associated with SU(3) algebra, which contains various SU(2) subalgebras. With this conception, both the free part and the interaction part of Hamiltonian can be written down in terms of the introduced two-mode exciton operators, which can be described by bosonic operators in the low-density limit. In sec. III, with rotating wave approximation (RWA) for the interaction $`(e/mc)𝐏𝐀`$, the coupled equations for the exciton-light field system are solved with the help of Weisskopf-Wigner approach (WWA) . In sec. IV, the explicit expressions of the electric field operators are derived for the monolayer case. The light intensity as well as the first- and second-order degree of quantum coherence are calculated to show certain new features of exciton fluorescence in a three-level crystal slab. We discuss the phenomenon of quantum beat in the time evolution of the light intensity for various initial exciton states. In sec. V, we consider the roles of both the non-RWA terms and the self-interaction term of photon. This consideration avoids “unphysical” roots of the characteristic equation when the non-perturbation approach given in Ref.  is used. The light intensity is calculated to compare with the results in sec. IV. The conclusions are presented in Sec. VI with some remarks.

## II two mode exciton system with SU(3) structure

We consider a plane crystal slab with a simple cubic structure which contains a stack of $`N`$ identical layers. V-type three-level atoms, as shown in Fig. 1, occupy $`N_T`$ lattice sites, where $`N_T=N_LN`$ and $`N_L`$ is the total lattice sites within each layer. The wave vectors of the excitons and light fields are all assumed perpendicular to the slab. We restrict ourselves to investigate only the low density excitons region.

The interaction $`\frac{e}{mc}𝐩𝐀`$ between the radiation field and the multi-atom system is written in the second quantization form

$`\widehat{H}_1`$ $`=`$ $`\mathrm{}{\displaystyle \underset{q;l,j}{}}g_1(q)|e_1_{lj}g|\left[\widehat{a}_q+\widehat{a}_q^{}\right]e^{iqla}`$ (2)
$`+\mathrm{}{\displaystyle \underset{q;l,j}{}}g_2(q)|e_2_{lj}g|\left[\widehat{a}_q+\widehat{a}_q^{}\right]e^{iqla}+\mathrm{H}.\mathrm{c}.,`$
where $`\widehat{a}_q`$ and $`\widehat{a}_q^{}`$ are the annihilation and creation operators of the photon with wave vector $`q`$ along $`z`$, respectively, $`j`$ and $`l`$ denote the $`j`$th lattice site in the $`l`$th layer. $`g_1(q)=\sqrt{\frac{2\pi \mathrm{\Omega }_1}{V\mathrm{}|q|c}}d_1`$, $`g_2(q)=\sqrt{\frac{2\pi \mathrm{\Omega }_2}{V\mathrm{}|q|c}}d_2`$ are the effective atom-photon coupling constants of $`|g|e_1`$ and $`|g|e_2`$ atomic transitions with transition frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, respectively. The corresponding transition dipole moments are represented by $`𝐝_1`$ and $`𝐝_2`$. We assume that $`𝐝_1`$ and $`𝐝_2`$ are parallel and lying in the slab plane.

Introduce the collective operators for the $`l`$th layer as

$`\widehat{A}^{(l)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_L}}}{\displaystyle \underset{j}{}}|e_1_{lj}g|,\widehat{A}^{(l)}={\displaystyle \frac{1}{\sqrt{N_L}}}{\displaystyle \underset{j}{}}|g_{lj}e_1|,`$ (3)
$`\widehat{B}^{(l)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_L}}}{\displaystyle \underset{j}{}}|e_2_{lj}g|,\widehat{B}^{(l)}={\displaystyle \frac{1}{\sqrt{N_L}}}{\displaystyle \underset{j}{}}|g_{lj}e_2|.`$ (4)
It is easy to prove that, two sets of the operators

$`E_1^{\left(l\right)}`$ $`=`$ $`{\displaystyle \underset{j}{}}|e_1_{lj}g|,F_1^{\left(l\right)}=E_1^{\left(l\right)},`$ (5)
and

$`F_2^{\left(l\right)}`$ $`=`$ $`{\displaystyle \underset{j}{}}|e_2_{lj}g|,E_2^{\left(l\right)}=F_2^{\left(l\right)},`$ (6)
just generate two SU(2) algebras not commuting with each other. This means the four collective operators $`E_1^{\left(l\right)}`$, $`E_2^{\left(l\right)}`$, $`F_1^{\left(l\right)}`$ and $`F_2^{\left(l\right)}`$ do not span a product algebra SU(2)$``$SU(2). A straightforward calculation in Appendix A checks that they satisfy SU(3) algebra. Actually, the above four collective operators define a spinor realization of SU(3) of $`N_L+1`$ dimensions. Furthermore, the unique number $`N_L`$, atom number, determines the dimensions $`N_L+1`$ of representations of the two SU(2) subalgebras. Since we understood the single-mode exciton in terms of the large $`N_L`$ limit of representations of SU(2), it is easy to prove for SU(3) case that, in a very large $`N_L`$ and the low density excitation region, the above defined collective operators $`\widehat{A}^{(l)}`$ and $`\widehat{B}^{(l)}`$ become the bosonic ones obeying the commutation relation

$$[\widehat{A}^{(l)},\widehat{A}^{(l^{})}]=\delta _{l,l^{}},[\widehat{B}^{(l)},\widehat{B}^{(l^{})}]=\delta _{l,l^{}},$$
(7)
where the ideal bosonic approximation is equivalent to the neglect of the phase space filling effect and the exciton-exciton interaction .

In terms of these collective operators, the two-mode Frenkel exciton operators with wave vectors $`k=\frac{2\pi m}{Na}`$ ($`m=\frac{1}{2}(N1),\frac{1}{2}(N3),\mathrm{},\frac{1}{2}(N1)`$), are just the discrete Fourier transformations for them

$`\widehat{A}_k`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{l}{}}e^{ikla}\widehat{A}^{(l)},`$ (8)
$`\widehat{B}_k`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{l}{}}e^{ikla}\widehat{B}^{(l)}.`$ (9)
In fact their conjugates $`A_k^{}`$ and $`B_k^{}`$ are just the generators for the quasi-spin wave states

$`|A_k=\widehat{A}_k^{}|g_1,g_2,\mathrm{},g_{N_T},`$ (10)
$`|B_k=\widehat{B}_k^{}|g_1,g_2,\mathrm{},g_{N_T}.`$ (11)
Since the operators $`\widehat{A}_k`$ and $`\widehat{B}_k`$ commute with each other for very large $`N_L`$ and low excitation, they form a two independent boson system. Then we obtain the interaction Hamiltonian for the the two-mode Frenkel exciton system coupled to a quantized electromagnetic field

$`\widehat{H}_1`$ $`=`$ $`\mathrm{}{\displaystyle \underset{q,k}{}}G_1(q)O(k+q)\left[\widehat{A}_k+\widehat{A}_k^{}\right]\left[\widehat{a}_q+\widehat{a}_q^{}\right]`$ (12)
$`+`$ $`\mathrm{}{\displaystyle \underset{q,k}{}}G_2(q)O(k+q)\left[\widehat{B}_k+\widehat{B}_k^{}\right]\left[\widehat{a}_q+\widehat{a}_q^{}\right],`$ (13)
where the coupling constants between the photons and excitons take the following form

$`G_1(q)=\sqrt{N_T}g_1(q)=\sqrt{N_T{\displaystyle \frac{2\pi \mathrm{\Omega }_1^2}{V\mathrm{}\omega _q}}}d_1,`$ (14)
$`G_2(q)=\sqrt{N_T}g_2(q)=\sqrt{N_T{\displaystyle \frac{2\pi \mathrm{\Omega }_2^2}{V\mathrm{}\omega _q}}}d_2,`$ (15)
here $`\omega _q=|q|c`$. The wave-vector matching factor  in Eq. (8) is

$$O(k+q)=\frac{1}{N}\underset{l}{}e^{i(k+q)la}=\frac{1}{N}\frac{\mathrm{sin}(\frac{k+q}{2}Na)}{\mathrm{sin}(\frac{k+q}{2}a)},$$
(16)
which is real and equal to 1 for $`k+q=0`$, and $`O(k+q)<1`$, for $`k+q0`$.

## III Equations of motion with rotating wave approximation

In this section, we consider the rotating wave approximation (RWA) to deal with the interaction $`\frac{e}{mc}𝐩𝐀`$ between the radiation field and the multi-atom system. After introducing the two-mode exciton operators in Eq. (6), the interaction Hamiltonian between the excitons and photons with RWA is obtained as

$`\widehat{H}_{RWA}`$ $`=`$ $`\mathrm{}{\displaystyle \underset{q,k}{}}G_1(q)O(k+q)\left[\widehat{A}_k\widehat{a}_q^{}+\widehat{A}_k^{}\widehat{a}_q\right]`$ (18)
$`+\mathrm{}{\displaystyle \underset{q,k}{}}G_2(q)O(k+q)\left[\widehat{B}_k\widehat{a}_q^{}+\widehat{B}_k^{}\widehat{a}_q\right],`$
where we have neglected the higher frequency (non-resonant) terms: $`\widehat{A}_k^{}\widehat{a}_q^{}`$, $`\widehat{A}_k\widehat{a}_q`$, $`\widehat{B}_k^{}\widehat{a}_q^{}`$, and $`\widehat{B}_k\widehat{a}_q`$.

The Heisenberg equations for the exciton and photon operators are

$`i{\displaystyle \frac{}{t}}\widehat{A}_k`$ $`=`$ $`\mathrm{\Omega }_1\widehat{A}_k+{\displaystyle \underset{q}{}}G_1(q)O(qk)\widehat{a}_q,`$ (19)
$`i{\displaystyle \frac{}{t}}\widehat{B}_k`$ $`=`$ $`\mathrm{\Omega }_2\widehat{B}_k+{\displaystyle \underset{q}{}}G_2(q)O(qk)\widehat{a}_q,`$ (20)
$`i{\displaystyle \frac{}{t}}\widehat{a}_q`$ $`=`$ $`\omega _q\widehat{a}_q+G_1(q){\displaystyle \underset{k}{}}O(kq)\widehat{A}_k`$ (22)
$`+G_2(q){\displaystyle \underset{k}{}}O(kq)\widehat{B}_k.`$
Taking the following transform

$`\widehat{A}_k\stackrel{~}{A}_k=\widehat{A}_ke^{i\mathrm{\Omega }_1t},\widehat{B}_k\stackrel{~}{B}_k=\widehat{B}_ke^{i\mathrm{\Omega }_2t},`$ (23)
to remove the fast varying factors, we obtain the formal equation for the exciton operator

$`{\displaystyle \frac{}{t}}\stackrel{~}{A}_k(t)`$ $`=`$ $`i{\displaystyle \underset{q}{}}G_1(q)O(qk)e^{i(\omega _q\mathrm{\Omega }_1)t}\widehat{a}_q(0)`$ (28)
$`{\displaystyle \underset{q,k^{}}{}}G_1(q)G_2(q)O(qk)O(k^{}q)`$
$`\times {\displaystyle _0^t}\stackrel{~}{B}_k^{}(t^{})e^{i(\omega _q\mathrm{\Omega }_1)t}e^{i(\omega _q\mathrm{\Omega }_2)t^{}}dt^{}`$
$`{\displaystyle \underset{q,k^{}}{}}G_1^2(q)O(qk)O(k^{}q)`$
$`\times {\displaystyle _0^t}\stackrel{~}{A}_k^{}(t^{})e^{i(\omega _q\mathrm{\Omega }_1)(tt^{})}dt^{}.`$
Here, the first term proportional to $`\widehat{a}_q(0)`$ is the so-called quantum noise term. The second term in the above equation corresponds to a multi-photon processes (MPP) including stimulated emission and absorption effects. Its contribution can be ignored since it is a higher order term from the standpoint of perturbation . The last term in Eq. (28) can be solved by using the Weisskopf-Wigner approximation (WWA), i.e., assuming that $`\stackrel{~}{A}_k^{}(t^{})`$ varies sufficiently slowly so that it can be factorized outside the integral. The remaining part of the time integral of the last term in Eq. (28) can be evaluated and get a Dirac $`\delta `$ function with variable $`(\omega _q\mathrm{\Omega }_1)`$ and a principal part $`𝒫\left[\frac{i}{\omega _q\mathrm{\Omega }_1}\right]`$ term, which contributes a frequency shift (Lamb shift).

The equation of motion for $`\stackrel{~}{B}_k(t)`$ can also be obtained in the similar way. By using the WWA, and neglecting the MMP, we can solve the equations for both of the two-mode excitons to obtain $`\stackrel{~}{A}_k(t)`$ and $`\stackrel{~}{B}_k(t)`$. Substituting them into Eq. (22), the explicit expression for $`\widehat{a}_q(t)`$ in terms of $`\widehat{a}_q(0)`$, $`\widehat{A}_k(0)`$ and $`\widehat{B}_k(0)`$ can be obtained. Finally we can get the positive-frequency part of the electric field

$`\widehat{E}^{(+)}(z,t)`$ $`=`$ $`i{\displaystyle \underset{q}{}}\sqrt{{\displaystyle \frac{2\pi \mathrm{}\omega _q}{V}}}\widehat{a}_q(t)e^{iqz}`$ (29)
$`=`$ $`i{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q\sqrt{{\displaystyle \frac{\mathrm{}\omega _qL}{2\pi A}}}\widehat{a}_q(t)e^{iqz},`$ (30)
where we have taken the photon normalization volume $`V`$ to be $`LA`$, with $`A`$ being the area of the crystal slab, and assumed that the slab is located at the middle of the volume. When $`L`$ is sufficiently large, the sum over $`q`$ has been replaced by an integral: $`_q\mathrm{}\frac{L}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑q\mathrm{}`$.

For an arbitrary initial state $`|\varphi =|\varphi _{ex}|\varphi _L`$ of the total system, the light intensity radiated from the two-mode exciton system is defined as

$$I(z,t)=\frac{c}{2\pi }\varphi |:\widehat{E}^{()}(z,t)\widehat{E}^{(+)}(z,t):|\varphi ,$$
(31)
where the symbol “: … :” means the normal product according to the exciton operators. $`z`$ is the position of a detector relative to the crystal slab. The first-order degree of coherence of fluorescence is given by

$$g^{(1)}(z;t,t+\tau )=\frac{\frac{c}{2\pi }\varphi |:\widehat{E}^{()}(z,t)\widehat{E}^{(+)}(z,t+\tau ):|\varphi }{\sqrt{I(z,t)I(z,t+\tau )}}.$$
(32)
Similarly, one can also define the second-order degree of coherence as

$`g^{(2)}(z;t,t+\tau )={\displaystyle \frac{\frac{c^2}{4\pi ^2}\varphi |:\widehat{E}^{()}(z,t)\widehat{E}^{()}(z,t+\tau )\widehat{E}^{(+)}(z,t+\tau )\widehat{E}^{(+)}(z,t):|\varphi }{I(z,t)I(z,t+\tau )}}.`$ (33)
## IV Quantum interference effect in light intensity

In this section, the features of exciton fluorescence in the case of monolayer will be studied in detail. From the definition of the wave vector of excitons, when $`N=1`$, it takes only one value $`k=0`$ and the wave-vector matching factor $`O(qk)=O(k^{}q)1`$. Thus with the help of WWA, we get the equations of motion for exciton modes $`A_0`$ and $`B_0`$ as

$`{\displaystyle \frac{}{t}}\stackrel{~}{A}_0(t)`$ $`=`$ $`i{\displaystyle \underset{q}{}}G_1(q)e^{i(\omega _q\mathrm{\Omega }_1)t}\widehat{a}_q(0){\displaystyle \frac{\eta _1}{2}}\stackrel{~}{A}_0(t),`$ (34)
$`{\displaystyle \frac{}{t}}\stackrel{~}{B}_0(t)`$ $`=`$ $`i{\displaystyle \underset{q}{}}G_2(q)e^{i(\omega _q\mathrm{\Omega }_2)t}\widehat{a}_q(0){\displaystyle \frac{\eta _2}{2}}\stackrel{~}{B}_0(t),`$ (35)
in which,

$`\eta _1`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{\Omega }_1d_1^2}{\mathrm{}a^2c}}=3\left({\displaystyle \frac{\pi \lambda _1^2}{a^2}}\right)\gamma _1,`$ (36)
$`\eta _2`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{\Omega }_2d_2^2}{\mathrm{}a^2c}}=3\left({\displaystyle \frac{\pi \lambda _2^2}{a^2}}\right)\gamma _2,`$ (37)
where, $`\gamma _1=\frac{4\mathrm{\Omega }_1^3d_1^2}{3\mathrm{}c^3}`$ is the radiative decay rate of an isolated atom from $`|e_1`$ to $`|g`$, and $`\gamma _2=\frac{4\mathrm{\Omega }_2^3d_2^2}{3\mathrm{}c^3}`$ is that of $`|e_2`$ to $`|g`$ , respectively. $`\lambda _j=\frac{c}{\mathrm{\Omega }_j}`$, for $`j=1\text{}2`$, denote the reduced wavelengthes. We find that, because of the implementation of WWA and the ignorance of multi-photon scattering processes (MPP) in the derivation of Eqs. (34) and (35), the two-mode excitons decay independently in an exponential rule. Even though our treatment seems to be bold, we can get a useful result that the decay rates of excitons in two dimensional crystal slab are proportional to the enhanced factor $`(\lambda /a)^2`$. In section V, we will consider a more complex case, in which a non-perturbation approach is used to restudy the radiative decay rates of the two-mode excitons without WWA, and the effects of MPP will be also investigated there.

Substituting the solutions of Eqs. (34) and (35) into Eq. (22), we get the explicit expression of the photon annihilation operator $`\widehat{a}_q(t)`$ for the monolayer case:

$`\widehat{a}_q(t)`$ $`=`$ $`u_q(t)\widehat{A}_0(0)+v_q(t)\widehat{B}_0(0)`$ (39)
$`+e^{i\omega _qt}\widehat{a}_q(0)+{\displaystyle \underset{q^{}}{}}w_{q,q^{}}(t)\widehat{a}_q^{}(0),`$
where

$`u_q(t)`$ $`=`$ $`G_1(q){\displaystyle \frac{e^{i\omega _qt}e^{\eta _1t/2}e^{i\mathrm{\Omega }_{01}t}}{\omega _q\mathrm{\Omega }_{01}+i\eta _1/2}},`$ (40)
$`v_q(t)`$ $`=`$ $`G_2(q){\displaystyle \frac{e^{i\omega _qt}e^{\eta _2t/2}e^{i\mathrm{\Omega }_{02}t}}{\omega _q\mathrm{\Omega }_{02}+i\eta _2/2}},`$ (41)
$`w_{q,q^{}}(t)`$ $`=`$ $`{\displaystyle \frac{G_1(q)G_1(q^{})e^{i\omega _qt}}{\omega _q^{}\mathrm{\Omega }_{01}+i\eta _1/2}}[{\displaystyle \frac{e^{i(\omega _q\mathrm{\Omega }_{01})t}e^{\eta _1t/2}1}{\omega _q\mathrm{\Omega }_{01}+i\eta _1/2}}`$ (43)
$`{\displaystyle \frac{e^{i(\omega _q\omega _q^{})t}1}{\omega _q\omega _q^{}}}]+\text{same with }12,`$
in which $`\mathrm{\Omega }_{0j}=\mathrm{\Omega }_j+\mathrm{\Delta }_j`$, for $`j=1\text{}2`$, are the renormalized physical frequencies, and $`\mathrm{\Delta }_j`$ are the lamb frequency shifts

$$\mathrm{\Delta }_j=𝒫_0^{\mathrm{}}𝑑\omega _q\frac{\rho (\omega _q)\left|G_j(q)\right|^2}{\omega _q\mathrm{\Omega }_j}.$$
(44)
Here “$`𝒫`$” denotes for the Cauthy principle part, and $`\rho (\omega _q)`$ is the density of states of the radiation field.

Therefore, one can obtain the positive frequency part of the electric field operator

$`\widehat{E}^{(+)}(z,t)`$ $`=`$ $`\widehat{E}_0^{(+)}(z,t)+\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_1\eta _1}{4Ac}}}F_A(z,t)\widehat{A}_0(0)`$ (46)
$`+\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_2\eta _2}{4Ac}}}F_B(z,t)\widehat{B}_0(0),`$
where

$`\widehat{E}_0^{(+)}(z,t)`$ $`=`$ $`i{\displaystyle \underset{q}{}}\sqrt{{\displaystyle \frac{2\pi \mathrm{}\omega _q}{V}}}[\widehat{a}_q(0)e^{i\omega _qt}`$ (48)
$`+{\displaystyle \underset{q^{}}{}}w_{q,q^{}}(t)\widehat{a}_q^{}(0)]e^{iqz}.`$
The first term in $`\widehat{E}_0^{(+)}(z,t)`$ is just the free varying photon field. The second term is proportional to the square of coupling constants. As mentioned in Eq. (17), the sum over the wave vector $`q`$ can be replaced by an integral, thus two time-dependent functions $`F_A(z,t)`$ and $`F_B(z,t)`$ in Eq. (30) are

$`F_A(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega _q2\mathrm{cos}\left(\omega _q{\displaystyle \frac{z}{c}}\right){\displaystyle \frac{e^{i\omega _qt}e^{i\omega _1t}}{\omega _q\omega _1}},`$ (49)
$`F_B(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega _q2\mathrm{cos}\left(\omega _q{\displaystyle \frac{z}{c}}\right){\displaystyle \frac{e^{i\omega _qt}e^{i\omega _2t}}{\omega _q\omega _2}}.`$ (50)
Here $`\omega _1=\mathrm{\Omega }_{01}i\eta _1/2`$, $`\omega _2=\mathrm{\Omega }_{02}i\eta _2/2`$. In the Ref. , similar expressions were presented. However, the authors of Ref.  did not calculate the integrals over $`d\omega _q`$, which may result in the light intensity going into the spacelike region. In this paper, we will further calculate the integral over $`d\omega _q`$  to give an explicit expression of $`\widehat{E}^{(+)}(z,t)`$.

$`F_A(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^t}𝑑t^{}e^{\eta _1t^{}/2}e^{i\mathrm{\Omega }_{01}t}{\displaystyle _0^{\mathrm{}}}𝑑\omega _q2\mathrm{cos}\left(\omega _q{\displaystyle \frac{z}{c}}\right)`$ (52)
$`\times e^{i(\omega _q\mathrm{\Omega }_{01})(t^{}t)}.`$
Replacing the integral variable $`\omega _q`$ by $`\omega _q\mathrm{\Omega }_{01}`$, we have

$`F_A(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^t}𝑑t^{}e^{\eta _1t^{}/2}e^{i\mathrm{\Omega }_{01}t}{\displaystyle _{\mathrm{\Omega }_{01}}^{\mathrm{}}}𝑑\omega _q`$ (54)
$`\times 2\mathrm{cos}\left[{\displaystyle \frac{\left(\omega _q+\mathrm{\Omega }_{01}\right)z}{c}}\right]e^{i\omega _q(t^{}t)}.`$
We assume that $`\mathrm{\Omega }_{01}`$ is much larger than all other quantities of the dimension of frequency , so the lower limit of the integration can be extended to $`\mathrm{}`$. Therefore, the above integral over $`d\omega _q`$ can be approximated by two $`\delta `$ functions with variables $`(t^{}t+\frac{z}{c})`$ and $`(t^{}t\frac{z}{c})`$. In the region $`z>0`$ outside the layer, performing the integral over $`t^{}`$, we get

$`F_A(z,t)`$ $`=`$ $`2e^{\eta _1(t\frac{z}{c})/2}e^{i\mathrm{\Omega }_{01}(t\frac{z}{c})}\mathrm{\Theta }\left(t{\displaystyle \frac{z}{c}}\right)`$ (55)
$`=`$ $`2e^{i\omega _1(t\frac{z}{c})}\mathrm{\Theta }\left(t{\displaystyle \frac{z}{c}}\right),`$ (56)
where $`\mathrm{\Theta }`$ is a Heaviside function. Similarly we can also obtain

$$F_B(z,t)=2e^{i\omega _2(t\frac{z}{c})}\mathrm{\Theta }\left(t\frac{z}{c}\right).$$
(57)
From Eqs. (56) and (57), one can find that there is no advanced propagator proportional to $`e^{i(\mathrm{})(t+\frac{z}{c})}`$ appeared in the positive frequency part of the electric field operator $`\widehat{E}^{(+)}(z,t)`$. So that the system is causal and well-behaved .

It should be pointed out that $`F_A(z,t)`$ and $`F_B(z,t)`$ in Eq. (56) and Eq. (57) show that the electric field of the fluorescence emitted from the two-mode excitons is a temporal damped plane wave, but does not decay with $`z`$. A realistic electric field, however, should decrease with $`z`$, which happens, especially when a light propagates from a medium with low refractive index to that of higher one. Besides, a realistic electric field should not simply be a plane wave. However, the theoretical calculations with this subtle consideration will be very complicated and do not provide us with necessary physics in general. Our ideal assumption of a plane wave propagating in $`z`$ may work well for longer wave length radiation near the crystal.

The vacuum state of the exciton may be defined as

$$\widehat{A}_0(0)|0=\widehat{B}_0(0)|0=0,$$
(58)
which represents that there is no exciton contained in the crystal slab. Then number states for the two-mode excitons are

$`|n_A={\displaystyle \frac{1}{\sqrt{n!}}}\left[\widehat{A}_0^{}(0)\right]^n|0,`$ (59)
$`|m_B={\displaystyle \frac{1}{\sqrt{m!}}}\left[\widehat{B}_0^{}(0)\right]^m|0.`$ (60)
The coherent states can be formally defined as eigenvectors of exciton annihilation operators

$`\widehat{A}_0(0)|\alpha _A=\alpha |\alpha _A,`$ (61)
$`\widehat{B}_0(0)|\beta _B=\beta |\beta _B.`$ (62)
It is noticed that this definition can work well only in the low excitation case. This is because the generic expansion of a coherent state concerns the Fock states with higher excitations. However, the coherent state with smaller average exciton numbers can still approximately describe the quantum coherence natures of the exciton systems. Actually, we can understand the Fock state or the coherent state in terms of the single atomic states (for the details please see the appendix A). In this sense, the initial state of Frenkel exciton can be written according to the initial preparations of the single atomic states. A simplest illustration is that the vacuum state of Frenkel exciton is just the state formed by all atoms in the ground state. According to the appendix A, the Fock state of Frenkel exciton is a symmetrized many-atom state with certain atoms in the excited state. Especially, the coherent state of Frenkel exciton is an atomic ( SU(2)) coherent state.

The fluorescence of the two-mode excitons will be studied in the following of this section. If one take the initial state of total system as $`|\varphi =|\varphi _{ex}|\{0\}`$, i.e., all modes of the light field are initially in vacuum state, $`\widehat{E}_0^{(+)}(z,t)`$ in the positive frequency part of the field operator of Eq. (46) gives zero contribution. Therefore, in the following calculations, we can neglect safely the “quantum noise”, the terms proportional to $`\widehat{a}_q(0)`$, and investigate only the influence of various initial exciton states on the fluorescence of the low density excitations in the crystal slab.

We find that the light field emitted from the two-mode exciton system contains two eigenmodes: $`\mathrm{\Omega }_{01}`$ and $`\mathrm{\Omega }_{02}`$. These two modes of the light field can give an interference if there are non-diagonal terms in the light intensity. To see this, we need calculate the light intensity firstly. Substituting Eq. (46) and its Hermitian conjugate into Eq. (18), the light intensity is

$`I(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_1{\displaystyle \frac{\eta _1}{A}}\left|F_A(z,t)\right|^2\widehat{A}_0^{}(0)\widehat{A}_0(0)+{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_2{\displaystyle \frac{\eta _2}{A}}\left|F_B(z,t)\right|^2\widehat{B}_0^{}(0)\widehat{B}_0(0)`$ (64)
$`+{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_3{\displaystyle \frac{\eta _3}{A}}[F_A^{}(z,t)F_B(z,t)\widehat{A}_0^{}(0)\widehat{B}_0(0)+C.c.],`$
where $`\mathrm{}`$ denotes ensemble average, $`\mathrm{\Omega }_3=\sqrt{\mathrm{\Omega }_1\mathrm{\Omega }_2}`$, and $`\eta _3=\sqrt{\eta _1\eta _2}`$. The first two terms in Eq. (64) represent the average intensities of the two eigenmodes: $`\mathrm{\Omega }_{01}`$ and $`\mathrm{\Omega }_{02}`$, respectively. The last term proportional to $`\frac{1}{8}\mathrm{}\mathrm{\Omega }_3\frac{\eta _3}{A}\left[\mathrm{}\right]`$, however, gives a contribution of interference of the light fields.

When the exciton system is initially in a separable state $`\rho (0)=\rho _A\rho _B`$ and both $`\rho _A`$ and $`\rho _B`$ can be diagonal in Fock representation, there are only two terms in the light intensity, which represent the contributions from the two-mode excitons, respectively. There is no non-diagonal term in the light intensity due to $`\widehat{A}_0^{}(0)\widehat{B}_0(0)=0`$. The light intensity is

$`I(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_1{\displaystyle \frac{\eta _1}{A}}n_A[4e^{\eta _1(tz/c)}`$ (66)
$`+4\chi e^{\eta _2(tz/c)}]\mathrm{\Theta }(tz/c),`$
where $`\chi =\chi _0m_B/n_A`$, $`\chi _0=\mathrm{\Omega }_2\eta _2/\mathrm{\Omega }_1\eta _1`$, and $`n_A`$ ($`m_B`$) is the initial mean number of $`A_0`$($`B_0`$)-mode excitons. Note that $`\chi _0`$ is only determined by the intrinsic properties of the 2D sample. For a fixed $`\chi _0`$, $`\chi `$ can be used to describe the degree of unsymmetrically excitation within the crystal slab initially. Eq. (66) is plotted in Fig. 2, and our results show that $`I(z,t)`$ increases abruptly after the propagation time $`t=2\pi /\mathrm{\Omega }_1`$ then decays exponentially and no interference pattern appears in this case. When $`m_B=0`$, i.e., the $`B_0`$-mode excitons are initially in vacuum state, as shown in Fig. 2 (the solid line), our result will go back to the two-level crystal atoms case  due to the using of RWA and neglecting MPP in the derivation of Eq. (34) and Eq. (35), i.e., the mode of exciton being initially in the vacuum state gives no contribution to the dynamics of the total system. Besides, from the dot line ($`m_B/n_A=1`$) and dashed-dot line ($`m_B/n_A=5`$) of Fig. 2, we find that, with the increase of $`\chi `$, the amplitude of the light intensity becomes higher, which shows clearly that the amplitude of the light intensity is the sum of the contributions of the two eigenmodes of the radiation emitted from the two-mode excitons.

However if $`\rho _A`$ and $`\rho _B`$ are non-diagonal in Fock representation, e.g., the exciton system initially being prepared in a factorized coherent states $`|\varphi _{ex}=|\alpha _A|\beta _B`$, the light intensity becomes

$`I(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_1{\displaystyle \frac{\eta _1}{A}}n_A\mathrm{\Theta }(tz/c)\{4e^{\eta _1(tz/c)}+4\chi e^{\eta _2(tz/c)}`$ (68)
$`+8\sqrt{\chi }\mathrm{cos}[(\mathrm{\Omega }_{01}\mathrm{\Omega }_{02})(tz/c)+\varphi ]e^{(\eta _1+\eta _2)(tz/c)/2}\},`$
where $`\varphi `$ is the phase difference between the two coherent states. We find that the last term $`\frac{1}{8}\mathrm{}\mathrm{\Omega }_3\frac{\eta _3}{A}\left\{\mathrm{}\right\}`$ in Eq. (68) gives a nonzero contribution. The detector at $`z`$ may register an fluorescence signal oscillating at frequency $`|\mathrm{\Omega }_{01}\mathrm{\Omega }_{02}|`$. The temporal interference phenomenon may be observed by using a broadband detector  because with which there is no way to know the photon received by the detector is emitted from which modes of excitons. This unknowing for the “which-way” detection will result in the so called quantum beat phenomenon.

Form Eq. (68), one can easily find that $`I(z,t)`$ is oscillating with the time evolution. The detector at $`z`$ will receive the first peak after the propagation time $`t=z/c`$ (In fig. 3, we take $`z=2\pi c/\mathrm{\Omega }_1`$. It is not a necessary choice.). The influences of $`\chi `$ on the time evolution of the light intensity are investigated in Fig. 3, where we have used $`\varphi =0`$, i.e., zero phase difference between the two coherent states. Our results show that if one of the two-mode excitons is initially in vacuum state (say $`\beta =0`$), we recover the normal exponential decay of the light intensity. Besides, the amplitudes of the light intensity become larger with the increase of $`\chi `$. The result that the light intensity at initial stage does not go down to zero (the dot line and the dashed-dot line) is the consequence of unsymmetrically excitation, i.e., the initially generation of the two-mode excitons with different mean numbers.

All our mentioned above are separable state cases, in which the density matrix of the two-mode exciton system can be factorized initially. On the other hand, if the system is initially in an entangled state, e.g., $`|\varphi _{ex}=\frac{1}{\sqrt{2}}\left(|0_A|1_B+|1_A|0_B\right)`$, one may also observe the quantum beat phenomenon.

The first-order degree of coherence of the light field $`|g^{(1)}(\tau )|`$ as a function of the time delay $`\tau `$ is calculated by substituting Eq. (46) and its Hermitian conjugate into Eq. (19). We find that if one of the two-mode excitons is initially in vacuum state, e.g., $`m_B=0`$, then $`|g^{(1)}(\tau )|=1`$ with regardless of the another mode state, which implies that single-mode exciton in the 2D crystal slab emits a complete coherent light. Mathematically, the complete coherent light is obtained when $`|g^{(1)}(\tau )|=1`$. This happens when

$$\widehat{E}^{()}(z,t)\widehat{E}^{(+)}(z,t+\tau )^{}(z,t)(z,t+\tau ),$$
(69)
where $`(z,t)=\widehat{E}^{(+)}(z,t)`$, i.e., two-time correlation function of the total electric field may be factorized as the ensemble averages of the electric fields. As a special case, when there is only one mode exciton being excited initially, one can easily find that the factorized condition is satisfied. The physical meaning of $`g^{(1)}(\tau )`$ can be understood by considering the visibility of the interference fringes, which is proportional to $`|g^{(1)}(\tau )|`$. A maximum visibility of the fringes is obtained when $`|g^{(1)}(\tau )|=1`$.

If both of the two-mode excitons are populated initially with $`\rho _{ex}=\rho _A\rho _B`$ and both $`\rho _A`$ and $`\rho _B`$ being diagonal in Fock representation, the factorized condition of the two-time correlation function is broken, which leads to $`0<|g^{(1)}(\tau )|<1`$. The magnitude of the first-order degree of coherence for this case is

$$|g^{(1)}(z,t;z,t+\tau )|_{z=ct}=\frac{\sqrt{1+2\chi e^{(\eta _2\eta _1)\tau /2}\mathrm{cos}(\mathrm{\Omega }_{01}\mathrm{\Omega }_{02})\tau +\chi ^2e^{(\eta _2\eta _1)\tau }}}{\sqrt{1+\chi +\chi e^{(\eta _2\eta _1)\tau }+\chi ^2e^{(\eta _2\eta _1)\tau }}}.$$
(70)
We find that $`|g^{(1)}(\tau )|`$ oscillates regularly with $`\tau `$, as shown in Fig. 4(a). The magnitude of the first-order degree of coherence also depends on $`\chi `$: the larger the degree of unsymmetrically excitation is, the smaller the oscillating amplitude of $`|g^{(1)}(\tau )|`$ is. When $`\chi \mathrm{}`$, $`|g^{(1)}(\tau )|`$ tends to $`1`$. For the case of the excitons being initially in a factorized coherent state, as shown in Fig. 4(b), $`|g^{(1)}(\tau )|`$ is equal to $`1`$ due to the satisfaction of the factorized condition.

The second degree of coherence $`|g^{(2)}(\tau )|`$ as a function of $`\tau `$ is also calculated. We find that, when one of the two-mode excitons is initially in a vacuum state (say $`m_B=0`$), then: (1) $`|g^{(2)}(\tau )|=2`$, for the case of the other mode exciton being in a chaotic state; (2) $`|g^{(2)}(\tau )|=11/n`$, for the case of the $`A_0`$-mode excitons being initially in a Fock state; (3) $`|g^{(2)}(\tau )|=1`$, for the coherent state case. If both of the two-mode excitons are populated initially, we find that $`|g^{(2)}(\tau )|`$ oscillates regularly with $`\tau `$, and the oscillating amplitude of $`|g^{(2)}(\tau )|`$ becomes smaller with the increase of $`\chi `$.

## V Non-perturbation approach without WWA

In general, the rotating-wave approximation and Weisskopf-Wigner approximation are used frequently in quantum optics to study the fluorescence emitted from a single atom system. However, for the case of exciton radiance from a multi-atom system, we do not know exactly if they are valid because of the super-radiance feature of excitons. In this section, we will restudy the fluorescence of Frenkel excitons in the V-typed crystal slab by using non-perturbation approach  without the rotating-wave approximation in the interaction Hamiltonian $`(e/mc)𝐏𝐀`$. The self-interaction term $`(e^2/2mc^2)𝐀^2`$ is also included in our model to avoid un-physical roots in the characteristic equations for the exciton dispersion relation . Both the stimulated emission and reabsorption effects are taken into account in the new theoretical treatment. We focus our attention on the radiative decay rates of the excitons and the light intensities for various exciton initial states to compare with our previous results in Sec. IV.

The interaction Hamiltonian $`H_1`$ between the Frenkel excitons and photons without RWA has been given in Eq. (8) in section II. Our further discussion will also take the self-interaction term $`\frac{e^2}{2mc^2}𝐀^2`$ into account. In second quantization, this self-interaction Hamiltonian is written as

$$\widehat{H}_2=\mathrm{}\underset{q,q^{};l}{}f(q,q^{})\left[\widehat{a}_q+\widehat{a}_q^{}\right]\left[\widehat{a}_q^{}+\widehat{a}_q^{}^{}\right]e^{i(q+q^{})la},$$
(71)
in which

$$f(q,q^{})=\frac{N_L\pi e^2}{mcV\sqrt{|qq^{}|}}\underset{\lambda ,\lambda ^{}}{}𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}},$$
(72)
where $`𝐞_{q\lambda }`$ is the unit polarization vector of the $`(q,\lambda )`$-mode photon with $`\lambda =1,2`$. The sum $`_{\lambda ,\lambda ^{}}𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}`$ had been given for the two-level atom case  with certain assumptions. Following them we can give a similar result for the V-type three-level case (the detailed calculations are presented in Appendix B):

$`f(q,q^{})`$ $`=`$ $`{\displaystyle \frac{2\pi N_L}{c\mathrm{}V\sqrt{|qq^{}|}}}\left[\mathrm{\Omega }_1d_1^2+\mathrm{\Omega }_2d_2^2\right]`$ (73)
$`=`$ $`{\displaystyle \frac{1}{N}}\left[{\displaystyle \frac{1}{\mathrm{\Omega }_1}}G_1(q)G_1(q^{})+{\displaystyle \frac{1}{\mathrm{\Omega }_2}}G_2(q)G_2(q^{})\right].`$ (74)
By making use of

$$\underset{l}{}e^{i(q+q^{})la}=N\underset{k}{}O(q^{}k)O(k+q),$$
(75)
we can get the self-interaction Hamiltonian of photon $`\widehat{H}_2`$. The Heisenberg equations for exciton and photon operators can be obtained from the total interaction Hamiltonian $`\widehat{H}_{int}=\widehat{H}_1+\widehat{H}_2`$ (Eq. (C1) in Appendix C). Performing half-side Fourier transformation (HSFT)

$$A_q(\omega )=_0^{\mathrm{}}\widehat{A}_q(t)e^{i\omega t}𝑑t$$
(76)
on both sides of the equations, four algebraic equations about $`A_k(\omega )`$, $`B_k(\omega )`$, $`a_q(\omega )`$, and $`a_q^{}(\omega )`$ are obtained as listed in Appendix C. Combining these equations we get

$`a_q(\omega )a_q^{}(\omega )`$ (77)
$`=`$ $`{\displaystyle \frac{2\omega }{\mathrm{\Omega }_1(\omega ^2q^2c^2)}}G_1(q){\displaystyle \underset{k}{}}O(kq)`$ (82)
$`\times \left[\omega \left[A_k(\omega )A_k^{}(\omega )\right]i\left[A_k(0)A_k^{}(0)\right]\right]`$
$`+{\displaystyle \frac{2\omega }{\mathrm{\Omega }_2(\omega ^2q^2c^2)}}G_2(q){\displaystyle \underset{k}{}}O(kq)`$
$`\times \left[\omega \left[B_k(\omega )B_k^{}(\omega )\right]i\left[B_k(0)B_k^{}(0)\right]\right]`$
$`+{\displaystyle \frac{i}{\omega |q|c}}a_q(0){\displaystyle \frac{i}{\omega +|q|c}}a_q^{}(0),`$
for the electric field of the fluorescence. In order to calculate it, one need to know $`A_k(\omega )A_k^{}(\omega )`$ and $`B_k(\omega )B_k^{}(\omega )`$, which are determined by $`a_q(\omega )+a_q^{}(\omega )`$. In appendix C, two coupled equations for $`A_k(\omega )A_k^{}(\omega )`$ and $`B_k(\omega )B_k^{}(\omega )`$ are presented in detail. We solve these two equations consistently and replace them into Eq. (82). Then we get the electric field by the inverse HSFT

$$E(z,t)=\frac{1}{2\pi }_{\mathrm{}+iϵ}^{\mathrm{}+iϵ}E(z,\omega )e^{i\omega t}𝑑\omega ,$$
(83)
where

$`E(z,\omega )`$ $`=`$ $`i{\displaystyle \underset{q}{}}\sqrt{{\displaystyle \frac{2\pi |q|c\mathrm{}}{V}}}\left[a_q(\omega )a_q^{}(\omega )\right]e^{iqz}`$ (84)
$`=`$ $`i{\displaystyle _{\mathrm{}}^{\mathrm{}}}\sqrt{{\displaystyle \frac{|q|c\mathrm{}L}{2\pi A}}}\left[a_q(\omega )a_q^{}(\omega )\right]e^{iqz}𝑑q,`$ (85)
for sufficiently large $`L`$.

From the explicit form of the electric field operator, the radiative decay rates of Frenkel excitons can be solved. For the sake of simplicity, we will restrict ourselves to the monolayer case $`N=1`$, in which only zero wave vector of excitons are involved in our model. Solving the two coupled equations for $`A_0(\omega )A_0^{}(\omega )`$ and $`B_0(\omega )B_0^{}(\omega )`$ (see Eqs. (C55) and (C60) in Appendix C), we get two independent equations

$`\{[\omega ^2\mathrm{\Omega }_1^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )][\omega ^2\mathrm{\Omega }_2^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )]`$ (87)
$`{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(3)}(\omega ){\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(3)}(\omega )\left\}\right[A_0(\omega )A_0^{}(\omega )]`$
$`=`$ $`i\left[\omega ^2\mathrm{\Omega }_2^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )\right]\left[(\omega +\mathrm{\Omega }_1)A_0(0)(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$ (90)
$`+2i\omega F_{00}^{(3)}(\omega )\left[(\omega +\mathrm{\Omega }_2)B_0(0)(\omega \mathrm{\Omega }_2)B_0^{}(0)\right]`$
$`i{\displaystyle \frac{2\omega }{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )\left(\omega ^2\mathrm{\Omega }_2^2\right)\left[A_0(0)A_0^{}(0)\right],`$
and

$`\{[\omega ^2\mathrm{\Omega }_2^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )][\omega ^2\mathrm{\Omega }_1^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )]`$ (92)
$`{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(3)}(\omega ){\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(3)}(\omega )\left\}\right[B_0(\omega )B_0^{}(\omega )]`$
$`=`$ $`i\left[\omega ^2\mathrm{\Omega }_1^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )\right]\left[(\omega +\mathrm{\Omega }_2)B_0(0)(\omega \mathrm{\Omega }_2)B_0^{}(0)\right]`$ (95)
$`+2i\omega F_{00}^{(3)}(\omega )\left[(\omega +\mathrm{\Omega }_1)A_0(0)(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$
$`i{\displaystyle \frac{2\omega }{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )(\omega ^2\mathrm{\Omega }_1^2)\left[B_0(0)B_0^{}(0)\right],`$
where

$$F_{00}^{(j)}(\omega )=i\frac{\mathrm{\Omega }_j\eta _j}{2\omega }\text{}\eta _j=\frac{af_j^2}{2c}\text{}f_j^2=\frac{8\pi \mathrm{\Omega }_jd_j^2}{\mathrm{}a^3},$$
(96)
for $`j=1\text{}2\text{}3`$. $`F_{00}^{(j)}(\omega )`$ for $`j=1\text{}2`$ give the exciton wave function’s overlap for the $`A_0`$\- and $`B_0`$-mode exciton, respectively, however $`F_{00}^{(3)}(\omega )`$ is that between the two modes. Note that MPP has been taken into account in Eqs. (90) and (95). Then, we get the same characteristic equations for decay rates and frequency shifts of the two exciton modes $`A_0`$ and $`B_0`$

$$\zeta (\omega )=0,$$
(97)
where

$`\zeta (\omega )`$ $`=`$ $`(\omega ^2\mathrm{\Omega }_1^2+i\eta _1\omega )(\omega ^2\mathrm{\Omega }_2^2+i\eta _2\omega )+\eta _1\eta _2\omega ^2.`$ (98)
The roots of Eq. (97) can be solved exactly,

$`\omega _1`$ $`=`$ $`\mathrm{\Omega }_{01}i\mathrm{\Gamma }_1/2,\omega _2=\mathrm{\Omega }_{01}i\mathrm{\Gamma }_1/2,`$ (100)
$`\omega _3`$ $`=`$ $`\mathrm{\Omega }_{02}i\mathrm{\Gamma }_2/2,\omega _4=\mathrm{\Omega }_{02}i\mathrm{\Gamma }_2/2,`$ (101)
which determine the poles of $`A_0(\omega )A_0^{}(\omega )`$ and $`B_0(\omega )B_0^{}(\omega )`$. The imaginary parts of the roots are the decay rates of the excitons, and the real parts are the renormalized physical frequencies. For the case $`\mathrm{\Omega }_1=\mathrm{\Omega }_2=\mathrm{\Omega }`$, i.e., the degenerate case, we get

$`\omega _1`$ $`=`$ $`\mathrm{\Omega },\omega _3=\mathrm{\Omega }_0i\left(\eta _1+\eta _2\right)/2,`$ (102)
$`\omega _2`$ $`=`$ $`\mathrm{\Omega },\omega _4=\mathrm{\Omega }_0i\left(\eta _1+\eta _2\right)/2,`$ (103)
where

$$\mathrm{\Omega }_0=\sqrt{\mathrm{\Omega }^2\frac{1}{4}\left(\eta _1+\eta _2\right)^2}\mathrm{\Omega }\left[1\frac{\left(\eta _1+\eta _2\right)^2}{8\mathrm{\Omega }^2}\right].$$
(104)
For the non-degenerate case $`\mathrm{\Omega }_1\mathrm{\Omega }_2`$, we assume that the physical roots of the characteristic equation are not far away from $`\pm \mathrm{\Omega }_j`$ due to the condition $`\eta _j\mathrm{\Omega }_j`$. Therefore, one can expand $`\omega `$ up to third-order of $`\eta _j/\mathrm{\Omega }_j`$, and get

$`\mathrm{\Omega }_{01}`$ $`=`$ $`\mathrm{\Omega }_1\left[1{\displaystyle \frac{\eta _1^2}{8\mathrm{\Omega }_1^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{\eta _1\eta _2}{\mathrm{\Omega }_1^2\mathrm{\Omega }_2^2}}\right],`$ (105)
$`\mathrm{\Omega }_{02}`$ $`=`$ $`\mathrm{\Omega }_2\left[1{\displaystyle \frac{\eta _2^2}{8\mathrm{\Omega }_2^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{\eta _1\eta _2}{\mathrm{\Omega }_2^2\mathrm{\Omega }_1^2}}\right],`$ (106)
$`\mathrm{\Gamma }_1`$ $`=`$ $`\eta _1\left[1{\displaystyle \frac{\eta _2\mathrm{\Omega }_1^2\eta _1\mathrm{\Omega }_2^2}{\left(\mathrm{\Omega }_1^2\mathrm{\Omega }_2^2\right)^2}}\eta _2\right],`$ (107)
$`\mathrm{\Gamma }_2`$ $`=`$ $`\eta _2\left[1{\displaystyle \frac{\eta _1\mathrm{\Omega }_2^2\eta _2\mathrm{\Omega }_1^2}{\left(\mathrm{\Omega }_2^2\mathrm{\Omega }_1^2\right)^2}}\eta _1\right],`$ (108)
where $`\eta _j`$ with $`j=1,2`$ are given in Eqs. (36) and (37).

We compare our results with the two-level lattice atoms case , in which the radiative decay rate of excitons is $`\eta =3\left(\frac{\pi \lambda ^2}{a^2}\right)\gamma `$, where $`\gamma `$ is the decay rate of an isolated lattice atom, $`\lambda =c/\mathrm{\Omega }`$ denotes the reduced wavelength, and physical frequency is $`\mathrm{\Omega }_0=\mathrm{\Omega }\left[1\frac{\eta ^2}{8\mathrm{\Omega }^2}\right]`$. However, we see from Eq. (108), for the three-level lattice atom case, the decay rates $`\mathrm{\Gamma }_j`$ are different from the two-level atom case $`\eta _j`$ by a third-order correction, which comes from MPP (including stimulation emission, absorption and other high-order processes). Here $`\mathrm{\Omega }_{01}\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_{02}\mathrm{\Omega }_2`$ denote the frequency shifts for the $`A_0`$-mode and $`B_0`$-mode exciton, respectively. For the two-level case the frequency shift is $`\frac{\eta ^2}{8\mathrm{\Omega }^2}`$. Whereas, for three-level case considered here, the phase shifts for the two-mode exciton are $`\frac{\eta _i^2}{8\mathrm{\Omega }_i^2}+\frac{1}{2}\frac{\eta _i\eta _j}{\mathrm{\Omega }_i^2\mathrm{\Omega }_j^2}`$, $`i\text{}j=1\text{}2`$. The last term in the frequency shifts is a second-order correction, which can be adjusted by tuning the energy spacing of the upper two levels of the lattice atoms.

In the remainder of this paper, we will present the explicit expressions of the electric field operator. The time evolution and spatial distribution of the light intensity are calculated to compare with the results of RWA in Sec. IV.

From Eqs. (90) and (95), the explicit expressions for $`A_0(\omega )A_0^{}(\omega )`$ (Eq. (C64)) and $`B_0(\omega )B_0^{}(\omega )`$ (Eq. (C68)) are obtained. Replacing them into Eq. (82), we get the explicit light field operator in terms of the initial exciton field operators

$`a_q(\omega )a_q^{}(\omega )`$ (109)
$`=`$ $`{\displaystyle \frac{2i\omega G_1(q)(\omega ^2\mathrm{\Omega }_2^2)}{(\omega ^2q^2c^2)\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_1)A_0(0)+(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$ (110)
$`+`$ $`{\displaystyle \frac{2i\omega G_2(q)(\omega ^2\mathrm{\Omega }_1^2)}{(\omega ^2q^2c^2)\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_2)B_0(0)+(\omega \mathrm{\Omega }_2)B_0^{}(0)\right],`$ (111)
where we have neglected the quantum noise terms proportional to $`a_q(0)`$ and $`a_q^{}(0)`$ in the above calculation. Substituting $`a_q(\omega )a_q^{}(\omega )`$ into Eq. (85) and choosing a proper contour integration in the upper complex $`q`$ plane, we get $`E(z,\omega )`$ in the positive $`z`$ region outside the crystal slab as the following,

$`E(z,\omega )`$ (113)
$`=`$ $`i\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_1\eta _1}{cA}}}{\displaystyle \frac{\omega ^2\mathrm{\Omega }_2^2}{\zeta (\omega )}}[(\omega +\mathrm{\Omega }_1)A_0(0)`$ (115)
$`+(\omega \mathrm{\Omega }_1)A_0^{}(0)]e^{i\frac{\omega }{c}z}`$
$`+`$ $`i\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_2\eta _2}{cA}}}{\displaystyle \frac{\omega ^2\mathrm{\Omega }_1^2}{\zeta (\omega )}}[(\omega +\mathrm{\Omega }_2)B_0(0)`$ (117)
$`+(\omega \mathrm{\Omega }_2)B_0^{}(0)]e^{i\frac{\omega }{c}z},`$
where $`A`$ is the area of the layer. As mentioned above, we choose the normalization volume for the photon as $`AL`$.

The electric field $`E(z,t)`$ in the positive $`z`$ region outside the crystal slab can be calculated by Eq. (83). For the case $`zct>0`$ ($`<0`$), the contour integral over $`d\omega `$ can be performed by choosing the integration path in the upper (lower) complex $`\omega `$ plane. The result is that $`E(z,t)=0`$, for $`zct>0`$; but for $`zct<0`$,

$$E(z,t)=^{(+)}(z,t)+\mathrm{H}.\mathrm{c}..$$
(119)
The explicit expression for $`^{(+)}(z,t)`$ in Appendix C shows that, in positive $`z`$ region, the electric field generated by the two-mode exciton is damped exponentially with two eigen-decay rates $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$. The corresponding eigen-modes are linear combinations of $`A_0`$\- and $`B_0`$-mode. The electric field in negative $`z`$ region can also be derived similarly and shows that it is a damped wave propagating in the negative $`z`$ direction. Since all the roots of the characteristic equation are in lower half plane of complex $`\omega `$ plane, then $`E(z,t)=0`$, for $`zct>0`$, so our result is reasonable and obeys the causal rule.

Similarly, we can also calculate the positive frequency part of the electric field operator, which is determined by $`a_q(\omega )`$. The explicit expression of $`a_q(\omega )`$ can be obtained by substituting Eqs. (C64) and (C68) into Eq. (C16) (see details in Appendix C.) and omitting quantum noise terms. We get

$`(\omega |q|c)a_q(\omega )`$ (120)
$`=`$ $`iG_1(q){\displaystyle \frac{\omega ^2\mathrm{\Omega }_2^2}{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_1)A_0(0)+(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$ (121)
$`+`$ $`iG_2(q){\displaystyle \frac{\omega ^2\mathrm{\Omega }_1^2}{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_2)B_0(0)+(\omega \mathrm{\Omega }_2)B_0^{}(0)\right].`$ (122)
Therefore, we get

$`E^{(+)}(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q{\displaystyle _{\mathrm{}+iϵ}^{\mathrm{}+iϵ}}𝑑\omega \sqrt{{\displaystyle \frac{|q|c\mathrm{}L}{2\pi A}}}a_q(\omega )e^{i(qz\omega t)}`$ (124)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑qe^{iqz}{\displaystyle _{\mathrm{}+iϵ}^{\mathrm{}+iϵ}}𝑑\omega {\displaystyle \frac{e^{i\omega t}}{\omega _q\omega }}`$ (129)
$`\times \{\sqrt{{\displaystyle \frac{c\mathrm{}\mathrm{\Omega }_1\eta _1}{4\pi A}}}{\displaystyle \frac{(\omega ^2\mathrm{\Omega }_2^2)}{\zeta (\omega )}}[(\omega +\mathrm{\Omega }_1)`$
$`\times A_0(0)+(\omega \mathrm{\Omega }_1)A_0^{}(0)]`$
$`+\sqrt{{\displaystyle \frac{c\mathrm{}\mathrm{\Omega }_2\eta _2}{4\pi A}}}{\displaystyle \frac{(\omega ^2\mathrm{\Omega }_1^2)}{\zeta (\omega )}}[(\omega +\mathrm{\Omega }_2)`$
$`\times B_0(0)+(\omega \mathrm{\Omega }_2)B_0^{}(0)]\}.`$
A straightforward calculations give the following results:

$$E^{(+)}(z,t)=0\text{, for }\text{t¡0}\text{,}$$
(130)
$`E^{(+)}(z,t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_1\eta _1}{4Ac}}}[F_A^{(+)}(z,t)A_0(0)`$ (132)
$`+F_A^{()}(z,t)A_0^{}(0)]`$
$`+`$ $`\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_2\eta _2}{4Ac}}}[F_B^{(+)}(z,t)B_0(0)`$ (134)
$`+F_B^{()}(z,t)B_0^{}(0)]\text{, for }t>0\text{,}`$
in which the time-dependent coefficients $`F_A^{(\pm )}(z,t)`$ and $`F_B^{(\pm )}(z,t)`$ are

$`F_A^{(\pm )}(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}d\omega _q2[{\displaystyle \frac{\left(\omega _q^2\mathrm{\Omega }_2^2\right)\left(\omega _q\pm \mathrm{\Omega }_1\right)e^{i\omega _qt}}{\left(\omega _q\omega _1\right)\left(\omega _q\omega _2\right)\left(\omega _q\omega _3\right)\left(\omega _q\omega _4\right)}}`$ (139)
$`{\displaystyle \frac{\left(\omega _1^2\mathrm{\Omega }_2^2\right)\left(\omega _1\pm \mathrm{\Omega }_1\right)e^{i\omega _1t}}{\left(\omega _q\omega _1\right)\left(\omega _1\omega _2\right)\left(\omega _1\omega _3\right)\left(\omega _1\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _2^2\mathrm{\Omega }_2^2\right)\left(\omega _2\pm \mathrm{\Omega }_1\right)e^{i\omega _2t}}{\left(\omega _q\omega _2\right)\left(\omega _2\omega _1\right)\left(\omega _2\omega _3\right)\left(\omega _2\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _3^2\mathrm{\Omega }_2^2\right)\left(\omega _3\pm \mathrm{\Omega }_1\right)e^{i\omega _3t}}{\left(\omega _q\omega _3\right)\left(\omega _3\omega _1\right)\left(\omega _3\omega _2\right)\left(\omega _3\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _4^2\mathrm{\Omega }_2^2\right)\left(\omega _4\pm \mathrm{\Omega }_1\right)e^{i\omega _4t}}{\left(\omega _q\omega _4\right)\left(\omega _4\omega _1\right)\left(\omega _4\omega _2\right)\left(\omega _4\omega _3\right)}}\left]\mathrm{cos}\right(\omega _q{\displaystyle \frac{z}{c}}),`$
$`F_B^{(\pm )}(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}d\omega _q2[{\displaystyle \frac{\left(\omega _q^2\mathrm{\Omega }_1^2\right)\left(\omega _q\pm \mathrm{\Omega }_2\right)e^{i\omega _qt}}{\left(\omega _q\omega _1\right)\left(\omega _q\omega _2\right)\left(\omega _q\omega _3\right)\left(\omega _q\omega _4\right)}}`$ (144)
$`{\displaystyle \frac{\left(\omega _1^2\mathrm{\Omega }_1^2\right)\left(\omega _1\pm \mathrm{\Omega }_2\right)e^{i\omega _1t}}{\left(\omega _q\omega _1\right)\left(\omega _1\omega _2\right)\left(\omega _1\omega _3\right)\left(\omega _1\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _2^2\mathrm{\Omega }_1^2\right)\left(\omega _2\pm \mathrm{\Omega }_2\right)e^{i\omega _2t}}{\left(\omega _q\omega _2\right)\left(\omega _2\omega _1\right)\left(\omega _2\omega _3\right)\left(\omega _2\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _3^2\mathrm{\Omega }_1^2\right)\left(\omega _3\pm \mathrm{\Omega }_2\right)e^{i\omega _3t}}{\left(\omega _q\omega _3\right)\left(\omega _3\omega _1\right)\left(\omega _3\omega _2\right)\left(\omega _3\omega _4\right)}}`$
$`{\displaystyle \frac{\left(\omega _4^2\mathrm{\Omega }_1^2\right)\left(\omega _4\pm \mathrm{\Omega }_2\right)e^{i\omega _4t}}{\left(\omega _q\omega _4\right)\left(\omega _4\omega _1\right)\left(\omega _4\omega _2\right)\left(\omega _4\omega _3\right)}}\left]\mathrm{cos}\right(\omega _q{\displaystyle \frac{z}{c}}),`$
where $`\omega _j`$, for $`j=1,2,3,4`$, are four roots of the characteristic equation of Eq. (97). We again need to carry out the integration over $`d\omega _q`$ in the above two equations. It is easy to prove that $`F_A^{(+)}(z,t)`$ can be simplified as

$`F_A^{(+)}(z,t)`$ $`=`$ $`{\displaystyle \frac{i}{\pi }}{\displaystyle _0^{\mathrm{}}}d\omega _q2\mathrm{cos}(\omega _q{\displaystyle \frac{z}{c}})[c_1{\displaystyle \frac{e^{i\omega _qt}e^{i\omega _1t}}{\omega _q\omega _1}}`$ (147)
$`+c_2{\displaystyle \frac{e^{i\omega _qt}e^{i\omega _2t}}{\omega _q\omega _2}}+c_3{\displaystyle \frac{e^{i\omega _qt}e^{i\omega _3t}}{\omega _q\omega _3}}`$
$`+c_4{\displaystyle \frac{e^{i\omega _qt}e^{i\omega _4t}}{\omega _q\omega _4}}],`$
where the four coefficients are

$`c_1`$ $`=`$ $`{\displaystyle \frac{\left(\omega _1^2\mathrm{\Omega }_2^2\right)\left(\omega _1+\mathrm{\Omega }_1\right)}{\left(\omega _1\omega _2\right)\left(\omega _1\omega _3\right)\left(\omega _1\omega _4\right)}},`$ (148)
$`c_2`$ $`=`$ $`{\displaystyle \frac{\left(\omega _2^2\mathrm{\Omega }_2^2\right)\left(\omega _2+\mathrm{\Omega }_1\right)}{\left(\omega _2\omega _1\right)\left(\omega _2\omega _3\right)\left(\omega _2\omega _4\right)}},`$ (149)
$`c_3`$ $`=`$ $`{\displaystyle \frac{\left(\omega _3^2\mathrm{\Omega }_2^2\right)\left(\omega _3+\mathrm{\Omega }_1\right)}{\left(\omega _3\omega _1\right)\left(\omega _3\omega _2\right)\left(\omega _3\omega _4\right)}},`$ (150)
$`c_4`$ $`=`$ $`{\displaystyle \frac{\left(\omega _4^2\mathrm{\Omega }_2^2\right)\left(\omega _4+\mathrm{\Omega }_1\right)}{\left(\omega _4\omega _1\right)\left(\omega _4\omega _2\right)\left(\omega _4\omega _3\right)}}.`$ (151)
The above integral over $`d\omega _q`$ can be done as Sec. IV. The explicit expressions of $`F_A^{()}(z,t)`$ and $`F_B^{(\pm )}(z,t)`$ can also be obtained with the same calculations. Our results also confirm the causal rule.

For an arbitrary exciton initial state $`|\varphi _{ex}`$, the light intensity is defined by Eq. (18). When the excitons are initially in a state with density matrix $`\rho (0)=\rho _A\rho _B`$, and both $`\rho _A`$ and $`\rho _B`$ diagonal in Fock representation, the light intensity is

$`I(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_1{\displaystyle \frac{\eta _1}{A}}n_A\left[\left|F_A^{(+)}(z,t)\right|^2+\left|F_A^{()}(z,t)\right|^2\right]`$ (152)
$`+`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_2{\displaystyle \frac{\eta _2}{A}}m_B\left[\left|F_B^{(+)}(z,t)\right|^2+\left|F_B^{()}(z,t)\right|^2\right],`$ (153)
which is shown in Fig. 5. When $`m_B=0`$, the solid line, our result will go back to the two-level lattice atoms case . From the dot line ($`m_B/n_A=1`$) and dashed-dot line ($`m_B/n_A=5`$) of Fig. 5, we find that with the increase of $`\chi `$, the amplitude of the light intensity becomes higher, which is the same with that of obtained by using RWA (Fig. 2). Contrary to the results of Sec. IV, however, when we consider it without RWA, the light intensity does not decay exponentially but in an irregular way due to the existence of counter-rotating terms in $`E^{(+)}(z,t)`$. Besides, it’s also deserved to mentioned that the contributions of the non-rotating terms do not appear as quivers presented in Ref. .

For the case that the two-mode excitons are initially in a factorized coherent states, the light intensity becomes

$`I(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_1{\displaystyle \frac{\eta _1}{A}}\left[\right|\alpha |^2|F_A^{(+)}(z,t)|^2+|\alpha |^2|F_A^{()}(z,t)|^2`$ (160)
$`+(\alpha ^{})^2F_A^{(+)}(z,t)F_A^{()}(z,t)+\alpha ^2F_A^{()}(z,t)F_A^{(+)}(z,t)]`$
$`+{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_2{\displaystyle \frac{\eta _2}{A}}\left[\right|\beta |^2|F_B^{(+)}(z,t)|^2+|\beta |^2|F_B^{()}(z,t)|^2`$
$`+\left(\beta ^{}\right)^2F_B^{(+)}(z,t)F_B^{()}(z,t)+\beta ^2F_B^{()}(z,t)F_B^{(+)}(z,t)]`$
$`+{\displaystyle \frac{1}{8}}\mathrm{}\mathrm{\Omega }_3{\displaystyle \frac{\eta _3}{A}}[\alpha ^{}\beta F_A^{(+)}(z,t)F_B^{(+)}(z,t)+\alpha ^{}\beta ^{}F_A^{(+)}(z,t)F_B^{()}(z,t)`$
$`+\alpha \beta F_A^{()}(z,t)F_B^{(+)}(z,t)+\alpha \beta ^{}F_A^{()}(z,t)F_B^{()}(z,t)+C.c.],`$
where the last term $`\frac{1}{8}\mathrm{}\mathrm{\Omega }_3\frac{\eta _3}{A}[\mathrm{}]`$ gives a temporal interference term. The effects of phase difference between the two coherent states on the light intensity are studied in Fig. 6. We find that the first peak (at $`t=z/c`$) of the light intensity becomes lower in amplitude with the increase of the phase difference from $`\varphi =0`$ to $`\varphi =\pi `$ (monotonic regime). In fact, the whole curve will be shifted left with the increase of the phase difference within the monotonic regime, which leads to the magnitude of the peak (at $`t=z/c`$) becomes lower (comparing Fig. 6(b) and Fig. 6(c) with Fig. 6(a)). When $`\varphi =\pi `$, $`I(z,t=z/c)`$ tends to zero, i.e., the first peak is disappeared (see Fig. 6(c)). Our results also show that both the phase difference and the degree of unsymmetrically excitation do not effect the oscillation frequency of the light intensity, which is determined by the exciton splitting. The oscillation behavior in the light intensity may take place as long as $`|\mathrm{\Omega }_{01}\mathrm{\Omega }_{02}|>\mathrm{\Gamma }_j`$, i.e., the exciton splitting is greater than the natural linewidth of exiton, so that one can observe the beating phenomenon within the lifetime of the exciton.

The spatial distribution of the light intensity is plotted in Fig. 7 for the case that the two-mode excitons are initially in a factorized Fock state. Our results show that $`I(z,t_0)`$ increases exponentially (see Eq. (66)) within the regions $`0<|z|<ct_0`$ for the case of rotating-wave approximation, and vanishes immediately as $`|z|`$ goes beyond $`ct_0`$. The solid lines in Fig. 7 are obtained without RWA and shows small-amplitude oscillations due to the contribution of counter-rotating terms. Compared with the results of Ref. , our results show that not only the electric field $`E(z,t)`$ but also the light intensity $`I(z,t)`$ do meet the requirement of causal rule.

## VI conclusion

In summary, we have studied the collective radiations of a collection of many V-type three-level atoms in a crystal slab. By introducing two-mode exciton operators in the large $`N`$ limits of the collective quasi-spin operators, these coherent radiations can be depicted as the fluorescence of low density Frenkel excitons. The exciton fluorescence exhibits the stronger coherence natures that the statistical characters of spectrum are identical from the initial to final stages. This is indeed different from the ensemble situation with free-moving atoms that the atoms need a finitely-long time to produce a cooperative radiation, the enhanced fluorescence.

As a main result of this paper, the occurrence of the quantum beat aroused from which quantum states of Frenkel exciton is investigated in this paper. Our results show that not all the quantum states of excitons may lead to the oscillating behavior in the light field radiated from the two-mode excitons. This quantum interference phenomenon may be observed when the two-mode excitons are initially in a factorized coherent state or an entangled state. We expect that our theoretical study of quantum beat would be helpful in practical experiment to measure quantum states of excitons. Further study of the anharmonic exciton-exciton interaction in the model is needed. It is also pointed out that, the algebraic consideration for the definition of the multi-mode excitons can be generalized to study other exciton system, e.g., the quasi-spin wave collective excitations of $`\mathrm{\Lambda }`$-type atom collection in lattice of crystal that can be used as a new type quantum memory.

###### Acknowledgements.

This work is supported by the NSF of China and the knowledged Innovation Program(KIP) of the Chinese Academy of Science. It is also founded by the National Fundamental Research Program of China with No 001GB309310. Y. X. Liu is supported by Japan Society for the Promotion of Science (JSPS). The authors would like to thank C. Q. Cao, W. M. Liu, D. L. Zhou and S. X. Yu for valuable discussions.

## A SU(3) algebra structure of exciton operators for many-atom system

From a point of view based on the representation of Lie algebra, this section describes the mathematical origin of definition of exciton operator. Physically, this description will clarify why the conception of the exciton based on the collective operators of many atoms can be valid only in the case of low excitation. We first discuss SU(2) algebra structure of excitonic operators for the two-level many-atom system since the SU(3) case for the 2D crystal slab containing V-type three-level atoms shares the same basic idea as that of SU(2). The detailed discussions for the SU(3) case will follow that of SU(2) in this appendix.

We consider an ensemble of $`N`$ two-level atoms with their ground states $`|g^j`$ and the excited ones $`|e^j`$. Since we can define the quasi-spin with the Pauli operators

$$\sigma _{}^i=(\sigma _+^i)^{}=|g^{ii}e|\text{}\sigma _z^i=|e^{ii}e||g^{ii}g|,$$
(A1)
the total quasi-angular momentum operators

$$\widehat{J}_{}=\underset{i=1}{\overset{N}{}}\sigma _{}^i,\widehat{J}_+=(\widehat{J}_{})^{},\widehat{J}_z=\frac{1}{2}\underset{i=1}{\overset{N}{}}\sigma _z^i,$$
(A2)
define a representation with the highest weight $`J=N/2`$. The $`(2J+1)`$-dimensional irreducible spinor representation of SU(2) in the symmetric subspace is embedded in the total Hilbert space of dimension $`2^N`$.

For a physics system with dynamical SU(2) symmetry, its Hamiltonian $`\widehat{H}=\widehat{H}(\widehat{J}_{},\widehat{J}_+,\widehat{J}_z)`$ is a functional of $`\widehat{J}_\pm `$ and $`\widehat{J}_z`$. Since the Casimir operator $`\widehat{J}^2`$ commutes with $`\widehat{J}_\pm `$ and $`\widehat{J}_z`$, the eigenvalue $`J(J+1)`$ will keep conservation in the time evolution. Here, $`J`$ can take one of the integers and half-integers $`\frac{N}{2},\frac{N}{2}1,\mathrm{},0`$. In a physical process, which one of these $`J`$ takes depends on the initial state of the atomic ensemble. The symmetric state $`|J=\frac{N}{2},M=\frac{N}{2}=_{i=1}^N|g^i`$ represents the “condensate” with all atoms filling in the ground state. It is very similar to the situation of an electronic system that the filling in ground state forms the Fermi surface. In this sense, we can introduce the atomic collective excitation operators

$$\widehat{B}^{}=\frac{1}{\sqrt{N}}\widehat{J}_+,\widehat{B}=\frac{1}{\sqrt{N}}\widehat{J}_{},$$
(A3)
which is very similar to the exciton operators for an electron-hole pair. Considering that

$$\widehat{J}_z=\frac{1}{2}\left[1\sqrt{(N+1)^24\widehat{J}_+\widehat{J}_{}}\right],$$
(A4)
as one solution for the basic angular momentum relation

$`\widehat{J}_x^2+\widehat{J}_y^2+\widehat{J}_z^2={\displaystyle \frac{N}{2}}\left({\displaystyle \frac{N}{2}}+1\right),`$ (A5)
$`\widehat{J}_+\widehat{J}_{}=\widehat{J}_x^2+\widehat{J}_y^2+\widehat{J}_z={\displaystyle \frac{N}{2}}\left({\displaystyle \frac{N}{2}}+1\right)\widehat{J}_z^2+\widehat{J}_z,`$ (A6)
the commutation relation for atomic collective excitation operators

$$\left[\widehat{B}\text{}\widehat{B}^{}\right]=\frac{2}{N}\widehat{J}_z,$$
(A7)
will be reduced to the bosonic relation $`[\widehat{B},\widehat{B}^{}]=1`$ for those angular momentum state $`|J,M`$ with much smaller $`M`$ in comparison to $`N=2J`$. That is to say, in the limit with very large $`N`$ and low excitation, the collective excitation behave as a boson and thus we call it atomic exciton. The detailed proof for this was given in ref.  by considering the physical realization of q-deformed boson algebra  for a very large, but finite $`N`$. The main point to prove that is to consider that

$$\frac{\widehat{J}_z}{N}=\frac{1}{2}\left[\frac{1}{N}\sqrt{\left(\frac{1}{N}+1\right)^2\frac{4}{N^2}\widehat{J}_+\widehat{J}_{}}\right],$$
(A8)
approaches $`1/2`$ for the infinite $`N`$ and the low excitation. Since the operators can make sense by acting on the symmetric space, only for those low excitation states $`|J,M`$ with very small $`M`$ can let the value of $`\widehat{J}_+\widehat{J}_{}/N^2`$ approach zero, so that the bosonic commutation relation is obtained. We can also prove that

$`\widehat{B}^{}\widehat{B}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}|e^{ii}e|+O\left({\displaystyle \frac{1}{N}}\right),`$ (A9)
which means that the free part of the many-atom Hamiltonian can rationally be described as a free boson in the large $`N`$ and the low excitation.

Because the “condensate” in a ground state plays the crucial role in defining the exciton operators, the introduction of exciton operators has to depend on the configuration of the atoms and thus on the form of interaction. This is just the line, along which we will define the exciton operators for the V-type atomic system with the interaction Hamiltonian

$`\widehat{H}_I`$ $`=`$ $`g_1\widehat{a}^{}{\displaystyle \underset{i=1}{\overset{N}{}}}|g^{ii}e_1|+g_2\widehat{a}^{}{\displaystyle \underset{i=1}{\overset{N}{}}}|g^{ii}e_2|+H.c..`$ (A10)
It is easy to prove that

$`E_1={\displaystyle \underset{i=1}{\overset{N}{}}}|e_1^{ii}g|,F_1={\displaystyle \underset{i=1}{\overset{N}{}}}|g^{ii}e_1|,`$ (A11)
$`E_2={\displaystyle \underset{i=1}{\overset{N}{}}}|g^{ii}e_2|,F_2={\displaystyle \underset{i=1}{\overset{N}{}}}|e_2^{ii}g|,`$ (A12)
$`H_1={\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(|e_1^{ii}e_1||g^{ii}g|\right),`$ (A13)
$`H_2={\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(|g^{ii}g||e_2^{ii}e_2|\right),`$ (A14)
generate a SU(3) algebra with the Cartan subalgebra spanned by $`H_1`$ and $`H_2`$. The basic commutation relations are

$`[H_1,E_1]`$ $`=`$ $`E_1,[H_1,F_1]=F_1,`$ (A15)
$`[H_2,E_2]`$ $`=`$ $`E_2,[H_2,F_2]=F_2,`$ (A16)
$`[H_2,E_1]`$ $`=`$ $`{\displaystyle \frac{E_1}{2}},[H_2,F_1]={\displaystyle \frac{F_1}{2}},`$ (A17)
$`[H_1,E_2]`$ $`=`$ $`{\displaystyle \frac{E_2}{2}},[H_1,F_2]={\displaystyle \frac{F_2}{2}}.`$ (A18)
Two sets of the operators {$`E_1`$, $`F_1`$, $`H_1`$} and {$`E_2`$, $`F_2`$ ,$`H_2`$} generate two non-commutation SU(2) subalgebras, respectively. Actually, the above four collective operators define a spinor realization for the symmetric representation of SU(3) in $`N_L+1`$ dimensional space. The atom number $`N_L`$ determines the dimensions $`N_L+1`$ of representations. In this sense, we can understand the two-mode excitation in terms of the large $`N_L`$ limit of representations of SU(3), which just corresponds to the low density excitation region. Thus we can define the collective operators

$`\widehat{A}={\displaystyle \frac{1}{\sqrt{N}}}F_1\text{}\widehat{A}^{}={\displaystyle \frac{1}{\sqrt{N}}}E_1,`$ (A19)
$`\widehat{B}={\displaystyle \frac{1}{\sqrt{N}}}E_2\text{}\widehat{B}^{}={\displaystyle \frac{1}{\sqrt{N}}}F_2.`$ (A20)
For a very large $`N`$ and low excitation, it is easy to prove that $`\widehat{A}`$ and $`\widehat{B}`$ commute with each other and obey the standard bosonic commutation relation.

Now we can consider the construction of Frenkel excitonic states. In this way, the initial conditions for Frenkel exciton can be given in terms of the single atom preparations. For example, when $`n=0`$, $`1`$, $`2`$, …, the Fock states of the $`A`$-mode excitons

$$|n_A=\frac{1}{\sqrt{n!}}\left(\widehat{A}^{}\right)^n|0=\frac{1}{N^{n/2}\sqrt{n!}}\left(\underset{i=1}{\overset{N}{}}|e_1^{ii}g|\right)^n|0$$
(A21)
take the symmetric excitation states

$`|0=|g,g,\mathrm{}g,`$ (A22)
$`|1_A={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{j=1}{\overset{N}{}}}|g,g,\mathrm{},e_1^j,\mathrm{},g,`$ (A23)
$`|2_A={\displaystyle \frac{1}{N\sqrt{2}}}{\displaystyle \underset{j,k=1}{\overset{N}{}}}|g,g,\mathrm{},e_1^j,..,e_1^k,\mathrm{},g,`$ (A24)
$`\mathrm{}.`$ (A25)
The second example is the coherent state of the $`A`$-mode Frenkel exciton

$`|\alpha _A`$ $``$ $`\mathrm{exp}(\alpha \widehat{A}^{})|0`$ (A26)
$`=`$ $`{\displaystyle \underset{j=1}{}}\left[\mathrm{cos}\theta |g^j+\mathrm{sin}\theta e^{i\varphi }|e_1^j\right],`$ (A27)
where $`\mathrm{tan}\theta =\frac{|\alpha |}{\sqrt{N}}`$ and $`\alpha =|\alpha |e^{i\varphi }`$. The coherent nature of this many-atomic state is reflected by the fact that both $`\theta `$ and $`\varphi `$ are independent of the index $`j`$ of atoms. The quantum states of $`B`$-mode excitons can also be constructed by using the same procedure as discussed above.

## B The self-interaction term of the light field

In this section we will calculate $`𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}`$ for the three-level case. For an atom with a complete set of eigenvectors $`\{|n\}`$, we have

$`𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}`$ $`=`$ $`g|𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}|g`$ (B1)
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{n}{}}g|𝐞_{q\lambda }𝐱|nn|𝐩𝐞_{q^{}\lambda ^{}}|g`$ (B3)
$`{\displaystyle \underset{n}{}}g|𝐞_{q\lambda }𝐩|nn|𝐱𝐞_{q^{}\lambda ^{}}|g`$
$`=`$ $`{\displaystyle \frac{2m}{\mathrm{}e^2}}{\displaystyle \underset{n}{}}\mathrm{\Omega }_ng|𝐞_{q\lambda }𝐝|nn|𝐝𝐞_{q^{}\lambda ^{}}|g,`$ (B4)
where we have used $`[𝐱,𝐩]=i\mathrm{}`$ and $`𝐩=\frac{m}{i\mathrm{}}\left(𝐱\widehat{H}_A\widehat{H}_A𝐱\right)`$. $`\widehat{H}_A`$ is the free atomic Hamiltonian and gives the eigenvalue equation: $`\widehat{H}_A|n=E_n|n`$.

For the V-type three-level case, we take the three-level approximation in the above equation as ref.  for the two level atom case, obtaining

$`𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}`$ $``$ $`{\displaystyle \frac{2m}{\mathrm{}e^2}}[\mathrm{\Omega }_1(𝐞_{q\lambda }𝐝_1)(𝐝_1𝐞_{q^{}\lambda ^{}})`$ (B6)
$`+\mathrm{\Omega }_2(𝐞_{q\lambda }𝐝_2)(𝐝_2𝐞_{q^{}\lambda ^{}})],`$
where $`\mathrm{\Omega }_1=(E_{e_1}E_g)/\mathrm{}`$ and $`\mathrm{\Omega }_2=(E_{e_2}E_g)/\mathrm{}`$ are atomic transition frequencies for $`|g|e_1`$ and $`|g|e_2`$, respectively. $`𝐝_1=e_1|𝐝|g=g|𝐝|e_1`$ and $`𝐝_2=e_2|𝐝|g=g|𝐝|e_2`$ are the corresponding transition dipole moments. Choosing $`𝐞_{q\lambda }`$ as the following

$`𝐞_{q1}𝐝_1`$ $`=`$ $`d_1\text{}𝐞_{q1}𝐝_2=d_2,`$ (B7)
$`𝐞_{q^{}2}𝐝_1`$ $`=`$ $`0\text{}𝐞_{q^{}2}𝐝_2=0,`$ (B8)
then we get

$$𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}\frac{2m}{\mathrm{}e^2}[\mathrm{\Omega }_1d_1^2\delta _{\lambda ,1}\delta _{\lambda ,\lambda ^{}}+\mathrm{\Omega }_2d_2^2\delta _{\lambda ,1}\delta _{\lambda ,\lambda ^{}}].$$
(B9)
Substituting $`_{\lambda ,\lambda ^{}}𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}`$ into $`f(q,q^{})`$, we get Eq. (45) in section V. It is noticed that, for the case with single direction polarization of light, $`_{\lambda ,\lambda ^{}}𝐞_{q\lambda }𝐞_{q^{}\lambda ^{}}=1`$ strictly. The cut-off of the complete relation for the sum by only three levels, however, will lead to the departure from $`1`$. Only with this cut-off approximation the exactly-solvable mode is built for the two mode excitons coupling to the quantized electromagnetic fields.

## C Detailed non-perturbation calculation

In this Appendix, we will add the necessary details and list the more expatiatory expressions for section V. We start from the total interaction Hamiltonian

$`\widehat{H}_{int}`$ $`=`$ $`\mathrm{}{\displaystyle \underset{q,k}{}}G_1(q)O(k+q)\left[\widehat{A}_k+\widehat{A}_k^{}\right]\left[\widehat{a}_q+\widehat{a}_q^{}\right]`$ (C1)
$`+`$ $`\mathrm{}{\displaystyle \underset{q,k}{}}G_2(q)O(k+q)\left[\widehat{B}_k+\widehat{B}_k^{}\right]\left[\widehat{a}_q+\widehat{a}_q^{}\right]`$ (C2)
$`+`$ $`\mathrm{}{\displaystyle \underset{q,q^{},k}{}}{\displaystyle \frac{1}{\mathrm{\Omega }_1}}G_1(q)G_1(q^{})O(q^{}k)O(k+q)`$ (C4)
$`\times \left[\widehat{a}_q+\widehat{a}_q^{}\right]\left[\widehat{a}_q^{}+\widehat{a}_q^{}^{}\right]`$
$`+`$ $`\mathrm{}{\displaystyle \underset{q,q^{},k}{}}{\displaystyle \frac{1}{\mathrm{\Omega }_2}}G_2(q)G_2(q^{})O(q^{}k)O(k+q)`$ (C6)
$`\times \left[\widehat{a}_q+\widehat{a}_q^{}\right]\left[\widehat{a}_q^{}+\widehat{a}_q^{}^{}\right].`$
which includes the non-RWA terms and the self-interaction of the light field.

The half-side Fourier transformation (HSFT) of the Heisenberg equations governed by the total Hamiltonian (C1) are

$`(\omega \mathrm{\Omega }_1)A_k(\omega )`$ $`=`$ $`{\displaystyle \underset{q}{}}G_1(q)O(qk)\left[a_q(\omega )+a_q^{}(\omega )\right]`$ (C8)
$`+iA_k(0),`$
$`(\omega \mathrm{\Omega }_2)B_k(\omega )`$ $`=`$ $`{\displaystyle \underset{q}{}}G_2(q)O(qk)\left[a_q(\omega )+a_q^{}(\omega )\right]`$ (C10)
$`+iB_k(0),`$
for excitons, and

$`(\omega |q|c)a_q(\omega )`$ (C11)
$`=`$ $`{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(\omega )A_k^{}(\omega )\right]`$ (C12)
$``$ $`i{\displaystyle \frac{1}{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(0)A_k^{}(0)\right]`$ (C13)
$`+`$ $`{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(\omega )B_k^{}(\omega )\right]`$ (C14)
$``$ $`i{\displaystyle \frac{1}{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(0)B_k^{}(0)\right]`$ (C15)
$`+`$ $`ia_q(0),`$ (C16)
$`(\omega +|q|c)a_q^{}(\omega )`$ (C17)
$`=`$ $`{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(\omega )A_k^{}(\omega )\right]`$ (C18)
$`+`$ $`i{\displaystyle \frac{1}{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(0)A_k^{}(0)\right]`$ (C19)
$``$ $`{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(\omega )B_k^{}(\omega )\right]`$ (C20)
$`+`$ $`i{\displaystyle \frac{1}{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(0)B_k^{}(0)\right]`$ (C21)
$`+`$ $`ia_q^{}(0),`$ (C22)
for photons. Here, we have eliminated $`_q^{}G(q^{})O(q^{}k)\left[a_q^{}(\omega )+a_q^{}^{}(\omega )\right]`$ in the derivation of Eqs. (C4) and (C5).

Combining the above two equations, we get

$`(\omega ^2q^2c^2)\left[a_q(\omega )+a_q^{}(\omega )\right]`$ (C23)
$`=`$ $`i(\omega +|q|c)a_q(0)+i(\omega |q|c)a_q^{}(0)`$ (C28)
$`+2|q|c{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(\omega )A_k^{}(\omega )\right]`$
$`i{\displaystyle \frac{2|q|c}{\mathrm{\Omega }_1}}G_1(q){\displaystyle \underset{k}{}}O(kq)\left[A_k(0)A_k^{}(0)\right]`$
$`+2|q|c{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(\omega )B_k^{}(\omega )\right]`$
$`i{\displaystyle \frac{2|q|c}{\mathrm{\Omega }_2}}G_2(q){\displaystyle \underset{k}{}}O(kq)\left[B_k(0)B_k^{}(0)\right].`$
Substituting Eq. (C6) into Eq. (C2) and (C3), after a straightforward calculation, we obtain two coupled equations for the Frenkel exciton operators

$`{\displaystyle \underset{k^{}}{}}\left[(\omega ^2\mathrm{\Omega }_1^2)\delta _{kk^{}}{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{kk^{}}^{(1)}(\omega )\right]`$ (C30)
$`\times \left[A_k^{}(\omega )A_k^{}^{}(\omega )\right]`$
$`=`$ $`i\left[(\omega +\mathrm{\Omega }_1)A_k(0)(\omega \mathrm{\Omega }_1)A_k^{}(0)\right]`$ (C31)
$`+`$ $`2i\omega {\displaystyle \underset{q}{}}G_1(q)O(qk)\left[{\displaystyle \frac{a_q(0)}{\omega |q|c}}+{\displaystyle \frac{a_q^{}(0)}{\omega +|q|c}}\right]`$ (C32)
$``$ $`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(1)}(\omega )\left[A_k^{}(0)A_k^{}^{}(0)\right]`$ (C33)
$``$ $`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(3)}(\omega )\left[B_k^{}(0)B_k^{}^{}(0)\right]`$ (C34)
$`+`$ $`{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(3)}(\omega )\left[B_k^{}(\omega )B_k^{}^{}(\omega )\right],`$ (C35)
and

$`{\displaystyle \underset{k^{}}{}}\left[(\omega ^2\mathrm{\Omega }_2^2)\delta _{k,k^{}}{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{kk^{}}^{(2)}(\omega )\right]`$ (C37)
$`\times \left[B_k^{}(\omega )B_k^{}^{}(\omega )\right]`$
$`=`$ $`i\left[(\omega +\mathrm{\Omega }_2)B_k(0)(\omega \mathrm{\Omega }_2)B_k^{}(0)\right]`$ (C38)
$`+`$ $`2i\omega {\displaystyle \underset{q}{}}G_2(q)O(qk)\left[{\displaystyle \frac{a_q(0)}{\omega |q|c}}+{\displaystyle \frac{a_q^{}(0)}{\omega +|q|c}}\right]`$ (C39)
$``$ $`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(2)}(\omega )\left[B_k^{}(0)B_k^{}^{}(0)\right]`$ (C40)
$``$ $`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(3)}(\omega )\left[A_k^{}(0)A_k^{}^{}(0)\right]`$ (C41)
$`+`$ $`{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}{\displaystyle \underset{k^{}}{}}F_{kk^{}}^{(3)}(\omega )\left[A_k^{}(\omega )A_k^{}^{}(\omega )\right],`$ (C42)
where only the initial photon operators and the exciton operators are concerned. Three factors introduced in Eqs. (c7) and (c8) are

$`F_{kk^{}}^{(1)}(\omega )`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle \frac{2|q|cG_1^2(q)O(qk)O(k^{}q)}{\omega ^2q^2c^2}},`$ (C43)
$`F_{kk^{}}^{(2)}(\omega )`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle \frac{2|q|cG_2^2(q)O(qk)O(k^{}q)}{\omega ^2q^2c^2}},`$ (C44)
$`F_{kk^{}}^{(3)}(\omega )`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle \frac{2|q|cG_1(q)G_2(q)O(qk)O(k^{}q)}{\omega ^2q^2c^2}},`$ (C45)
which represent the overlap of exciton wave functions with different wave vectors. We take the photon normalization volume $`V`$ to be $`AL`$ where $`A`$ is the area of the crystal slab, and place the slab at the middle of the volume. When $`L`$ is sufficient large, the sum over $`q`$ can be replaced to a integral: $`_q\mathrm{}\frac{L}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑q\mathrm{}`$. Thus

$`F_{kk^{}}^{(i)}(\omega )`$ $`=`$ $`{\displaystyle \frac{Na\mathrm{\Omega }_if_i^2}{4\pi c^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q{\displaystyle \frac{O(qk)O(k^{}q)}{q^2(\frac{\omega }{c})^2}},`$ (C46)
By carrying out the integrations as in Ref. , the explicit expressions are

$`F_{kk^{}}^{(i)}(\omega )`$ $`=`$ $`{\displaystyle \frac{af_i^2\mathrm{\Omega }_i}{8Nc\omega }}[{\displaystyle \frac{\mathrm{sin}\frac{k^{}+\frac{\omega }{c}}{2}Na}{\mathrm{sin}\frac{k^{}+\frac{\omega }{c}}{2}a}}{\displaystyle \frac{e^{i\frac{k+\frac{\omega }{c}}{2}Na}}{\mathrm{sin}\frac{k+\frac{\omega }{c}}{2}a}}`$ (C47)
$``$ $`{\displaystyle \frac{\mathrm{sin}\frac{k^{}\frac{\omega }{c}}{2}Na}{\mathrm{sin}\frac{k^{}\frac{\omega }{c}}{2}a}}{\displaystyle \frac{e^{i\frac{k\frac{\omega }{c}}{2}Na}}{\mathrm{sin}\frac{k\frac{\omega }{c}}{2}a}}`$ (C48)
$`+`$ $`{\displaystyle \frac{\mathrm{sin}\frac{kk^{}}{2}Na}{\mathrm{sin}\frac{kk^{}}{2}a}}{\displaystyle \frac{\mathrm{sin}\frac{\omega }{c}a}{\mathrm{sin}(\frac{k+\frac{\omega }{c}}{2}a)\mathrm{sin}(\frac{k\frac{\omega }{c}}{2}a)}}].`$ (C49)
When $`N\mathrm{}`$, the first two terms tend to zero, thus

$$\underset{N\mathrm{}}{lim}F_{kk^{}}^{(i)}(\omega )=\frac{af_i^2\mathrm{\Omega }_i}{8c\omega }\frac{\mathrm{sin}\frac{\omega }{c}a}{\mathrm{sin}(\frac{k+\frac{\omega }{c}}{2}a)\mathrm{sin}(\frac{k\frac{\omega }{c}}{2}a)}\delta _{k,k^{}}.$$
(C50)
The above results will be used to determine $`a_q(\omega )a_q^{}(\omega )`$ explicitly in Sec. V.

For the single lattice layer case, ignoring the quantum noise terms proportional to $`a_q(0)`$ or $`a_q^{}(0)`$ induced by the background light field, we obtain the coupled equations of the exciton operators

$`\left[\omega ^2\mathrm{\Omega }_2^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )\right]\left[B_0(\omega )B_0^{}(\omega )\right]`$ (C51)
$`=`$ $`i\left[(\omega +\mathrm{\Omega }_2)B_0(0)(\omega \mathrm{\Omega }_2)B_0^{}(0)\right]`$ (C55)
$`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}F_{00}^{(2)}(\omega )\left[B_0(0)B_0^{}(0)\right]`$
$`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}F_{00}^{(3)}(\omega )\left[A_0(0)A_0^{}(0)\right]`$
$`+{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(3)}(\omega )\left[A_0(\omega )A_0^{}(\omega )\right],`$
and

$`\left[\omega ^2\mathrm{\Omega }_1^2{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )\right]\left[A_0(\omega )A_0^{}(\omega )\right]`$ (C56)
$`=`$ $`i\left[(\omega +\mathrm{\Omega }_1)A_0(0)(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$ (C60)
$`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_1}}F_{00}^{(1)}(\omega )\left[A_0(0)A_0^{}(0)\right]`$
$`2i{\displaystyle \frac{\omega }{\mathrm{\Omega }_2}}F_{00}^{(3)}(\omega )\left[B_0(0)B_0^{}(0)\right]`$
$`+{\displaystyle \frac{2\omega ^2}{\mathrm{\Omega }_2}}F_{00}^{(3)}(\omega )\left[B_0(\omega )B_0^{}(\omega )\right].`$
These two equations lead to two decoupled equations Eq. (90) and Eq. (95) in Sec. V for $`A_0(\omega )A_0^{}(\omega )`$ and $`B_0(\omega )B_0^{}(\omega )`$. The solutions of Eq. (90) and Eq. (95) are

$`A_0(\omega )A_0^{}(\omega )`$ (C61)
$`=`$ $`i{\displaystyle \frac{\omega ^2\mathrm{\Omega }_2^2+i\eta _2\omega }{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_1)A_0(0)(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$ (C64)
$`+{\displaystyle \frac{\mathrm{\Omega }_3\eta _3}{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_2)B_0(0)(\omega \mathrm{\Omega }_2)B_0^{}(0)\right]`$
$`{\displaystyle \frac{\eta _1}{\zeta (\omega )}}\left(\omega ^2\mathrm{\Omega }_2^2\right)\left[A_0(0)A_0^{}(0)\right],`$
and

$`B_0(\omega )B_0^{}(\omega )`$ (C65)
$`=`$ $`i{\displaystyle \frac{\omega ^2\mathrm{\Omega }_1^2+i\eta _1\omega }{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_2)B_0(0)(\omega \mathrm{\Omega }_2)B_0^{}(0)\right]`$ (C68)
$`+{\displaystyle \frac{\mathrm{\Omega }_3\eta _3}{\zeta (\omega )}}\left[(\omega +\mathrm{\Omega }_1)A_0(0)(\omega \mathrm{\Omega }_1)A_0^{}(0)\right]`$
$`{\displaystyle \frac{\eta _2}{\zeta (\omega )}}\left(\omega ^2\mathrm{\Omega }_1^2\right)\left[B_0(0)B_0^{}(0)\right],`$
where $`\zeta (\omega )`$ is defined in Eq. (LABEL:eq:zeta). The above two equations determine $`a_q(\omega )a_q^{}(\omega )`$ and give the non-zero electric field $`E(z,t)`$ (Eq. (65) in Sec. V)

$`E(z,t)=^{(+)}(z,t)+\mathrm{H}.\mathrm{c}.,`$
for $`zct<0`$, where

$`^{(+)}(z,t)`$ (C69)
$`=`$ $`\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_1\eta _1}{cA}}}\left[\left(\omega _1+\mathrm{\Omega }_1\right)A_0(0)+\left(\omega _1\mathrm{\Omega }_1\right)A_0^{}(0)\right]`$ (C77)
$`\times {\displaystyle \frac{(\omega _1^2\mathrm{\Omega }_2^2)e^{i\omega _1(t\frac{z}{c})}}{\left(\omega _1\omega _2\right)\left(\omega _1\omega _3\right)\left(\omega _1\omega _4\right)}}`$
$`+\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_1\eta _1}{cA}}}\left[\left(\omega _3+\mathrm{\Omega }_1\right)A_0(0)+\left(\omega _3\mathrm{\Omega }_1\right)A_0^{}(0)\right]`$
$`\times {\displaystyle \frac{(\omega _3^2\mathrm{\Omega }_2^2)e^{i\omega _3(t\frac{z}{c})}}{\left(\omega _3\omega _1\right)\left(\omega _3\omega _2\right)\left(\omega _3\omega _4\right)}}`$
$`+\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_2\eta _2}{cA}}}\left[\left(\omega _1+\mathrm{\Omega }_2\right)B_0(0)+\left(\omega _1\mathrm{\Omega }_2\right)B_0^{}(0)\right]`$
$`\times {\displaystyle \frac{(\omega _1^2\mathrm{\Omega }_1^2)e^{i\omega _1(t\frac{z}{c})}}{\left(\omega _1\omega _2\right)\left(\omega _1\omega _3\right)\left(\omega _1\omega _4\right)}}`$
$`+\sqrt{{\displaystyle \frac{\pi \mathrm{}\mathrm{\Omega }_2\eta _2}{cA}}}\left[\left(\omega _3+\mathrm{\Omega }_2\right)B_0(0)+\left(\omega _1\mathrm{\Omega }_2\right)B_0^{}(0)\right]`$
$`\times {\displaystyle \frac{(\omega _3^2\mathrm{\Omega }_1^2)e^{i\omega _3(t\frac{z}{c})}}{\left(\omega _3\omega _1\right)\left(\omega _3\omega _2\right)\left(\omega _3\omega _4\right)}},`$
here $`\omega _j`$, for $`j=1,2,3,4`$, are four roots of the characteristic equation of Eq. (97). Eq. (C77) shows that, in positive $`z`$ region, the electric field generated by the two-mode exciton is damped exponentially with two eigendecay rates $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$. Unlike the results of Sec. IV, the corresponding eigenmodes, however, are two linear combinations of the $`A_0`$ and $`B_0`$ modes. This is because we include non-RWA terms and MPP in the derivation of Eq. (90) and Eq. (95) in Sec. V.