Trapping Effects in Quantum Atomic Arrays

Quantum emitters, particularly atomic arrays with subwavelength lattice constant, have been proposed to be an ideal platform for studying the interplay between photons and electric dipoles. In this work, motivated by the recent experiment [1], we develop a microscopic quantum treatment using annihilation and creation operator of atoms in deep optical lattices. Using a diagrammatic approach on the Keldysh contour, we derive the cooperative scattering of the light and obtain the general formula for the $S$ matrix. We apply our method to study the trapping effect, which is beyond previous treatment with spin operators. If the optical lattices are formed by light fields with magical wavelength, the result matches previous results using spin operators. When there is a mismatch between the trapping potentials for atoms in the ground state and the excited state, atomic mirrors become imperfect, with multiple resonances in the optical response. We further study the effect of recoil for large but finite trapping frequency. Our results are consistent with existing experiments.

E 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " q u J I P S D L 6 w A 4 x 5 W 3 7 z o H g 2 o j V 7 A = " > A A A C z 3 i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w V d I q 6 E o K I r h s w b Z C W 8 p k O m 1 D 8 y K Z K K V U 3 P o D b v W v x D / Q v / D O m I J a R C c k O X P u P W f m 3 m u H r h N L y 3 r N G A u L S 8 s r 2 d X c 2 v r G 5 l Z + e 6 c R B 0 n E R Z 0 H b h B d 2 y w W r u O L u n S k K 6 7 D S D D P d k X T H p 2 r e P N G R L E T + F d y H I q O x w a + 0 3 c 4 k 0 S 1 2 x 6 T Q 7 s / u e h a 0 2 6 + Y B U t v c x 5 U E p B A e m q B v k X t N F D A I 4 E H g R 8 S M I u G G J 6 W i j B Q k h c B x P i I k K O j g t M k S N t Q l m C M h i x I / o O a N d K W Z / 2 y j P W a k 6 n u P R G p D R x Q J q A 8 i L C 6 j R T x x P t r N j f v C f a U 9 1 t T H 8 7 9 f K I l R g S + 5 d u l v l f n a p F o o 9 T X Y N D N Y W a U d X x 1 C X R X V E 3 N 7 9 U J c k h J E 7 h H s U j w l w r Z 3 0 2 t S b W t a v e M h 1 / 0 5 m K V X u e 5 i Z 4 V 7 e k A Z d + j n M e N M r F 0 l G x X D s u V M 7 S U W e x h 3 0 c 0 j x P U M E l q q i T d 4 h H P O H Z q B m 3 x p 1 x / 5 l q Z F L N L r 4 t 4 + E D 9 l O U D Q = = < / l a t e x i t > p T E 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " M H B p V / F X y u X Z T c x F E y C + W y x C 6 c g = " > A A A C 2 X i c j V H L S s N A F D 2 N r 1 p f 9 b F z E y y C q 5 J W Q
In most of these works, atoms are treated as point-like with no motional degree of freedom.The evolution of the system is described by using non-Hermitian Hamiltonian or Lindblad master equation [16,17], with spin degree of freedom σ − im = |r m , g r m , e i |.Here |r m , g is the s-wave ground state for the atom at position r m .|r m , e i is the p-wave excited labeled by the dipole moment d = d e i of the corresponding transition g → e i .However, in real experiments, the system consists of atoms moving in optical lattices.For deep optical lattices, although atoms are trapped near the potential minimum, the wave function for the motional degree of freedom may still play a role.Moreover, the consequence of fractional filling has been studied in the experiment.It is difficult to analyze the absence of an atom in the spin-operator language, and consequently, theoretical predictions for the fractional filling case are still absent.
In this work, we overcome this difficulty by using a microscopic model for the coupled system consisting of atoms in deep optical lattices and photons.After making plausible assumptions, we derive the cooperative response of the system using a diagrammatic approach on the Keldysh contour.By summing up bubble diagrams with dressed Green's function, we obtain neat results for the transmission coefficient and the reflection coefficient, with the contribution from the motional wave function.Our result matches the previous analysis for unit filling when the potential of the excited state atoms are the same as that of the ground state.We study the effect of the discrepancy of optical lattices for the ground state and excited state atoms, where the transition of the internal state can be accompanied by transitions in the motional degree of freedom.In particular, we find that multiple resonances can exist in the response function.The cooperative linewidth is linear in n, consistent with the experimental observation and previous theoretical predictions [38].We further study the effect of recoil for large but finite trapping frequency.
2 Diagrammatic Approach to Quantum Atomic Arrays

