Trapping Effects in Quantum Atomic Arrays

Pengfei Zhang

Institute for Quantum Information and Matter and Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA

*<EMAIL_ADDRESS>
August 27, 2024

§ ABSTRACT

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 <cit.>, 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 $\bm{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.

§ INTRODUCTION

The ability to coherently storing photons and controlling their interaction with quantum matters is of vital importance for quantum science. Although single atoms and photons usually interact less efficiently, ensembles of atoms can show a cooperative response of photons. As an example, superradiance can be realized when the radiations between atoms interfere constructively <cit.>. Recently, atomic arrays with subwavelength lattice structures are found to be an ideal platform where electric dipole-dipole interactions between atoms are mediated by photons <cit.>. The analysis shows the atomic arrays exhibit subradiance and are nearly perfect mirrors for a wide range of incident angles <cit.>, as observed in recent experiments <cit.>. Later, there are many theoretical studies on the fruitful physics in atomic arrays <cit.>. For example, there are proposals for realizing non-trivial topology in atomic arrays <cit.>, controlling atom-photon interaction using atomic arrays <cit.>, and efforts in understanding their subradiant behaviors and ability of photon storage <cit.>.

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

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 <cit.>, with spin degree of freedom $\sigma^{-}_{im}=|\mathbf{r}_m,g\rangle \langle \mathbf{r}_m,e_i|$. Here $|\mathbf{r}_m, g\rangle$ is the $s$-wave ground state for the atom at position $\mathbf{r}_m$. $|\mathbf{r}_m, e_i\rangle$ is the $p$-wave excited labeled by the dipole moment $\mathbf{d}=d~\mathbf{e}_i$ of the corresponding transition $g\rightarrow 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 <cit.>. We further study the effect of recoil for large but finite trapping frequency.

§ DIAGRAMMATIC APPROACH TO QUANTUM ATOMIC ARRAYS

§.§ Model

We consider coupled systems with atoms and photons. Th Hamiltonian reads

\begin{equation}\label{eq:Htot}
\end{equation}

Here the first term is the Hamiltonian of the electromagnetic field

\begin{equation}\label{eq:HEM}
H_{\text{EM}}=\int d\mathbf{r}~\left(\frac{\epsilon_0}{2}\mathbf{E}(\mathbf{r})^2+\frac{\mu_0}{2}\mathbf{H}(\mathbf{r})^2\right),
\end{equation}

The second term describes the motion of atoms in optical lattices

\begin{equation}
\begin{aligned}
H_{\text{A}}&=\int d\mathbf{r}~ \sum_i\psi^\dagger_{e_i}(\mathbf{r})\left( \omega_0 -\frac{\nabla^2}{2}+V_{e_i}(\mathbf{r})\right)\psi_{e_i}(\mathbf{r})+\int d\mathbf{r}~\psi^\dagger_{g}(\mathbf{r})\left( -\frac{\nabla^2}{2}+V_g(\mathbf{r})\right)\psi_{g}(\mathbf{r}),
\end{aligned}
\end{equation}

$V_{g/{e_i}}(\mathbf{r})$ describes the optical lattice potential for ground/excited-state atoms. We have set $\hbar=1$ and $m=1$. We assume each site is occupied by at most one atom, which corresponds to choosing fermionic commutation relation $\{\psi^\dagger_{a}(\mathbf{r}),\psi_{b}(\mathbf{r}')\}=\delta_{ab}\delta(\mathbf{r}-\mathbf{r}')$. The last term describes the interaction between atoms

\begin{equation}\label{eq:Hint}
H_{\text{int}}=-\int d\mathbf{r}~\sum_{i}\left(P^+_i(\mathbf{r})+P^-_i(\mathbf{r})\right)\mathbf{e}_i\cdot \mathbf{E}(\mathbf{r}),
\end{equation}

here $\mathbf{e}_i$ is the unit polarization vector along the $i$ direction and

\begin{equation}
\begin{aligned}
\end{aligned}
\end{equation}

The full Hamiltonian (<ref>) 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 $\sigma\ll a_0$, where $a_0$ is the lattice constant <cit.>. Assuming the excited state is also deeply trapped, we expand

\begin{equation}\label{eq:approx_psi}
\begin{aligned}
\psi_g(\mathbf{r})&\approx\sum_{n,a}\varphi_a(\mathbf{r}-\mathbf{r}_n)\psi_{g}^a(\mathbf{r}_n),\ \ \ \ \ \
\psi_{e_i}(\mathbf{r})&\approx\sum_{n,a}\varphi_{i,a}'(\mathbf{r}-\mathbf{r}_n)\psi_{e_i}^a(\mathbf{r}_n).
\end{aligned}
\end{equation}

Here $\varphi_a(\mathbf{r})/\varphi_{i,a}'(\mathbf{r})$ is the motional wave function for ground/excited-state atoms near the local minimum $\mathbf{r}_n=\mathbf{0}$ with the energy $\varepsilon_a/\varepsilon_{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 $\mathbf{r}_n=a_0(n_1\mathbf{e}_1+n_2\mathbf{e}_2)$, where we use a single index $n$ to represent $(n_1,n_2)$ for conciseness. The commutation relation for $\psi_\eta^a(\mathbf{r}_n)$ now becomes $\{\psi_{\eta}^{a,\dagger}(\mathbf{r}_m),\psi_{\xi }^b(\mathbf{r}_n)\}=\delta_{\eta\xi}\delta_{ab}\delta_{mn}$. Using (<ref>), the Hamiltonian $H_\text{A}$ and $H_\text{int}$ can be simplified. We have

\begin{equation}\label{eq:HAmodel}
\begin{aligned}
\end{aligned}
\end{equation}

and (<ref>) becomes

\begin{equation}\label{eq:Hintmodel}
H_{\text{int}}=-\sum_{n,i} \left(p_i^+(\mathbf{r_n})~\mathbf{e}_i\cdot \mathbf{E}(\mathbf{r_n})+\text{H.C.}\right),
\end{equation}

\begin{equation}\label{eq:p+}
\begin{aligned}
p^{+}_i(\mathbf{r_n})&=d\sum_{ab}\int d\mathbf{r}~\varphi_{i,a}'(\mathbf{r})^*\varphi_b(\mathbf{r})~\psi^{a,\dagger}_{e_i}(\mathbf{r}_n)\psi_{g}^b(\mathbf{r}_n).
\end{aligned}
\end{equation}

Equation (<ref>), (<ref>) and (<ref>) describe the dynamics of the atomic array. Initially, we prepare all atoms in the $s$-wave internal ground state $|g\rangle$ with motional ground state $\varphi_0(\mathbf{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 $\omega$, which is near-resonant with $\delta\equiv\omega-\omega_0\ll \omega,\omega_0$ [As a result, we will not distinguish $\omega$ and $\omega_0$ unless they combine into $\delta$.]. The electric field reads $\mathbf{E}_0(\mathbf r)=\mathbf{E}_0 e^{i\mathbf{k} \cdot \mathbf{r}}$ with $c|\mathbf{k}|=\omega$. 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

\begin{equation}
\delta H=-\sum_{n,i}\left(p^+_i(\mathbf{r_n})e^{-i\omega t}\mathbf{e}_i\cdot \mathbf{E}_0(\mathbf{r}_n)+\text{H.C.}\right).
\end{equation}

We assume the field strength $\mathbf{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

\begin{equation}
\mathbf{E}_{\text{tot}}(\omega,\mathbf{r})=\mathbf{E}_0(\mathbf{r})+\left<\mathbf{E}(\omega,\mathbf{r})\right>.
\end{equation}

Far from the atomic array, when only a single diffraction order exists, we expect

\begin{equation}
\mathbf{E}_{\text{tot}}(\omega,\mathbf{r})=\left(\mathbf{1}e^{ik_z z}+\mathbf{S}(\omega,\mathbf{k}_\parallel)e^{ik_z|z|}\right)\cdot \mathbf{E}_0 e^{i{\mathbf{k}}_\parallel\cdot \mathbf{r}_\parallel},
\end{equation}

and $\mathbf{S}(\omega,\mathbf{k}_\parallel)$ is the corresponding $\mathbf{S}$ matrix.

§.§ Diagrammatic Expansion

The contribution to the scattered light $\left<\mathbf{E}(\omega,\mathbf{r})\right>$ can be efficiently organized using the path-integral formulism. In particular, we work on the Keldysh contour <cit.>, 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

\begin{equation}\label{eq:Ep}
\begin{aligned}
\left<\mathbf{E(\omega,\mathbf{r})}\right>&=
\begin{tikzpicture}[scale=1.3, baseline={([yshift=-10pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[thick,snake arrow] (20pt,0pt) -- (-20pt,0pt);
\node[left] at (-10pt,10pt) {\scriptsize $\mathbf E(\omega,\mathbf{r})$};
\node[right] at (10pt,10pt) {\scriptsize $\mathbf{p}^-(\omega,\mathbf{r}_n)$};
%\draw[thick,mid arrow] (20pt,0pt) .. controls (8pt,14pt) and (-8pt,14pt) .. (-20pt,0pt);
\end{tikzpicture}
\\
&=-\sum_{n} \mathbf{G}_\text{R}^\mathbf{E}(\omega,\mathbf{r}-\mathbf{r}_n)\cdot \left<\mathbf{p}^-(\omega,\mathbf r_n)\right>.
\end{aligned}
\end{equation}

Here we use the wavy line to represent the propagation of photons. $\mathbf{G}_\text{R}^\mathbf{E}$ is the retarded Green's function matrix of $\mathbf{E}$ in free space defined as

\begin{equation}
\mathbf{G}_\text{R}^\mathbf{E}(t,\mathbf{r})\equiv-i\theta(t)\left<[E(t,\mathbf{r}),E(0,\mathbf{0})]\right>_{{d=0}}.
\end{equation}

In frequency and momentum space, we have

\begin{equation}\label{eq:GRtoG}
\begin{aligned}
\tilde{\mathbf{G}}_\text{R}^\mathbf{E}(\omega,\mathbf{k})&=(\epsilon_0\mathbf{1}+\epsilon_0\mathbf{k}\times \mathbf{k}\times/\omega^2)^{-1}=-\frac{\omega^2}{\epsilon_0}\tilde{\mathbf{G}}(\omega,\mathbf{k}).
\end{aligned}
\end{equation}

Here $\tilde{\mathbf{G}}(\omega,\mathbf{k})$ is the standard dyadic Green’s function <cit.>. 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 $\mathbf{p}^-$ is related to the incident light $\mathbf{E}_0$ by the Kubo formula <cit.>

\begin{equation}\label{eq:pE0}
\begin{aligned}
\left<\mathbf{p^-(\omega,\mathbf{r}_n)}\right>=&
\begin{tikzpicture}[scale=1.3, baseline={([yshift=-10pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[thick] (20pt,0.8pt) -- (-20pt,0.8pt);
\draw[thick] (20pt,-0.8pt) -- (-20pt,-0.8pt);
\node[left] at (-10pt,10pt) {\scriptsize $\mathbf p^-(\omega,\mathbf{r}_n)$};
\node[right] at (7pt,10pt) {\scriptsize $\mathbf{p}^+(-\omega,\mathbf{r}_m)$};
%\draw[thick,mid arrow] (20pt,0pt) .. controls (8pt,14pt) and (-8pt,14pt) .. (-20pt,0pt);
\end{tikzpicture}\\
=&-\sum_m~ \mathbf{G}_\text{R}^\mathbf{p}(\omega,\mathbf{r}_{nm})\cdot \mathbf{E}_0(\omega,\mathbf{r}_m).
\end{aligned}
\end{equation}

We use the double solid line for the retarded Green's function for dipole momentums $\mathbf{G}_\text{R}^\mathbf{p}(\omega,\mathbf{r})$ and $\mathbf{r}_{nm}\equiv\mathbf{r}_{n}-\mathbf{r}_{m}$. This is consistent with the semi-classical analysis <cit.>. The remaining task is to derive approximate formula for $\mathbf{G}_\text{R}^\mathbf{p}(\omega,\mathbf{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_\text{R}^{e_i}(\omega,\mathbf{r},\mathbf{r}')$ by emission and absorption of photons. As we will see, since the wave function for ground-state atoms is localized, only $G_\text{R}^{e_i}(\omega,\mathbf{r},\mathbf{r}')$ with $\mathbf{r}\approx \mathbf{r}\approx \mathbf{r}_n$ contributes to the light scattering. The bare Green's function near $\mathbf{r}_n=\mathbf{0}$ reads

\begin{equation}
G_\text{R}^{0,e_i}(q_0,\mathbf{r},\mathbf{r}')\approx \sum_a\frac{\varphi_{i,a}'(\mathbf{r})\varphi_{i,a}'(\mathbf{r}')^*}{q_0-\omega_0-\varepsilon'_{i,a}+i0^+}.
\end{equation}

The Schwinger-Dyson equation reads $(G_\text{R}^{e_i})^{-1}=(G_\text{R}^{0,e_i})^{-1}-\Sigma_\text{R}^{e_i}$, with the self-energy $\Sigma_\text{R}^{e_i}$

\begin{equation}\label{eq:Sigma}
\begin{aligned}
\Sigma_\text{R}^{e_i}(q_0,\mathbf{r},\mathbf{r}')=&
\begin{tikzpicture}[scale=1.5, baseline={([yshift=-1.8pt]current bounding box.center)}]
\filldraw (-12pt,0pt) circle (1pt);
\filldraw (12pt,0pt) circle (1pt);
\draw[thick,mid arrow] (24pt,0pt) -- (12pt,0pt);
\draw[thick,mid arrow] (-12pt,0pt) -- (-24pt,0pt);
\draw[thick,snake arrow] (12pt,0pt) -- (-12pt,0pt);
\node[right] at (14pt,-6pt) {\scriptsize $e_i$};
\node at (0pt,-6pt) {\scriptsize $\mathbf{E}$};
\draw[thick,mid arrow] (12pt,0pt) .. controls (4pt,9.5pt) and (-4pt,9.5pt) .. (-12pt,0pt);
\end{tikzpicture}\\
\approx&-\frac{\omega^2d^2}{\epsilon_0} G_{ii}(\omega,\mathbf{0})\sum_a(1-n_a)\varphi_a(\mathbf{r})\varphi_a^*(\mathbf{r}').
\end{aligned}
\end{equation}

Here $n_{a\geq 1}=0$ and $n_0=n$ is equal to the filling fraction. The appearance of $G_{ii}(\omega,\mathbf{0})=\mathbf{e}_i\cdot\mathbf{G}(\omega,\mathbf{0})\cdot \mathbf{e}_i$ owes to the approximation in (<ref>) by using $\mathbf{E}(\mathbf{r}_n)$ instead of $\mathbf{E}(\mathbf{r})$. The real-part of $\mathbf{G}(\omega,\mathbf{0})$ contributes to the lamb shift, which can be absorbed in the definition of $\omega_0$. As a result, we only keep the imaginary part $\mathbf{G}(\omega,\mathbf{0})=i\omega/6\pi\times \mathbf{1}$. We also assume $\delta,\varepsilon_{a},\varepsilon_{i,a}'\ll \omega$, and the resonance frequency $\omega$ is much larger than the loop frequency, which is an analogy of the Markovian approximation in the master equation <cit.>. 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 $\omega$ is known to be $\gamma =\omega_0^3d^2/3\pi\epsilon_0$. This leads to

\begin{equation}\label{eq:Sigmadelta}
\Sigma_\text{R}^{e_i}(q_0,\mathbf{r},\mathbf{r}')\approx-\frac{i\gamma}{2}\left[\delta(\mathbf{r}-\mathbf{r}')-n\varphi_0(\mathbf{r})\varphi_0^*(\mathbf{r}')\right],
\end{equation}

where we have used the completeness of local wave functions $\sum_a\varphi_a(\mathbf{r})\varphi_a^*(\mathbf{r}')=\delta(\mathbf{r}-\mathbf{r}')$.

Having obtained the dressed Green's function, we turn to the calculation of $\mathbf{G}_\text{R}^\mathbf{p}(\omega,\mathbf{r})$. Motivated by the standard Random Phase Approximation (RPA) in interacting fermions <cit.>, we consider the diagrams

\begin{equation}\label{eq:bubbles}
\begin{aligned}
\mathbf{G}_\text{R}^\mathbf{p}\ =&\
\begin{tikzpicture}[scale=1., baseline={([yshift=-3pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (20pt,0pt) .. controls (8pt,14pt) and (-8pt,14pt) .. (-20pt,0pt);
\draw[thick,mid arrow] (-20pt,0pt) .. controls (-8pt,-14pt) and (8pt,-14pt) .. (20pt,0pt);
\node[right] at (9pt,10pt) {\scriptsize $e_i$};
\node[right] at (9pt,-10pt) {\scriptsize $g$};
\node at (0pt,0pt) {\scriptsize $\delta_{mn}$};
\end{tikzpicture}
+\begin{tikzpicture}[scale=1., baseline={([yshift=-3pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (20pt,0pt) .. controls (8pt,14pt) and (-8pt,14pt) .. (-20pt,0pt);
\draw[thick,mid arrow] (-20pt,0pt) .. controls (-8pt,-14pt) and (8pt,-14pt) .. (20pt,0pt);
\node[right] at (9pt,10pt) {\scriptsize $e_i$};
\node[right] at (9pt,-10pt) {\scriptsize $g$};
\draw[thick,snake arrow] (40pt,0pt) -- (20pt,0pt);
\filldraw (40pt,0pt) circle (1pt);
\filldraw (80pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (80pt,0pt) .. controls (68pt,14pt) and (52pt,14pt) .. (40pt,0pt);
\draw[thick,mid arrow] (40pt,0pt) .. controls (52pt,-14pt) and (68pt,-14pt) .. (80pt,0pt);
\node[right] at (69pt,10pt) {\scriptsize $e_j$};
\node[right] at (69pt,-10pt) {\scriptsize $g$};
\node at (0pt,0pt) {\scriptsize $\mathbf{r_n}$};
\node at (60pt,0pt) {\scriptsize $\mathbf{r_m}$};
\end{tikzpicture}...
\end{aligned}
\end{equation}

Note that in our diagrammatic approach, $\mathbf{r}_n$ can be equal to $\mathbf{r}_m$, which is important as we will see later. The thick solid line represents the normalized Green's function $G_\text{R}^{e_i}$. The first bubble diagram, which is an elementary building block, is given by

\begin{equation}\label{eq:Pi}
\begin{aligned}
i[\bm{\Pi}_\text{R}(\omega)]_{ij}=&\frac{d^2}{2}\int \frac{d q_0}{2\pi}d\mathbf{r'}d\mathbf{r}~\left[G_\text{R}^{e_i}(q_0+\omega,\mathbf{r},\mathbf{r}')G_\text{K}^{g}(q_0,\mathbf{r}',\mathbf{r})+G_\text{K}^{e_i}(q_0+\omega,\mathbf{r},\mathbf{r}')G_\text{A}^{g}(q_0,\mathbf{r}',\mathbf{r})\right]\delta_{ij}.
\end{aligned}
\end{equation}

Here we write $\bm{\Pi}_\text{R}(\omega)$ as a diagonal matrix for later convenience. $G^\eta_\text{A}$ is the advanced Green's function. $G^\eta_{\mathbf{K}}=G^\eta_{\text{R}}\circ F_\eta-F_\eta\circ G^\eta_\text{A}$ is the Keldysh Green's function, and $F_\eta=(1-2n_\eta)$ is the quantum distribution function <cit.>. Here we use $\circ$ to represent the inner product of functions by integration. This leads to

\begin{equation}\label{eq:Pisim}
\begin{aligned}
\bm{\Pi}_{\text{R}}(\omega)_{ii}=&d^2n\int d\mathbf{r'}d\mathbf{r}~G_\text{R}^{e_i}(\omega+\varepsilon_0,\mathbf{r},\mathbf{r}')\varphi_0(\mathbf{r})^*\varphi_0(\mathbf{r}')=d^2n ~\varphi_0^*~\circ G_\text{R}^{e_i}\circ\varphi_0.
\end{aligned}
\end{equation}

It can be further simplified by noticing that

\begin{equation}
\begin{aligned}
\varphi_0^*~\circ G_\text{R}^{e_i}&\circ\varphi_0=\sum_a\frac{(\varphi_0\circ \varphi_{i,a}')^*~\varphi_{i,a}'\circ\varphi_0}{\delta+\varepsilon_0-\varepsilon_{i,a}'+\frac{i\gamma}{2}}
+\frac{i\gamma n}{2}\left(\sum_a\frac{(\varphi_0\circ \varphi_{i,a}')^*~\varphi_{i,a}'\circ\varphi_0}{\delta+\varepsilon_0-\varepsilon_{i,a}'+\frac{i\gamma}{2}}\right)^2+...
\end{aligned}
\end{equation}

As a result, we have

\begin{equation}\label{eq:pimain}
\begin{aligned}
\bm{\Pi}_{\text{R}}(\omega)_{ii}=&\frac{d^2n}{\pi_i(\omega)^{-1}-i\frac{\gamma n}{2}},\ \ \ \ \ \  \pi_i(\omega)\equiv\sum_a\frac{(\varphi_0\circ \varphi_{i,a}')^*~\varphi_{i,a}'\circ\varphi_0}{\delta+\varepsilon_0-\varepsilon_{i,a}'+\frac{i\gamma}{2}}.
\end{aligned}
\end{equation}

Then we can sum over the diagrams with multiple bubbles in (<ref>). This gives

\begin{equation}\label{eq:222}
\begin{aligned}
\end{aligned}
\end{equation}

Since the summation in (<ref>) is descrete, the Fourier transform here is defined as

\begin{equation}\label{eq:discreteFT}
\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{k}_\parallel)=\sum_n\bm{G}_\text{R}^\mathbf{E}(\omega,\mathbf{k}_\parallel)e^{-i\mathbf{k}_\parallel\cdot \mathbf{r}_n}.
\end{equation}

In particular, the denominator of (<ref>) is a generalization of the corresponding result under non-Hermitian Hamiltonians, which is $\omega\mathbf{1}-\mathcal{H}_{\text{eff}}$. As we will see later, (<ref>) 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 (<ref>), we obtain the relation between $\left<\mathbf p^-\right>$ and $\mathbf{E}_0$ in momentum space as

\begin{equation}\label{eq:alpha}
\begin{aligned}
\left<\mathbf p^-(\omega,\mathbf{k}_\parallel)\right>&=\bm{\alpha}(\omega,\mathbf{k}_\parallel)\cdot\mathbf{E}_0(\mathbf{k}_\parallel),\\
\bm{\alpha}(\omega,\mathbf{k}_\parallel)^{-1}&=-\bm{\Pi}_\text{R}(\omega)^{-1}+\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{k}_\parallel),
\end{aligned}
\end{equation}

Finally, for a single diffraction order, using (<ref>), $\bm{\alpha}$ is related to the $\bm{S}$ matrix as <cit.>

\begin{equation}\label{eq:S}
\begin{aligned}
\mathbf{S}(\omega,\mathbf{k}_\parallel)&=\frac{i\omega^2}{2a^2\epsilon_0k_z} \bm{\mathcal{P}}(\omega,{\mathbf{k}}_\parallel)\cdot \bm{\alpha}(\omega,\mathbf{k}_\parallel).
\end{aligned}
\end{equation}

Here $\mathcal{P}_{ij}(\omega,{\mathbf{k}}_\parallel)=\delta_{ij}-\xi_{ij}\frac{k_ik_j}{\omega^2}$. $\xi_{ij}=-1$ if only one of $(i,j)$ is in $z$ direction, and at the same time $z<0$. In other cases, $\xi_{ij}=1$.

Further simplification is possible for the normal incident case as in the experiment <cit.>, where we have $k_x=k_y=0$, $k_z=\omega$. Since, $\mathbf{E}_0$ lies in the $x$-$y$ plane, we have $\mathcal{P}=\mathbf{1}$. Moreover, due to the rotational symmetry, $\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{0})$ is also a diagonal matrix. Following the convention <cit.>, we define $\bm{\Delta}(\mathbf{k}_\parallel)=-\frac{3\pi\gamma}{\omega}\sum_{n\neq 0}\text{Re}~\bm{G}(\omega,\mathbf{r}_n)e^{-i\mathbf{r_n}\cdot \mathbf{k}_\parallel}$ and $\bm{\Gamma}(\mathbf{k}_\parallel)=\frac{6\pi\gamma}{\omega}\sum_{n\neq 0}\text{Im}~\bm{G}(\omega,\mathbf{r}_n)e^{-i\mathbf{r_n}\cdot \mathbf{k}_\parallel}$, which also become scalars $\Delta$ and $\Gamma$ in the $x$-$y$ plane for normal incident light. In particular, it is known that $\Gamma+\gamma=\gamma\frac{3\pi}{a^2\omega^2}$ <cit.>. As a result, the $\mathbf{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 $\Pi_\text{R}(\omega)\equiv\Pi_{\text{R},i_0i_0}(\omega)$, which is related to the local response $\pi(\omega)\equiv \pi_{i_0}(\omega)$. We also drop the $i_0$ index in $\varepsilon_{a}'\equiv\varepsilon_{i_0,a}'$ and $\varphi_{a}'\equiv\varphi_{i_0,a}'$ for conciseness. From now on, we focus on this normal incident case unless mentioned otherwise.

Using these definitions, we have

\begin{equation}\label{eq:S_gen}
\begin{aligned}
\pi(\omega)=\sum_a\frac{(\varphi_0\circ \varphi_{a}')^*~\varphi_{a}'\circ\varphi_0}{\delta+\varepsilon_0-\varepsilon_{a}'+\frac{i\gamma}{2}},\ \ \ \alpha&=-\frac{6\pi\epsilon_0 }{\omega^3}\frac{n\gamma/2}{\pi(\omega)^{-1}-n\Delta+in\Gamma/2}, \ \ \
\end{aligned}
\end{equation}

Since $\pi(\omega)$ is independent of $n$, the cooperative linewidth is linear in filling fraction $n$, consistent with the previous work <cit.>.

Equation (<ref>) and (<ref>) 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 <cit.>.

§ TRAPPING EFFECTS IN QUANTUM ATOMIC ARRAYS

In this section, we analyze (<ref>) and (<ref>) in several limits. We first consider the case with $V_{e_{i_0}}(\mathbf{r})=V_g(\mathbf{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 <cit.> [We thank the Referee for bringing this work to our attention.]. We then study the effect of trapping mismatch, which leads imperfectness of mirrors and multiple resonances. We finally go beyond (<ref>) and (<ref>) by analyzing the recoil effect.

§.§ Magic Wavelength

Numerical result for the fractional-filling effect with normal incident light with $\omega a_0=2\pi \times0.68$. Here we take $V_{e_{i_0}}(\mathbf{r})=V_g(\mathbf{r})$. (a). The reflection coefficient $R(\omega)$ as a function of detuning $\delta-n\Delta$ for different filling fraction $n$. (b). The transmission coefficient $T(\omega)$ as a function of detuning $\delta-n\Delta$ for different filling number $n$. (c). The filling-normalized absorptance $A$ and reflectance $\tilde{R}$, together with $T+R$ at cooperative resonance $\delta=n\Delta$, as a function of filling fraction $n$.

We begin with the special case $V_{e_{i_0}}(\mathbf{r})=V_g(\mathbf{r})$, which is valid when the optical lattice is formed by a light with the magic wavelength <cit.>.  In this case, only the $a=0$ term in (<ref>) contributes and the transition of internal state does not couple to the motional degree of freedom. Moreover, there is no dependence on the detailed shape of the potential. This leads to

\begin{equation}\label{eq:Magic_pi}
\begin{aligned}
\pi(\omega)\equiv \frac{1}{\delta+\frac{i\gamma}{2}},\ \ \ \ \ \  \alpha&=-\frac{6\pi\epsilon_0 }{\omega_0^3}\frac{n\gamma/2}{\delta-n\Delta+i(\gamma+n\Gamma)/2},\ \ \ \ \ \
\end{aligned}
\end{equation}

Here we have assumed the normal incidence for the probe light. The cooperative linewidth becomes $\gamma+n\Gamma$, consistent with  the experimental observation and numerical simulation in <cit.>, and the previous work <cit.>. 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

\begin{equation}
\begin{aligned}
T&=\frac{(\delta-n\Delta)^2+(1-n)^2\gamma^2/4}{(\delta-n\Delta)^2+(\gamma+n\Gamma)^2/4},\ \ \ \ \ \
\end{aligned}
\end{equation}

The filling-normalized absorptance $A=(1-T)/n$ and reflectance $\tilde{R}=R/n$ can the be computed straightforwardly.

We plot the numerical result for $\omega a_0=2\pi \times0.68$ as in the experiment <cit.> for various $n$ in Figure <ref>, where we have $\Delta/\gamma\approx 0.18$ and $\Gamma/\gamma\approx-0.48$. All above results reduces to the semi-classical results using spin operators when $n=1$, where the frequency shift is $\Delta$ and the linewidth becomes $\gamma+\Gamma$. On the other hand, for $n\rightarrow 0$, we get back to the single-atom response with natural linewidth $\gamma$.
As observed in the experiment <cit.>, generally, we have $T+R< 1$. This is due to the fact that the self-energy of the excited state (<ref>) contains the contribution of spontaneous emission of photons in arbitrary directions with random phases, which can not be observed by averaged $\mathbf{E}_{\text{tot}}$. However, the corresponding contribution exists if we measure energy density of electromagnetic field $\left<E^2(\mathbf{r})\right>$ <cit.>.
The filling-normalized absorptance $A$ show a weak dependence of $n$, while $\tilde{R}$ vanishes as $n \rightarrow 0$, consistent with the experimental observation and numerical simulation in <cit.>.

Numerical result for the trapping effect with normal incident light with $\omega a_0=2\pi \times0.68$. We fix $\omega^h_g=\gamma/4$ and consider $\omega^h_e<\gamma$. (a). The reflection coefficient $R(\omega)$ as a function of detuning for $n=1$ with different $\omega_e^h$. (b). The transmission coefficient $T(\omega)$ as a function of detuning for $n=1$ with different $\omega_e^h$. (c). The fitted $A_{\text{res}}$ and $\tilde{R}_{\text{res}}$ as a function of $n$ for different $\omega_e^h/\omega_g^h$. Here the dashed lines corresponds to $A_{\text{res}}$. (d). The fitted linewidth $\Gamma_{\text{cor}}$ as a function of $n$ for different $\omega_e^h$.

Finally we comment on the relation between our results and the Swchinger boson/fermion representation of spins <cit.>. Using (<ref>), we find the bubble reads $\Pi_\text{R}(\omega)={d^2n}/{(\delta+i\frac{\gamma(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\gamma n /2$ indeed cancels out with the corresponding contribution of the Green's function of photons $\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{k}_\parallel)$ in (<ref>) due to (Recall that the definition of $\Gamma$ does not contain $\mathbf{r_n}=0$, which is equal to $\gamma$.)

\begin{equation}
\bm{\Pi}_\text{R}(\omega)^{-1}-\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{k}_\parallel)=\frac{1}{d^2n}\left(\pi(\omega)^{-1}-\frac{i\gamma n}{2}+n\Delta+in\frac{\gamma+\Gamma}{2}\right)=\frac{1}{d^2n}\left(\pi(\omega)^{-1}+n\Delta+in\frac{\Gamma}{2}\right).
\end{equation}

This cancellation can be dated back to the cancellation between diagrams. Let's consider diagrams with one internal photons. Before contractions between $\psi_g$ and $\psi_g^\dagger$, it takes the form

\begin{equation}
\begin{tikzpicture}[scale=1., baseline={([yshift=-3pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[thick,mid arrow] (20pt,0pt) -- (-20pt,0pt);
%\draw[thick,mid arrow] (-20pt,0pt) .. controls (-8pt,-14pt) and (8pt,-14pt) .. (20pt,0pt);
\draw[thick,snake arrow] (60pt,0pt) -- (20pt,0pt);
\filldraw (60pt,0pt) circle (1pt);
\filldraw (100pt,0pt) circle (1pt);
\draw[thick,mid arrow] (100pt,0pt) -- (60pt,0pt);
%\draw[thick,mid arrow] (40pt,0pt) .. controls (52pt,-14pt) and (68pt,-14pt) .. (80pt,0pt);
\node at (0pt,10pt) {\scriptsize $\mathbf{r_n}, e_{i_0}$};
\node at (80pt,10pt) {\scriptsize $\mathbf{r_n}, e_{i_0}$};
\node at (-20pt,-10pt) {\scriptsize $\psi_g^\dagger(1)$};
\node at (20pt,-10pt) {\scriptsize $\psi_g(1)$};
\node at (60pt,-10pt) {\scriptsize $\psi_g^\dagger(2)$};
\node at (100pt,-10pt) {\scriptsize $\psi_g(2)$};
\end{tikzpicture}
\end{equation}

If we contract $\psi_g(1)$ with $\psi_g^\dagger(2)$ and $\psi_g(2)$ with $\psi_g^\dagger(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 $\psi_g^\dagger (2)$, there is already no occupation due to $\psi_g(2)$, we also need to take into account the contribution by contracting $\psi_g(1)$ with $\psi_g^\dagger(1)$ and $\psi_g(2)$ with $\psi_g^\dagger(2)$. More explicitly, we have $n=\left<\psi^\dagger\psi\psi^\dagger\psi\right>=\left<\psi^\dagger\psi\right>\left<\psi^\dagger\psi\right>+\left<\psi\psi^\dagger\right>\left<\psi^\dagger\psi\right>=n^2+n(1-n)=n$. However, the new diagram is just the bubble diagram, which has been taken into account in our diagrammatic expansion (<ref>). 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 $\Pi_\text{R}(\omega)$, and restrict the summation to $\mathbf{r}\neq 0$ in (<ref>). This is then consistent with rules in <cit.> using Schwinger particle 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 <cit.>. In this case, the corresponding contribution in the self-energy should be kept. We give an example in Appendix <ref>.

§.§ Trapping Mismatch

Numerical result for the trapping effect with normal incident light with $\omega a_0=2\pi \times0.68$. We fix $\omega^h_g=\gamma/4$ and consider $\omega^h_e\geq \gamma$ in (a-c). (a). The reflection coefficient $R(\omega)$ as a function of detuning for $n=1$ with different $\omega_e^h$. (b). The transmission coefficient $T(\omega)$ as a function of detuning for $n=1$ with different $\omega_e^h$. (c). The fitted linewidth $\Gamma_{\text{cor}}$ as a function of $n$ for different $\omega_e^h$. (d). The fitted center of the second peak $\Delta \omega$ as a function of $\omega_e^h$ for different filling fraction $n$. Here the dashed line corresponds to $\Delta \omega=2\omega_e^h$.

In this section, we discuss the effect of having $V_{e_{i_0}}(\mathbf{r})\neq V_g(\mathbf{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\rangle$ and excited-state atoms $|e_{i_0}\rangle$ have a trapping frequency $\omega^h_g$ and $\omega^h_e$ correspondingly. The motional ground state wave function $\varphi_0(\mathbf{r})$ reads

\begin{equation}
\varphi_0(\mathbf{r})=\left(\frac{\omega_g^h}{\pi}\right)^\frac{3}{4}e^{-\frac{\omega_g^h r^2}{2}}.
\end{equation}

Under this approximation, the analytical expression for $\pi(\omega)$ can be obtained by relating $\pi(\omega)$ to the single-particle propagator in harmonic traps. The results are presented in Appendix <ref>. It contains multiple resonances near $\delta=(3/2+2n)\omega_e^h-3\omega_g^h/2$, broadened by the natural lifetime $\gamma$ of the excited state. For $\omega_e^h \gtrsim \gamma$, this leads to different peaks in the spectral $-\text{Im}~\pi(\omega)/\pi$. For $\omega_e^h \lesssim \gamma$, different resonances merges, and only a single peak exists.

Numerical results for the effect of the trapping mismatch with normal incident light and $\omega a_0=2\pi \times0.68$. We fix small excited state trapping frequency $\omega^h_e=\gamma/20$. (a). The fitted $A_{\text{res}}$ as a function of $\omega_g^h$ for different $n$. (b). The fitted linewidth $\Gamma_{\text{cor}}$ as a function of $\omega_g^h$ for different $n$.

The parameters in the experiment <cit.> corresponds to $\omega_g^h<\gamma$, but at the same order $\sim$ MHz. We plot our results (<ref>) for different $\omega_e^h/\gamma$, $\omega_g^h/\gamma$ and $n$ in Figure <ref> and <ref>. We first fix $\omega_g^h/\gamma=1/4$ and study the effect of $\omega_e^h\neq \omega_g^h$ for small $\omega_e^h<\gamma$. As shown in Figure <ref> (a) and (b), both reflection coefficient $R(\omega)$ and transmission coefficient $T(\omega)$ show a single peak near $\delta-n\Delta-3(\omega_e^h-\omega_g^h)/2=0$. For either $\omega_e^h>\omega_g^h$ or $\omega_g^h>\omega_e^h$, the atomic mirror becomes imperfect with $\text{max}~R<1$ and $\text{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(\omega)$ near $\delta=n\Delta+3(\omega_e^h-\omega_g^h)/2$ as

\begin{equation}
\end{equation}

and define $T_{\text{res}}=T(n\Delta+3(\omega_e^h-\omega_g^h)/2+\delta_0)$. $\tilde{R}_{\text{res}}$ and $A_{\text{res}}$ can then be computed correspondingly using $R_{\text{res}}$ and $T_{\text{res}}$. The numerical results in (c-d) show $\tilde{R}_{\text{res}}$ and $A_{\text{res}}$ also decreases when $\omega_e^h\neq \omega_g^h$. The cooperative linewidth $\Gamma_{\text{cor}}$ is linear in $n$, with similar slope for different $\omega_e^h<\gamma$.

We then consider larger $\omega_e^h\geq \gamma$ in Figure <ref>. Now as shown in Figure <ref> (a) and (b), multiply peaks appear in both reflection coefficient $R(\omega)$ and transmission coefficient $T(\omega)$. The center of peaks locates near energy $2n\omega_e^h$, where the transition from $|g\rangle$ to $|e_{i_0}\rangle$ is accompanied with the excitation of motional degree of freedom. As an example, we fit the position of the second peak $\Delta \omega$, and plot it as a function of $\omega_e^h$ in (d). Similar to the small $\omega_e^h$ case, the cooperative linewidth $\Gamma_{\text{cor}}$ is still linear in $n$. However, their slope show dependence on $\omega^h_e$.

We finally study $A_{\text{res}}$ and $\Gamma_{\text{cor}}$ as a function of $\omega_g^h$. We fix a small $\omega_e^h=\gamma/20$, as an analogy of the anti-trapped excited state in experiment <cit.>, and tune $\omega_g^h$. As shown in Figure <ref>, we find when $\omega_g^h$ becomes larger, the absorption rate decreases and the cooperative linewidth becomes larger. This is due to the increase of the trapping mismatch for small $\omega_e^h$. For small $\omega_g^h$, the decrease in $A_{\text{res}}$ and the increase of the decay rate show quadratic dependence, while for large $\omega_g^h$, the dependence becomes linear. This is a close analogy of the experimental observation in <cit.>.

§.§ Recoil Effects

Now we go beyond the limit of $\sigma \ll a_0$ and consider corrections to the leading order of $\eta=\sigma/a_0$ due to the recoil of atoms. The recoil effect has been discussed in previous works <cit.> 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 (<ref>) without the approximation (<ref>) as explained in the Appendix <ref>. However, the same result can also be obtained using direct physical intuition. Our result (<ref>) for $\pi(\omega)$ has a simple physical meaning: it measures the optical response of a single particle in the harmonic trap, with a lift time $\gamma$ for the excited state particle. Without the recoil effect, we take the inner product between wavefunctions $\varphi_0(\mathbf{r})$ and $\varphi_a'(\mathbf{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 $\mathbf{k}_1$ and emitting a photon with momentum $\mathbf{k}_2$. Physically, we expect the local response $\pi(\omega)$ takes the form

\begin{equation}\label{eq:pik}
\begin{aligned}
\pi_{\mathbf{k}_1\mathbf{k}_2}(\omega)
&=\sum_a\int d\mathbf{r}d\mathbf{r}'~\varphi_0(\mathbf{r})^*\varphi_a'(\mathbf{r})e^{-i\mathbf{k}_2\cdot \mathbf{r}}\frac{1}{\delta+\varepsilon_0-\varepsilon_{a}'+\frac{i\gamma}{2}}e^{i\mathbf{k}_1\cdot \mathbf{r}'}\varphi_a'(\mathbf{r}')^*\varphi_0(\mathbf{r}')
\end{aligned}
\end{equation}

Since the photons can propagate in any direction within $x-y$ plane, we further average over the direction of the momentum $\mathbf{k}_i=k (\cos\theta_i,\sin\theta_i,0)$:

\begin{equation}\label{eq:recoil}
\begin{aligned}
\pi(\omega)&=\sum_a\int d\mathbf{r}d\mathbf{r}'\frac{d\theta_1}{2\pi}\frac{d\theta_2}{2\pi}~\varphi_0(\mathbf{r})^*\varphi_a'(\mathbf{r})e^{-i\mathbf{k}_2\cdot \mathbf{r}}\frac{1}{\delta+\varepsilon_0-\varepsilon_{a}'+\frac{i\gamma}{2}}e^{i\mathbf{k}_1\cdot \mathbf{r}'}\varphi_a'(\mathbf{r}')^*\varphi_0(\mathbf{r}')
\end{aligned}
\end{equation}

We focus on the case with $\omega_g^h=\omega_e^h$. For small $\eta=\sigma/a_0=k/\sqrt{\omega_g^h}$, we can expand $\pi(\omega)$ to obtain

\begin{equation}
\pi(\omega)=\frac{1}{ \delta +i \gamma/2 }-\frac{\eta^2/2}{ \delta +i \gamma/2 }+\frac{1}{8} \eta^4 \left(\frac{1}{ \delta +i \gamma/2 }+\frac{1/2}{  \delta -2 \omega_g^h+i \gamma/2}\right)+O(\eta^6).
\end{equation}

To the leading order $O(\eta^2)$, we find only a renormalization of the residue near $\delta+i\gamma/2=0$. At the sub-leading order, $O(\eta^4)$ we see a non-trivial contribution where atoms are excited due to the recoil even without trapping mismatch $\omega^h_e=\omega^h_g$, because the recoil increases the energy of atoms. However, when we further plot the reflection and transmission coefficients (see Figure <ref>), we find no visible peaks at $\delta-n\Delta=2 \omega^h_e$ even for $\omega_g^h>\gamma$ due to the suppression of $\eta^4$. We find $R(\omega)$ and $T(\omega)$ only show weak dependence on $\eta$ for small $\eta$. In particular, the parameter regime in the experiment <cit.> corresponds to $\eta\sim 0.1$, which is almost the same as $\eta=0$. This justifies our discussions in previous sections.

Numerical results for the recoil effect with normal incident light and $\omega a_0=2\pi \times0.68$. We fix small excited state trapping frequency $\omega^h_g=2\gamma$ and vary $\eta=\sigma/a_0$. (a). The reflection coefficient $R(\omega)$ as a function of detuning for $n=1$ with different $\eta$. (b). The transmission coefficient $T(\omega)$ as a function of detuning for $n=1$ with different $\eta$.

§ 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 $\omega^g_e \sim \gamma$. 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 <cit.>. Recently, there are also experimental studies on the Pauli blocking of light scattering in degenerate fermions <cit.>. The diagrammatic approach developed here can also be applied to study the optical response of degenerate fermion gases.

§ ACKNOWLEDGMENTS

We especially thank Yu Chen and Jianwen Jie for helpful discussions. We thank the Referee for bringing several related works to our attention, and the suggestion of studying the recoil effect.
P.Z. acknowledges support from the Walter Burke Institute for Theoretical Physics at Caltech.

§ HOMOGENEOUS MIRRORS

Now we consider the case where the atom gas is homogeneous in $x-y$ plane at $z=0$, similarly to <cit.>. The model of the full system reads

\begin{equation}
\begin{aligned}
S=&\int dtd\mathbf{r}~ \sum_{i}\psi^\dagger_{e_i}(t,\mathbf{r})\left( i\partial_t-\omega_0 +\frac{\nabla^2}{2}\right)\psi_{e_i}(t,\mathbf{r})+\psi^\dagger_{g}(t,\mathbf{r})\left( i\partial_t +\frac{\nabla^2}{2}\right)\psi_{g}(t,\mathbf{r})\\
&+\int dtd\mathbf{r}~\sum_{i}|d|\left(\psi^\dagger_{e_i}(t,\mathbf{r})\psi_{g}(t,\mathbf{r})+\psi^\dagger_{g}(t,\mathbf{r})\psi_{e_i}(t,\mathbf{r})\right)\mathbf{e}_i\cdot \mathbf{E}(t,\mathbf{r}).
\end{aligned}
\end{equation}

The first step is again to determine the renormalization of the excited state Green's function. The self-energy reads

\begin{equation}
\begin{aligned}
\Sigma^{e_i}_{\text{R}}(\tilde{\omega},\mathbf{q}_\parallel)_{ij}=&\begin{tikzpicture}[scale=1.5, baseline={([yshift=-1.8pt]current bounding box.center)}]
\filldraw (-12pt,0pt) circle (1pt);
\filldraw (12pt,0pt) circle (1pt);
\draw[thick,mid arrow] (24pt,0pt) -- (12pt,0pt);
\draw[thick,mid arrow] (-12pt,0pt) -- (-24pt,0pt);
\draw[thick,snake arrow] (12pt,0pt) -- (-12pt,0pt);
\node[right] at (10pt,-6pt) {\scriptsize $e$};
\node at (0pt,-6pt) {\scriptsize $\mathbf{E}$};
\draw[thick,mid arrow] (12pt,0pt) .. controls (4pt,9.5pt) and (-4pt,9.5pt) .. (-12pt,0pt);
\end{tikzpicture}
=-\int \frac{d\mathbf{k}'}{(2\pi)^2}\frac{\omega^2d^2}{\epsilon_0}\mathbf{e}_i\cdot \mathbf{G}(\omega,\mathbf{q}_\parallel-\mathbf{k}')\cdot\mathbf{e}_j (1-n_F(\varepsilon_{k'}))\\
=&-i\frac{\gamma}{2}+\int \frac{d\mathbf{k}'}{(2\pi)^2}\frac{\omega^2d^2}{\epsilon_0}\mathbf{e}_i\cdot \mathbf{G}(\omega,\mathbf{q}_\parallel-\mathbf{k}')\cdot\mathbf{e}_j n_F(\varepsilon_{k'})\equiv-i\frac{\gamma}{2}\delta_{ij}+\delta \mathcal{H}(\mathbf{q}_\parallel)_{ij}.
\end{aligned}
\end{equation}

This leads to

\begin{equation}
G^{e_i}_{\text{R}}(\tilde{\omega},\mathbf{q})=\frac{1}{\tilde \omega-\omega_0-\varepsilon_q-\Sigma^{e_i}_{\text{R}}(\tilde{\omega},\mathbf{q}_\parallel)}=\frac{1}{(\tilde \omega-\omega_0-\varepsilon_q+i\frac{\gamma}{2})\mathbf{1}-\delta \mathcal{H}(\mathbf{q})}.
\end{equation}

Here $\varepsilon_q=q^2/2$. We then compute the $\Pi_\text{R}$. The result reads

\begin{equation}
\bm{\Pi}_\text{R}(\omega,\mathbf{k}_\parallel)=\int\frac{d\mathbf{q}}{(2\pi)^2}\frac{d^2n_F(\varepsilon_q)}{(\delta+\varepsilon_q-\varepsilon_{q+k}+i\gamma/2)\mathbf{1}-\delta \mathcal{H}(\mathbf{q}+\mathbf{k}_\parallel)}.
\end{equation}

Summing up the contribution from photons, we find

\begin{equation}
\bm{\alpha}(\omega,\mathbf{k}_\parallel)=\frac{1}{-\bm{\Pi}_\text{R}(\omega,\mathbf{k}_\parallel)^{-1}-{ \omega^2} \mathbf{G}(\omega,\mathbf{k}_\parallel,z=0)/\epsilon_0}.
\end{equation}

Here since the system is homogeneous, the Fourier transform is

\begin{equation}
\mathbf{G}(\omega,\mathbf{k}_\parallel,z)=\int d\mathbf{r}_\parallel~\mathbf{G}(\omega,\mathbf{r})e^{-i\mathbf{k}_\parallel\cdot \mathbf{r}_\parallel}=\frac{i}{2k_z}e^{ik_z|z|}\bm{\mathcal{P}}(\omega,\mathbf{k}_\parallel).
\end{equation}

The relation between $\mathbf{E}_{\text{tot}}$ and $\alpha$ becomes

\begin{equation}
\begin{aligned}
\mathbf{E}_{\text{tot}}&=\mathbf{E}_0(\mathbf{r})+\frac{\omega^2}{\epsilon_0}\int d\mathbf{r}' \mathbf{G}(\omega,\mathbf{r}-\mathbf{r}')\cdot \bm{\alpha}(\omega,\mathbf{k}_\parallel)\cdot \mathbf{E}_0 e^{i\mathbf{k}_\parallel \cdot \mathbf{r}'}\\
&=\mathbf{E}_0(\mathbf{r})+\frac{\omega^2}{\epsilon_0} \mathbf{G}(\omega,\mathbf{k}_\parallel,z)\cdot \bm{\alpha}(\omega,\mathbf{k}_\parallel)\cdot \mathbf{E}_0 e^{i\mathbf{k}_\parallel \cdot \mathbf{r}}.
\end{aligned}
\end{equation}

The $S$ matrix reads

\begin{equation}
\mathbf{S}(\omega,\mathbf{k}_\parallel)=\frac{i\omega^2}{2\epsilon_0 k_z}\bm{\mathcal{P}}(\omega,\mathbf{k}_\parallel)\cdot \bm{\alpha}(\omega,\mathbf{k}_\parallel).
\end{equation}

Again we consider the normal incident light. We further assume the density of the system is low as in <cit.>. We have

\begin{equation}
\begin{aligned}
\Pi_\text{R}(\omega,\mathbf{0})_{ij}=&\int\frac{d\mathbf{q}}{(2\pi)^2}\frac{d^2n_F(\varepsilon_q)}{(\delta+i\gamma/2)\mathbf{1}-\delta \mathcal{H}(\mathbf{q})}\\
=&\frac{n_{\text{2D}}d^2}{\delta+i\gamma/2}\delta_{ij}+\frac{\omega^2d^4}{\epsilon_0(\delta+i\gamma/2)^2}\int\frac{d\mathbf{q}d\mathbf{q}'}{(2\pi)^4}\left[\mathbf{e}_i\cdot \mathbf{G}(\omega,\mathbf{q}-\mathbf{q}',z=0)\cdot\mathbf{e}_j\right]n_F(\varepsilon_{q'})n_F(\varepsilon_{q})\\
=&\left(\frac{n_{\text{2D}}d^2}{\delta+i\gamma/2}+\delta \Pi_\text{R}(\omega)\right)\delta_{ij}.
\end{aligned}
\end{equation}

We find

\begin{equation}
\alpha(\omega,\mathbf{0})=\frac{1}{-\frac{\delta+i\gamma/2}{d^2n_{\text{2D}}}+\frac{(\delta+i\gamma/2)^2}{d^4n_{\text{2D}}^2}\delta \Pi_\text{R}(\omega)-\frac{i{ \omega}}{2\epsilon_0}}, \ \ \ \ \ \ \ \ \ \ S(\omega)=\frac{\frac{i{ \omega}}{2\epsilon_0}}{-\frac{\delta+i\gamma/2}{d^2n_{\text{2D}}}+\frac{(\delta+i\gamma/2)^2}{d^4n_{\text{2D}}^2}\delta \Pi_\text{R}(\omega)-\frac{i{ \omega}}{2\epsilon_0}}.
\end{equation}

Here we have

\begin{equation}
\begin{aligned}
\frac{(\delta+i\gamma/2)^2}{d^4n_{\text{2D}}^2}\delta \Pi_\text{R}(\omega)
=\frac{\omega^2}{\epsilon_0n_{\text{2D}}^2}\int\frac{d\mathbf{q}d\mathbf{q}'}{(2\pi)^4}\left[\mathbf{e}_i\cdot \mathbf{G}(\omega,\mathbf{q}-\mathbf{q}',z=0)\cdot\mathbf{e}_i\right]n_F(\varepsilon_{q'})n_F(\varepsilon_{q}),
\end{aligned}
\end{equation}

As a result, . This takes the similar form as results in <cit.> for the 3D case. This is the contribution from the density-density correlation in free fermion gases. Finally, we have

\begin{equation}
S(\omega)=\frac{\frac{i{ \omega}}{2\epsilon_0}d^2n_{\text{2D}}}{-\delta-i\gamma/2+\frac{(\delta+i\gamma/2)^2}{d^2n_{\text{2D}}}\delta \Pi_\text{R}(\omega)-\frac{i{ \omega}}{2\epsilon_0}d^2n_{\text{2D}}},
\end{equation}

which means $\delta \Pi_\text{R}(\omega)$ effectively shifts the resonant energy and the decay rate.

§ THE ANALYTICAL FORMULA FOR TEXT

In this Appendix, we present detailed derivation of the analytical formula for $\pi(\omega)$. We trick is to use the transformation to the time domain

\begin{equation}
\begin{aligned}
\pi(\omega)&=\sum_a\int d\mathbf{r}d\mathbf{r}'~\varphi_0(\mathbf{r})^*\varphi_a'(\mathbf{r})\frac{1}{\delta+\varepsilon_0-\varepsilon_{a}'+\frac{i\gamma}{2}}\varphi_a'(\mathbf{r}')^*\varphi_0(\mathbf{r}')\\
&=-\sum_a\int d\mathbf{r}d\mathbf{r}'\int_0^\infty d\tau~\varphi_0(\mathbf{r})^*\varphi_a'(\mathbf{r})e^{(\delta+\varepsilon_0-\varepsilon_{a}'+\frac{i\gamma}{2})\tau}\varphi_a'(\mathbf{r}')^*\varphi_0(\mathbf{r}')\\
&=-\int d\mathbf{r}d\mathbf{r}'\int_0^\infty d\tau~e^{(\delta+\varepsilon_0+\frac{i\gamma}{2})\tau}\varphi_0(\mathbf{r})^*K_{\omega^h_e}(\tau,\mathbf{r},\mathbf{r}')\varphi_0(\mathbf{r}').
\end{aligned}
\end{equation}

Here we have assumed the integral over $\tau$ is convergent by restricting the $\delta+\varepsilon_0<3\omega_e^h/2$. After the integration, analytical continuation can be applied to release this restriction. Here $K_{\omega^h_e}(\tau,\mathbf{r},\mathbf{r}')$ is the imaginary time Green's function in a harmonic trap with trapping frequency $\omega^h_e$. We have

\begin{equation}
K_{\omega^h_e}(\tau,\mathbf{r},\mathbf{r}')=\left(\frac{\omega^h_e}{2\pi\sinh \omega^h_e\tau}\right)^{\frac{3}{2}}\exp\left(-\frac{\omega^h_e}{2}\left[(r^2+r'^2)\coth\omega^h_e\tau-\frac{2\mathbf{r}\cdot\mathbf{r}'}{\sinh \omega^h_e\tau}\right]\right).
\end{equation}

The integral over $\mathbf{r}$ and $\mathbf{r}'$ can be carried out first. We find

\begin{equation}
\begin{aligned}
\pi(\omega)&=-2\sqrt{2}\int_0^\infty d\tau~e^{a_0\tau}\left(\frac{\omega^h_e\omega^h_g}{2\omega^h_e\omega^h_g\cosh \omega^h_e\tau+[(\omega^h_e)^2+(\omega^h_g)^2]\sinh \omega^h_e\tau}\right)^{\frac{3}{2}}.
\end{aligned}
\end{equation}

Here we have defined $a_0=\left(\delta+\frac{3\omega^h_g}{2}+\frac{i\gamma}{2}\right)$ for conciseness. Then the integral over $\tau$ gives

\begin{equation}
\begin{aligned}
\pi(\omega)&=\frac{\sqrt{2}}{\omega^h_e}\left(\frac{2\omega^h_e\omega^h_g}{(\omega^h_e+\omega^h_g)^2}\right)^{\frac{3}{2}}\frac{q^{-p-1} \left((-2 p (q-1)-q+2) B_q\left(p+1,\frac{1}{2}\right)-(2 p+3)
\end{aligned}
\end{equation}

where $p\equiv -\frac{2a_0+\omega^h_e}{4\omega^h_e}$ and $q\equiv \frac{(\omega^h_e+\omega^h_g)^2}{(\omega^h_e-\omega^h_g)^2}$ . $B_z(a,b)$ is the incomplete beta function defined as $B_z(a,b)=\int_0^z t^{a-1}(1-t)^{b-1}dt$. For $\omega^h_e=\omega^h_g$, one can check that above result can be simplified as $\pi(\omega)^{-1}=a_0-3\omega^h_e/2$.

§ THE DERIVATION OF THE RECOIL EFFECTS

Here we give the derivation of (<ref>) using the diagrammatic approach. For simplicity, we directly take diagrams under the rule of Schwinger bosons/fermions as discussed section <ref>. As a result, the self-energy of excited state $\Sigma_\text{R}^{e_i}(q_0)$ is just a constant $-i\gamma /2$, and the summation in (<ref>) is restricted by $\mathbf{r}_m \neq \mathbf{r}_n$.

To determine $\pi(\omega)$, we now examine a single diagram

\begin{equation}
\begin{tikzpicture}[scale=1., baseline={([yshift=-3pt]current bounding box.center)}]
\filldraw (-20pt,0pt) circle (1pt);
\filldraw (20pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (20pt,0pt) .. controls (8pt,14pt) and (-8pt,14pt) .. (-20pt,0pt);
\draw[thick,mid arrow] (-20pt,0pt) .. controls (-8pt,-14pt) and (8pt,-14pt) .. (20pt,0pt);
\node[right] at (9pt,10pt) {\scriptsize $e_i$};
\node[right] at (9pt,-10pt) {\scriptsize $g$};
\draw[thick,snake arrow] (40pt,0pt) -- (20pt,0pt);
\filldraw (40pt,0pt) circle (1pt);
\filldraw (80pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (80pt,0pt) .. controls (68pt,14pt) and (52pt,14pt) .. (40pt,0pt);
\draw[thick,mid arrow] (40pt,0pt) .. controls (52pt,-14pt) and (68pt,-14pt) .. (80pt,0pt);
\node[right] at (69pt,10pt) {\scriptsize $e_j$};
\node[right] at (69pt,-10pt) {\scriptsize $g$};
\node at (0pt,0pt) {\scriptsize $\mathbf{r_n}$};
\node at (60pt,0pt) {\scriptsize $\mathbf{r_m}$};

\draw[thick,snake arrow] (100pt,0pt) -- (80pt,0pt);
\filldraw (100pt,0pt) circle (1pt);
\filldraw (140pt,0pt) circle (1pt);
\draw[line width=0.5mm,mid arrow] (140pt,0pt) .. controls (128pt,14pt) and (112pt,14pt) .. (100pt,0pt);
\draw[thick,mid arrow] (100pt,0pt) .. controls (112pt,-14pt) and (128pt,-14pt) .. (140pt,0pt);
\node[right] at (129pt,10pt) {\scriptsize $e_k$};
\node[right] at (129pt,-10pt) {\scriptsize $g$};
\node at (120pt,0pt) {\scriptsize $\mathbf{r_p}$};
\draw[thick,snake arrow] (-20pt,0pt) -- (-40pt,0pt);
\draw[thick,snake arrow] (160pt,0pt) -- (140pt,0pt);
\end{tikzpicture}
\end{equation}

This corresponds to

\begin{equation}
\begin{aligned}
i\int d\mathbf{r}_i\sum_{m}\bm{\Pi}_\text{R}(\omega,\mathbf{r}',\mathbf{r}_3)\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{nm}+\mathbf{r}_3-\mathbf{r}_2)\bm{\Pi}_\text{R}(\omega,\mathbf{r}_3,\mathbf{r}_2)\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{mp}+\mathbf{r}_2-\mathbf{r}_1)\bm{\Pi}_\text{R}(\omega,\mathbf{r}_1,\mathbf{r}),
\end{aligned}
\end{equation}

where we have

\begin{equation}
\begin{aligned}
\bm{\Pi}_\text{R}(\omega,\mathbf{r},\mathbf{r}'')_{ii}=d^2n~\pi_i(\omega,\mathbf{r},\mathbf{r}'')
&=\sum_a\varphi_0(\mathbf{r})^*\varphi_{i,a}'(\mathbf{r})\frac{d^2 n}{\delta+\varepsilon_0-\varepsilon_{i,a}'+\frac{i\gamma}{2}}\varphi_{i,a}'(\mathbf{r}')^*\varphi_0(\mathbf{r}').
\end{aligned}
\end{equation}

We have separated the integral over the full space into a summation over lattice sites $\mathbf{r}_n$, and an integral near each sites $\mathbf{r}_i$. For small $\sigma\ll a_0$, the dominate contribution comes from $\mathbf{r}_{mp}\gg \mathbf{r}_2-\mathbf{r}_1$ and $\mathbf{r}_{nm}\gg \mathbf{r}_3-\mathbf{r}_2$. Moreover, for normal (or nearly normal) incident light, we are probing the system with small $\mathbf{k}_\parallel$, which comes from contributions at large $\mathbf{r}_{nm}$ and $\mathbf{r}_{mp}$. We then expand $\bm{\tilde{\mathcal{G}}}_\text{R}^\mathbf{E}$ and take the standard approximation at long distance <cit.>:

\begin{equation}
\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{mp}+\mathbf{r}_2-\mathbf{r}_1)\approx\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{mp})\exp(ik \hat{\mathbf{r}}_{mp}\cdot(\mathbf{r}_2-\mathbf{r}_1)).
\end{equation}

Using this expression, we find

\begin{equation}
\begin{aligned}
i\int d\mathbf{r}_3d\mathbf{r}_1\sum_{m}\bm{\Pi}_\text{R}(\omega,\mathbf{r}',\mathbf{r}_3)e^{ik \hat{\mathbf{r}}_{nm}\cdot\mathbf{r}_3}\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{nm})\bm{\Pi}_\text{R}^{nmp}(\omega)\bm{{\mathcal{G}}}_\text{R}^\mathbf{E}(\omega,\mathbf{r}_{mp})e^{-ik \hat{\mathbf{r}}_{mp}\cdot\mathbf{r}_1}\bm{\Pi}_\text{R}(\omega,\mathbf{r}_1,\mathbf{r})
\end{aligned}
\end{equation}

Here we have defined

\begin{equation}
\bm{\Pi}_\text{R}^{nmp}(\omega)_{ii}=\sum_a\int d\mathbf{r}d\mathbf{r}'~\varphi_0(\mathbf{r})^*\varphi_{i,a}'(\mathbf{r})e^{-ik\hat{\mathbf{r}}_{nm}\cdot\mathbf{r}}\frac{d^2 n}{\delta+\varepsilon_0-\varepsilon_{i,a}'+\frac{i\gamma}{2}}e^{ik\hat{\mathbf{r}}_{mp}\cdot\mathbf{r}'}\varphi_{i,a}'(\mathbf{r}')^*\varphi_0(\mathbf{r}')
\end{equation}

This already takes the form of (<ref>), with the direction of the photons determined by $\hat{\mathbf{r}}_{nm}$ and $\hat{\mathbf{r}}_{mp}$. Note that due to the rotation symmetry of the isotropic harmonic trap, it only depends on $\hat{\mathbf{r}}_{nm} \cdot \hat{\mathbf{r}}_{mp}$. Finally, we need to perform the Fourier transform by summing up $m,n,p$. For large ${\mathbf{r}}_{nm}$ and ${\mathbf{r}}_{mp}$, we approximate this summation as a average over $\hat{\mathbf{r}}_{nm} \cdot \hat{\mathbf{r}}_{mp}$. Finally, we obtain the “naive” formula (<ref>) by focusing on incident light polarized in the direction of $i_0$.