Model
We consider coupled systems with atoms and photons.Th Hamiltonian reads Here the first term is the Hamiltonian of the electromagnetic field The second term describes the motion of atoms in optical lattices V g/e i (r) describes the optical lattice potential for ground/excited-state atoms.We have set = 1 and m = 1.We assume each site is occupied by at most one atom, which corresponds to choosing fermionic commutation relation {ψ † a (r), ψ b (r )} = δ ab δ(r − r ).The last term describes the interaction between atoms here e i is the unit polarization vector along the i direction and The full Hamiltonian (1) describes general model with interaction between atoms and light to the order of electric dipole transition.For atomic arrays, the ground state particle is always tightly trapped near the local minimum of optical lattices, with a spread of wave function σ a 0 , where a 0 is the lattice constant [1].Assuming the excited state is also deeply trapped, we expand Here ϕ a (r)/ϕ i,a (r) is the motional wave function for ground/excited-state atoms near the local minimum r n = 0 with the energy ε a /ε i,a .We neglect the tunneling between different sites, which is suppressed expoenentially and as a result the local motional wave function conincident with the Wannier function.We have r n = a 0 (n 1 e 1 + n 2 e 2 ), where we use a single index n to represent (n 1 , n 2 ) for conciseness.The commutation relation for ψ a η (r n ) now becomes {ψ a, † η (r m ), ψ b ξ (r n )} = δ ηξ δ ab δ mn .Using (6), the Hamiltonian H A and H int can be simplified.We have and ( 4) becomes with Equation ( 2), ( 7) and ( 8) describe the dynamics of the atomic array.Initially, we prepare all atoms in the s-wave internal ground state |g with motional ground state ϕ 0 (r).The number of atoms in the excited states are suppressed due to the violation of energy conservation.We further add an external probe light, at fixed frequency ω, which is near-resonant with δ ≡ ω − ω 0 ω, ω 01 .The electric field reads E 0 (r) = E 0 e ik•r with c|k| = ω.We take c = 1 from now on for conciseness.This probe corresponds to the incident light in the scattering experiment.Its coupling to atoms reads We assume the field strength E 0 is weak and the response can be analyzed using the linear response theory.The total electric field including the incident light and the scattered light then reads Far from the atomic array, when only a single diffraction order exists, we expect and S(ω, k ) is the corresponding S matrix.

Diagrammatic Expansion
The contribution to the scattered light E(ω, r) can be efficiently organized using the path-integral formulism.In particular, we work on the Keldysh contour [42], which contains a forwardly evolving branch and a backwardly evolving branch, corresponding to e −iHt and e iHt in the Heisenberg evolution.It is one of standard techniques for analyzing quantum many-body dynamics and systems with disorders.
The expectation of fluctuation field becomes non-zero due to the coupling to atoms.Diagrammatically, we have Here we use the wavy line to represent the propagation of photons.G E R is the retarded Green's function matrix of E in free space defined as In frequency and momentum space, we have Here G(ω, k) is the standard dyadic Green's function [16,43].Note that we have added an additional tilde for the Green's function of photons in momentum space to avoid possible confusion.The local dipole moment p − is related to the incident light E 0 by the Kubo formula [44] p We use the double solid line for the retarded Green's function for dipole momentums G p R (ω, r) and r nm ≡ r n − r m .This is consistent with the semi-classical analysis [16].The remaining task is to derive approximate formula for G p R (ω, r), which includes renormalization due to the coupling with photons.In this work, we take diagrams with single excitation which conserves the total energy.We first consider the correction of the excited state Green's function G e i R (ω, r, r ) by emission and absorption of photons.As we will see, since the wave function for ground-state atoms is localized, only G e i R (ω, r, r ) with r ≈ r ≈ r n contributes to the light scattering.The bare Green's function near r n = 0 reads G 0,e i R (q 0 , r, r ) ≈ The Schwinger-Dyson equation reads (G Here n a≥1 = 0 and n 0 = n is equal to the filling fraction.The appearance of G ii (ω, 0) = e i • G(ω, 0) • e i owes to the approximation in (8) by using E(r n ) instead of E(r).The real-part of G(ω, 0) contributes to the lamb shift, which can be absorbed in the definition of ω 0 .As a result, we only keep the imaginary part G(ω, 0) = iω/6π × 1.We also assume δ, ε a , ε i,a ω, and the resonance frequency ω is much larger than the loop frequency, which is an analogy of the Markovian approximation in the master equation [16].Note that in the path-integral approach, the Green's function is defined by adding an addtional particle on top of the many-body system with filling n, as a result, the Pauli exclusion principle exists and contributes to the (1 − n a ) factor.The natural linewidth of a transition with frequency ω is known to be γ = ω 3 0 d 2 /3π 0 .This leads to where we have used the completeness of local wave functions a ϕ a (r)ϕ * a (r ) = δ(r − r ).Having obtained the dressed Green's function, we turn to the calculation of G p R (ω, r).Motivated by the standard Random Phase Approximation (RPA) in interacting fermions [45], we consider the diagrams Note that in our diagrammatic approach, r n can be equal to r m , which is important as we will see later.The thick solid line represents the normalized Green's function G e i R .The first bubble diagram, which is an elementary building block, is given by Here we write Π R (ω) as a diagonal matrix for later convenience.
A is the Keldysh Green's function, and F η = (1−2n η ) is the quantum distribution function [42].Here we use • to represent the inner product of functions by integration.This leads to It can be further simplified by noticing that As a result, we have Then we can sum over the diagrams with multiple bubbles in (20).This gives Since the summation in ( 16) is descrete, the Fourier transform here is defined as In particular, the denominator of ( 25) is a generalization of the corresponding result under non-Hermitian Hamiltonians, which is ω1 − H eff .As we will see later, (21) takes such a form only for unit filling and V e i (r) = V g (r).This implies the breakdown of non-Hermitian Hamiltonian description for general setups.Then, using the relation (15), we obtain the relation between p − and E 0 in momentum space as Finally, for a single diffraction order, using (13), α is related to the S matrix as [16] S Here ) is in z direction, and at the same time z < 0. In other cases, ξ i j = 1.
Further simplification is possible for the normal incident case as in the experiment [1], where we have k x = k y = 0, k z = ω.Since, E 0 lies in the x-y plane, we have P = 1.Moreover, due to the rotational symmetry, GE R (ω, 0) is also a diagonal matrix.Following the convention [16], we define ∆(k which also become scalars ∆ and Γ in the x-y plane for normal incident light.In particular, it is known that Γ + γ = γ 3π a 2 ω 2 [16].As a result, the S matrix is diagonal and there is no mixing between contributions from different excited states e i .For the incident light polarized in i 0 direction, the only relevant response function is Π R (ω) ≡ Π R,i 0 i 0 (ω), which is related to the local response π(ω) ≡ π i 0 (ω).We also drop the i 0 index in ε a ≡ ε i 0 ,a and ϕ a ≡ ϕ i 0 ,a for conciseness.From now on, we focus on this normal incident case unless mentioned otherwise.
Using these definitions, we have (29) Since π(ω) is independent of n, the cooperative linewidth is linear in filling fraction n, consistent with the previous work [38].
Equation ( 24) and ( 29) together determine the cooperative optical response of the atomic array.In the next sections, we first validate our path-integral approach by showing that our result is consistent with previous literatures when the optical lattice is formed by a light with the magic wavelength.In this case, the trapping poential of the excited state is the same as that of the ground state.We then consider the effect of the trapping mismatch, which have been observed in the recent experimental realization of the atomic array [1].

Trapping Effects in Quantum Atomic Arrays
In this section, we analyze ( 24) and ( 29) in several limits.We first consider the case with V e i 0 (r) = V g (r), and show that the result then matches the spin operator result.We also add comments on the relation between our approach and the Schwinger boson representation of spins [38,39] 2 .We then study the effect of trapping mismatch, which leads imperfectness of mirrors and multiple resonances.We finally go beyond (24) and ( 29) by analyzing the recoil effect.

Magic Wavelength
We begin with the special case V e i 0 (r) = V g (r), which is valid when the optical lattice is formed by a light with the magic wavelength [46].In this case, only the a = 0 term in (24) contributes and the transition of internal state does not couple to the motional degree of freedom.Moreover, there is no r Z o V W q l 1 9 5 l q 5 T L N N r 4 t 6 + E D Q 5 m P z A = = < / l a t e x i t > (a) dependence on the detailed shape of the potential.This leads to

< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 c C 8 T K 4 h p + h J T z I e F r h w z w C y G f Y = " >
Here we have assumed the normal incidence for the probe light.The cooperative linewidth becomes γ + nΓ, consistent with the experimental observation and numerical simulation in [1], and the previous work [38].For n < 1, we find |S | < 1 even at the resonance and the mirror becomes imperfect.The transmission coefficient T = |1 + S | 2 and reflection coefficient R = |S | 2 are found to be The filling-normalized absorptance A = (1 − T )/n and reflectance R = R/n can the be computed straightforwardly.
We plot the numerical result for ωa 0 = 2π × 0.68 as in the experiment [1] for various n in Figure 2, where we have ∆/γ ≈ 0.18 and Γ/γ ≈ −0.48.All above results reduces to the semi-classical results using spin operators when n = 1, where the frequency shift is ∆ and the linewidth becomes γ + Γ.On the other hand, for n → 0, we get back to the single-atom response with natural linewidth γ.As observed in the experiment [1], generally, we have T + R < 1.This is due to the fact that the self-energy of the excited state (18) contains the contribution of spontaneous emission of photons in arbitrary directions with random phases, which can not be observed by averaged E tot .However, the corresponding contribution exists if we measure energy density of electromagnetic field E 2 (r) [34].The filling-normalized absorptance A show a weak dependence of n, while R vanishes as n → 0, consistent with the experimental observation and numerical simulation in [1].
Finally we comment on the relation between our results and the Swchinger boson/fermion representation of spins [38,39].Using (30), we find the bubble reads Π R (ω) = d 2 n/(δ + i γ(1−n) 2 ).As we mentioned in the last section, the factor of (1−n) exist due to the Pauli exclusion principle.This seems to be unphysical since after the ground state particle being excited on some site, no Pauli exclusion factor is needed.However, the contribution from −iγn/2 indeed cancels out with the corresponding contribution of the Green's function of photons GE R (ω, k ) in (29) due to (Recall that the definition of (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " a T N H   Γ does not contain r n = 0, which is equal to γ.) This cancellation can be dated back to the cancellation between diagrams.Let's consider diagrams with one internal photons.Before contractions between ψ g and ψ † g , it takes the form If we contract ψ g (1) with ψ † g (2) and ψ g (2) with ψ † g (1), this leads to the diagram in self-energies which contains the unwilling factor of (1 − n).To realize the fact that when the ground state particle is excited by ψ † g (2), there is already no occupation due to ψ g (2), we also need to take into account the contribution by contracting ψ g (1) with ψ † g (1) and ψ g (2) with ψ † g (2).More explicitly, we have However, the new diagram is just the bubble diagram, which has been taken into account in our diagrammatic expansion (20).As a result, by summing up self-energies and bubbles, we find the correct result.
Due to this cancellation, we can alternatively drop the factor of (1 − n) in Π R (ω), and restrict the summation to r 0 in (26).This is then consistent with rules in [38,39] using Schwinger particle (b)

< l a t e x i t s h a 1 _ b a s e 6 4 = " a T N H
r Z o V W q l 1 9 5 l q 5 T L N N r 4 t 6 + E D Q 5 m P z A = = < / l a t e x i t > (a) representation, where atom on the same site can not appear twice.On the other hand, if the atoms are not trapped in optical lattices, the density-density correlation indeed plays an role [47].In this case, the corresponding contribution in the self-energy should be kept.We give an example in Appendix A.

Trapping Mismatch
In this section, we discuss the effect of having V e i 0 (r) V g (r).To make this problem analytically solvable, we expand the potential of near the minimum of each site and use the approximation of 3D isotropic harmonic potential.Ground-state atoms |g and excited-state atoms |e i 0 have a trapping frequency ω h g and ω h e correspondingly.The motional ground state wave function ϕ 0 (r) reads Under this approximation, the analytical expression for π(ω) can be obtained by relating π(ω) to the single-particle propagator in harmonic traps

(b)
< l a t e x i t s h a 1 _ b a s e 6 4 = " a T N H  The parameters in the experiment [1] corresponds to ω h g < γ, but at the same order ∼ MHz.We plot our results (29) for different ω h e /γ, ω h g /γ and n in Figure 3 and 4. We first fix ω h g /γ = 1/4 and study the effect of ω h e ω h g for small ω h e < γ.As shown in Figure 3 (a) and (b), both reflection coefficient R(ω) and transmission coefficient T (ω) show a single peak near δ − n∆ − 3(ω h e − ω h g )/2 = 0.For either ω h e > ω h g or ω h g > ω h e , the atomic mirror becomes imperfect with max R < 1 and min T > 0. Motivated by the experimental result, we study the the cooperative linewidth of the atomic array by fitting the numerical result for R(ω) near δ = n∆ + 3(ω h e − ω h g )/2 as and define T res = T (n∆ + 3(ω h e − ω h g )/2 + δ 0 ).Rres and A res can then be computed correspondingly using R res and T res .The numerical results in (c-d) show Rres and A res also decreases when ω h e ω h g .The cooperative linewidth Γ cor is linear in n, with similar slope for different ω h e < γ.We then consider larger ω h e ≥ γ in Figure 4. Now as shown in Figure 4 (a) and (b), multiply peaks appear in both reflection coefficient R(ω) and transmission coefficient T (ω).The center of peaks locates near energy 2nω h e , where the transition from |g to |e i 0 is accompanied with the excitation of motional degree of freedom.As an example, we fit the position of the second peak ∆ω, and plot it as a function of ω h e in (d).Similar to the small ω h e case, the cooperative linewidth Γ cor is still linear in n.However, their slope show dependence on ω h e .We finally study A res and Γ cor as a function of ω h g .We fix a small ω h e = γ/20, as an analogy of the anti-trapped excited state in experiment [1], and tune ω h g .As shown in Figure 5, we find when ω h g becomes larger, the absorption rate decreases and the cooperative linewidth becomes larger.This is due to the increase of the trapping mismatch for small ω h e .For small ω h g , the decrease in A res and the increase of the decay rate show quadratic dependence, while for large ω h g , the dependence becomes linear.This is a close analogy of the experimental observation in [1].

Recoil Effects
Now we go beyond the limit of σ a 0 and consider corrections to the leading order of η = σ/a 0 due to the recoil of atoms.The recoil effect has been discussed in previous works [40,41] using the Lindblad master equation.
We adopt the isotropic harmonic trap approximation used in the last subsection.In our approach, the recoil effects can be analyzed by using (4) without the approximation (8) as explained in the Appendix C.However, the same result can also be obtained using direct physical intuition.Our result (24) for π(ω) has a simple physical meaning: it measures the optical response of a single particle in the harmonic trap, with a lift time γ for the excited state particle.Without the recoil effect, we take the inner product between wavefunctions ϕ 0 (r) and ϕ a (r) to determine the transition rate for the motional degree of freedom.To take the recoil effect into account, we consider absorbing a photon with momentum k 1 and emitting a photon with momentum k 2 .Physically, we expect the local response π(ω) takes the form Since the photons can propagate in any direction within x − y plane, we further average over the direction of the momentum k i = k(cos θ i , sin θ i , 0): We focus on the case with ω h g = ω h e .For small η = σ/a 0 = k/ ω h g , we can expand π(ω) to obtain To the leading order O(η 2 ), we find only a renormalization of the residue near δ + iγ/2 = 0.At the sub-leading order, O(η 4 ) we see a non-trivial contribution where atoms are excited due to the recoil even without trapping mismatch ω h e = ω h g , because the recoil increases the energy of atoms.However, when we further plot the reflection and transmission coefficients (see Figure 6), we find no visible peaks at δ − n∆ = 2ω h e even for ω h g > γ due to the suppression of η 4 .We find R(ω) and T (ω) only show weak dependence on η for small η.In particular, the parameter regime in the experiment [1] corresponds to η ∼ 0.1, which is almost the same as η = 0.This justifies our discussions in previous sections.

Summary and Outlook
In this work, we study quantum atomic arrays using a microscopic model with atoms in optical lattices.We take a diagrammatic approach with PRA-like diagrams and obtain concise results for transmission and reflection coefficients.We find trapping mismatch can result in the imperfectness of mirrors.Multiple peaks exist when the local trapping frequency of the excited state ω g e ∼ γ.We also study the trapping frequency effects on the cooperative lifetime, and the effect of recoil for large but finite trapping frequency.
Our results can be tested in the experimental platforms similar to that in [1].Recently, there are also experimental studies on the Pauli blocking of light scattering in degenerate fermions [48,49].The diagrammatic approach developed here can also be applied to study the optical response of degenerate fermion gases.

A Homogeneous Mirrors
Now we consider the case where the atom gas is homogeneous in x−y plane at z = 0, similarly to [49].The model of the full system reads The first step is again to determine the renormalization of the excited state Green's function.The self-energy reads This leads to Here ε q = q 2 /2.We then compute the Π R .The result reads Summing up the contribution from photons, we find Here since the system is homogeneous, the Fourier transform is The relation between E tot and α becomes The S matrix reads Again we consider the normal incident light.We further assume the density of the system is low as in [47].We have Π R (ω, 0) i j = dq (2π) 2 d 2 n F (ε q ) (δ + iγ/2)1 − δH(q) = n 2D d 2 δ + iγ/2 δ i j + ω 2 d 4 0 (δ + iγ/2) 2 dqdq (2π) 4 e i • G(ω, q − q , z = 0) • e j n F (ε q )n F (ε q ) We find α(ω, 0) = 1 Here we have As a result, .This takes the similar form as results in [47] for the 3D case.This is the contribution from the density-density correlation in free fermion gases.Finally, we have which means δΠ R (ω) effectively shifts the resonant energy and the decay rate.
B The Analytical Formula for π(ω) In this Appendix, we present detailed derivation of the analytical formula for π(ω).We trick is to use the transformation to the time domain (51) Here we have assumed the integral over τ is convergent by restricting the δ + ε 0 < 3ω h e /2.After the integration, analytical continuation can be applied to release this restriction.Here K ω h e (τ, r, r ) is the imaginary time Green's function in a harmonic trap with trapping frequency ω h e .We have q −p−1 (−2p(q − 1) − q + 2)B q p + 1, 1 2 − (2p + 3)B q p + 1, 3

C The Derivation of the Recoil Effects
Here we give the derivation of (37) using the diagrammatic approach.For simplicity, we directly take diagrams under the rule of Schwinger bosons/fermions as discussed section 3.1.As a result, the self-energy of excited state Σ e i R (q 0 ) is just a constant −iγ/2, and the summation in ( 20) is restricted by r m r n .
To determine π(ω), we now examine a single diagram We have separated the integral over the full space into a summation over lattice sites r n , and an integral near each sites r i .For small σ a 0 , the dominate contribution comes from r mp r 2 − r 1 and r nm r 3 − r 2 .Moreover, for normal (or nearly normal) incident light, we are probing the system with small k , which comes from contributions at large r nm and r mp .We then expand GE R and take the standard approximation at long distance [50]:

Figure 1 :
Figure 1: Schematics of the model considered in this work: the atomic array in the optical lattices at fractional filling.

Figure 2 :
Figure 2: Numerical result for the fractional-filling effect with normal incident light with ωa 0 = 2π × 0.68.Here we take V e i 0 (r) = V g (r).(a).The reflection coefficient R(ω) as a function of detuning δ − n∆ for different filling fraction n. (b).The transmission coefficient T (ω) as a function of detuning δ − n∆ for different filling number n. (c).The filling-normalized absorptance A and reflectance R, together with T + R at cooperative resonance δ = n∆, as a function of filling fraction n.
n s m R s 3 o l I B e S U o b e 6 S J K E 8 S 1 q f Z J p 4 a Z 8 3 + 5 j 0 x n v p u Y / r 7 m d e I W I U b Y v / S T T P / q 9 O 1K P R x Y m o Q V F N s G F 0 d y 1 x S 0 x V 9 c / t L V Y o c Y u I 0 7 l F c E m Z G O e 2 z b T S J q V 3 3 1 j P x N 5 O p W b 1 n W W 6 K d 3 1 L G r D 7 c 5 y z o H l Q c Q 8 r B + d H p e p p N u o 8 d r C L M s 3 z G F X U U E e D v A d 4 x B O er Z o V W q l 1 9 5 l q 5 T L N N r 4 t 6 + E D Q 5 m P z A = = < / l a t e x i t > (a)< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 c C 8 T K 4 h p + h J T z I e F r h w z w C y G f Y = " >A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o C s p u O m y o n 1 A L T K Z T m t o m o R k o p Q i + A N u9 d P E P 9 C / 8 M 4 4 B b W I T k h y 5 t x 7 z s y 9 1 4 s D P 5 W O 8 5 q z 5 u Y X F p f y y a 9 y g u C T O j n P b Z N p r E 1 K 5 7 6 5 n 4 m 8 n U r N 6 z L D f F u 7 4 l D d j 9 O c 5 Z 0 D y o u I e V g / O j U v U 0 G 3 U e O 9 h F m e Z 5 j C p q q K N B 3 g M 8 4 g n P V s 0 K r d S 6 + 0 y 1 c p l m G 9 + W 9 f A B R f q P z Q = = < / l a t e x i t > (d) < l a t e x i t s h a 1 _ b a s e 6 4 = "

Figure 3 :
Figure 3: Numerical result for the trapping effect with normal incident light with ωa 0 = 2π × 0.68.We fix ω h g = γ/4 and consider ω h e < γ. (a).The reflection coefficient R(ω) as a function of detuning for n = 1 with different ω h e .(b).The transmission coefficient T (ω) as a function of detuning for n = 1 with different ω h e .(c).The fitted A res and Rres as a function of n for different ω h e /ω h g .Here the dashed lines corresponds to A res .(d).The fitted linewidth Γ cor as a function of n for different ω h e .

Figure 4 :
Figure 4: Numerical result for the trapping effect with normal incident light with ωa 0 = 2π × 0.68.We fix ω h g = γ/4 and consider ω h e γ in (a-c).(a).The reflection coefficient R(ω) as a function of detuning for n = 1 with different ω h e .(b).The transmission coefficient T (ω) as a function of detuning for n = 1 with different ω h e .(c).The fitted linewidth Γ cor as a function of n for different ω h e .(d).The fitted center of the second peak ∆ω as a function of ω h e for different filling fraction n.Here the dashed line corresponds to ∆ω = 2ω h e .
. The results are presented in Appendix B. It contains multiple resonances near δ = (3/2 + 2n)ω h e − 3ω h g /2, broadened by the natural lifetime γ of the excited state.For ω h e γ, this leads to different peaks in the spectral −Im π(ω)/π.For ω h e γ, different resonances merges, and only a single peak exists.
e 2 z b T S J q V 3 31 j P x N 5 O p W b 1 n W W 6 K d 3 1 L G r D 7 c 5 y z o H l Q c Q 8 r B + d H p e p p N u o 8 d r C L M s 3 z G F X U U E e D v A d 4 x B O e r Z o V W q l 1 9 5 l q 5 T L N N r 4 t 6 + E D Q 5 m P z A = = < / l a t e x i t > (a)< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 c C 8 T K 4 h p + h J T z I e F r h w z w C y G f Y = " >A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o C s p u O m y o n 1 A L T K Z T m t o m o R k o p Q i + A N u9 d P E P 9 C / 8 M 4 4 B b W I T k h y 5 t x 7 z s y 9 1 4 s D P 5 W O 8 5 q z 5 u Y X F p f y y

Figure 5 :
Figure 5: Numerical results for the effect of the trapping mismatch with normal incident light and ωa 0 = 2π × 0.68.We fix small excited state trapping frequency ω h e = γ/20.(a).The fitted A res as a function of ω h g for different n. (b).The fitted linewidth Γ cor as a function of ω h g for different n.

) 2 .
B z (a, b) is the incomplete beta function defined as B z (a, b) = z 0 t a−1 (1 − t) b−1 dt.For ω h e = ω h g , one can check that above result can be simplified as π(ω) −1 = a 0 − 3ω h e /2.