Introduction

When light interacts with an ensemble of atoms, a spontaneous emission with an enhanced rate emerges due to strong atom-atom correlations induced by photons rescattering in the medium. This is so-called superradiance1,2,3, which can be manifested by resonant dipole-dipole interactions (RDDI)4,5,6. The RDDI feature a long-range effect that connects every pair of the atoms, and thus enable the light-matter interacting system collectively coupled. This collective coupling is responsible for coherent light scattering7,8,9,10,11,12,13,14. On the other hand in a dense medium with RDDI, an emission with a reduced decay rate, the subradiance15,16,17, can also show up as an afterglow which follows the initialized superradiance in a cloud17. Many other proposals to prepare and manipulate the subradiant states can be found in free-standing atomic arrays18,19,20,21,22,23,24,25,26,27,28, atoms in a cavity29, and metamolecules30.

In addition to RDDI in a free space or three-dimensional (3D) reservoir, this near-resonantly driven dipole-dipole interaction can be also established in a confined one-dimensional (1D) atom-nanophotonic waveguide system31,32,33,34,35, which shows infinite-range couplings in sinusoidal forms and is observed recently in two atomic clouds above the optical nanofibers36. By structuring 1D reservoir37,38,39, the confined 1D system can further allow an effective uni-directional decay channel, which breaks the time reversal symmetry, and construct a chiral quantum optical network40,41. In such confined systems, the guided mode of light can be stored via electromagnetically induced transparency42, and long-range quantum magnetism can be simulated with tunable spin-spin interactions43 where a pseudo spin is presented by two atomic internal states.

The RDDI involve the coherent and dissipative parts, which respectively determine the frequency shift and line width of the radiation. For the dissipative parts, Dicke’s regime can be reached in short distance of \(\xi \ll 1\) where ξ represents a dimensionless scale of mutual separation. In opposite regime of long distance that \(\xi \gg 1\) for both coherent and dissipative parts, 1D RDDI sustain the interaction strength as cos(ξ) and sin(ξ) respectively, in contrast to \( \sim \,\mathrm{1/}\xi \) of 3D RDDI. Another stark contrast of 1D and 3D RDDI lies in the short distance asymptotics of the coherent parts, which are \( \sim \xi \) and \( \sim \mathrm{1/}{\xi }^{3}\) respectively. The divergence of 3D RDDI in the limit of ξ → 0 indicates a breakdown of quantum optical treatment, which means a lack of complete and genuine description of induced dipoles in short distance. On the other hand in 1D confinement, ξ in general can not be made infinitely small, while its interaction strength of sin(ξ) already diminishes as \(\xi \ll 1\). Other than these well-known 3D and 1D RDDI, two-dimensional (2D) RDDI is less explored since it requires intricately designed 2D reservoir. Here we consider a setting of 2D optical lattices of atoms inside a planar cavity with perfect mirrors. This way, the atoms can only dissipate by spontaneous emissions into this 2D reservoir. Other potential settings for this confined 2D dissipative bath can be a 2D lattice of coupled ring resonators44 or 2D arrays of superconducting artificial atoms. We note that 2D scattering by point scatterers has been explored using 2D coupled dipole equations45, and single-photon superradiance can be engineered by modeling the distance-dependence of RDDI in various dimensional electromagnetic reservoirs46.

In this article, we show the distinctive feature of 2D RDDI, where different atomic polarizations display significantly distinct long range behavior. We investigate the collective radiation properties under a symmetric state of single photon excitation. This state also allows subradiant decays at some specific or much longer ξ, showing long-range atom-atom correlations. We further study the time dynamics of phase-imprinted subradiant states by applying spatially dependent phases on the 2D atomic arrays. The radiation properties of these potentially controllable single-excitation subradiant states highly depend on the 2D lattice structures. Thus, it allows engineering of light-matter interactions, and promises applications in quantum light storage and state manipulations.

Collective Properties from RDDI in a Confined Two-Dimensional Reservoir

RDDI in a confined two-dimensional reservoir

We follow the general formalism of RDDI in 3D free space5, with more details in Methods, and its Hamiltonian reads

$$H=\sum _{\mu =1}^{N}\,\hslash {\omega }_{e}{\hat{\sigma }}_{\mu }^{\dagger }{\hat{\sigma }}_{\mu }-\sum _{\mu =1}^{N}\,\sum _{q}\,{g}_{q}({e}^{i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }-i{\omega }_{q}t}{\hat{a}}_{q}+{e}^{-i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }+i{\omega }_{q}t}{\hat{a}}_{q}^{\dagger })({\hat{\sigma }}_{\mu }+{\hat{\sigma }}_{\mu }^{\dagger }),$$
(1)

where the atomic raising and lowering operators are \({\hat{\sigma }}_{\mu }^{\dagger }\equiv |e{\rangle }_{\mu }\langle g|\) and \({\hat{\sigma }}_{\mu }=({\hat{\sigma }}_{\mu }^{\dagger }{)}^{\dagger }\) respectively, with μ th dipole considered. We consider a system of N two-level quantum emitters with |g〉 and |e〉 for the ground and excited states respectively. The quantized bosonic fields \({\hat{a}}_{q}\) satisfy the commutation relations \([{\hat{a}}_{q},{\hat{a}}_{q^{\prime} }^{\dagger }]={\delta }_{q,q^{\prime} }\), and the coupling constant \({g}_{q}\equiv d/\hslash \sqrt{\hslash {\omega }_{q}\mathrm{/(2}{\varepsilon }_{0}V)}({\overrightarrow{\varepsilon }}_{q}\cdot \hat{d})\) involves a dipole moment d with its unit direction \(\hat{d}\), two possible polarizations of the fields \({\overrightarrow{\varepsilon }}_{q}\) with the modes q, and a quantization volume V.

We then consider a confined two-dimensional (2D) reservoir, where 2D lattice array of N two-level atoms are situated. From equation (27), the 2D reservoir has a quantization area A, and we obtain the 2D RDDI of Jμ,ν in polar coordinates,

$$\begin{array}{rcl}{J}_{\mu ,\nu } & = & {\int }_{0}^{\infty }\,\frac{q{\bar{g}}_{q}^{2}A}{{\mathrm{(2}\pi )}^{2}}dq\,{\int }_{0}^{2\pi }\,d\theta \mathrm{[1}-{(\hat{{\bf{q}}}\cdot \hat{{\bf{p}}})}^{2}]{e}^{i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu ,\nu }}[\pi \delta ({\omega }_{q}-{\omega }_{e})\\ & & +\pi \delta ({\omega }_{q}+{\omega }_{e})+i{\mathscr{P}}{({\omega }_{e}-{\omega }_{q})}^{-1}-i{\mathscr{P}}{({\omega }_{q}+{\omega }_{e})}^{-1}],\end{array}$$
(2)

where rμ,ν = rμ − rν, and \(\hat{{\bf{p}}}\) denotes the excitation polarization. We calculate the real part of Jμ,ν first, and \(\hat{{\bf{q}}}\) in general has a polar angle θ to the \(\hat{y}\) on the \(\hat{x}-\hat{y}\) plane. Without loss of generality, we assume rμ,ν along \(\hat{y}\) and \(\hat{{\bf{p}}}\) with a polar angle θ′ to \(\hat{y}\). We obtain

$$\begin{array}{rcl}{\rm{R}}e[{J}_{\mu ,\nu }] & = & {\int }_{0}^{\infty }\,\frac{q{\bar{g}}_{q}^{2}A}{{\mathrm{(2}\pi )}^{2}}dq\,{\int }_{0}^{2\pi }\,d\theta \mathrm{[1}-{(\cos \theta cos\theta ^{\prime} +\sin \theta \sin \theta ^{\prime} )}^{2}]{e}^{i\xi \cos \theta }\pi \delta ({\omega }_{q}-{\omega }_{e}),\\ & = & {{\rm{\Gamma }}}_{2D}\frac{1}{\pi }{\int }_{0}^{2\pi }\,d\theta \mathrm{[1}-{(\cos \theta \cos \theta ^{\prime} +\sin \theta \sin \theta ^{\prime} )}^{2}]{e}^{i\xi \cos \theta },\end{array}$$
(3)

where the intrinsic decay constant for the 2D reservoir is Γ2D ≡ \({k}_{L}|{\partial }_{\omega }q(\omega {)|}_{\omega ={\omega }_{e}}{\bar{g}}_{{k}_{L}}^{2}A\mathrm{/4}\) with the coupling strength \({\bar{g}}_{{k}_{L}}\), the inverse group velocity ∂ωq(ω), and the quantization area A. The dimensionless atomic separation is ξ ≡ kL|rμ − rν| with the near-resonant excitation wave number kL = ωe/c. Integrating out the polar angles, we obtain

$${\rm{Re}}[{J}_{\mu ,\nu }]=2{{\rm{\Gamma }}}_{2D}[{J}_{0}(\xi )-\frac{{J}_{1}(\xi )}{\xi }+{(\hat{{\bf{p}}}\cdot {\hat{{\bf{r}}}}_{\mu ,\nu })}^{2}{J}_{2}(\xi )]\equiv {{\rm{\Gamma }}}_{2D}f(\xi ),$$
(4)

where \({\hat{{\bf{r}}}}_{\mu ,\nu }\) = (rμ − rν)/|rμ − rν|. The above results can be derived from the following integrals,

$${\int }_{0}^{2\pi }\,{e}^{ia\cos \theta }\,d\theta =2\pi {J}_{0}(|a|),$$
(5)
$${\int }_{0}^{2\pi }\,{\cos }^{{\rm{2}}}\,\theta {e}^{ia\cos \theta }d\theta =2\pi (\frac{{J}_{1}(a)}{a}-{J}_{2}(a)),$$
(6)
$${\int }_{0}^{2\pi }\,{\rm{s}}{\rm{i}}{{\rm{n}}}^{2}\theta {e}^{ia\cos \theta }d\theta =2\pi \frac{{J}_{1}(|a|)}{|a|},$$
(7)
$${\int }_{0}^{2\pi }\,\sin \,\theta \,\cos \,\theta {e}^{ia\cos \theta }d\theta =0,$$
(8)

where Jn(a) are the Bessel functions of the first kind.

The Re[Jμ,ν] and Im[Jμ,ν] should satisfy the Kramers-Kronig relation, and therefore we derive Jμ,ν,

$${J}_{\mu ,\nu }={{\rm{\Gamma }}}_{2D}f(\xi )-\frac{i{\mathscr{P}}}{2\pi }{\int }_{0}^{\infty }\,d\omega ({k}_{L}|{\partial }_{\omega }q(\omega )|{\bar{g}}_{q}^{2}A\mathrm{/2)}f(\omega |{{\bf{r}}}_{\mu }-{{\bf{r}}}_{\nu }|/c)(\frac{1}{\omega -{\omega }_{e}}+\frac{1}{\omega +{\omega }_{e}}),$$
(9)
$$=\,{{\rm{\Gamma }}}_{2D}[f(\xi )+ig(\xi )],$$
(10)

where

$$f(\xi )\equiv 2[{J}_{0}(\xi )-\frac{{J}_{1}(\xi )}{\xi }+{(\hat{{\bf{p}}}\cdot {\hat{{\bf{r}}}}_{\mu ,\nu })}^{2}{J}_{2}(\xi )],$$
(11)
$$g(\xi )\equiv 2{Y}_{0}(\xi )-2\frac{{Y}_{1}(\xi )}{\xi }+\mathrm{2(}\hat{{\bf{p}}}\cdot {\hat{{\bf{r}}}}_{\mu ,\nu }{)}^{2}{Y}_{2}(\xi )-\frac{4}{\pi {\xi }^{2}}\mathrm{[1}-\mathrm{2(}\hat{{\bf{p}}}\cdot {\hat{{\bf{r}}}}_{\mu ,\nu }{)}^{2}],$$
(12)

and Yn(ξ) are the Bessel functions of the second kind. The above g(ξ) can be derived by using the following integrals,

$${\mathscr{P}}{\int }_{0}^{\infty }\,da\frac{{J}_{0}(a)}{a\mp b}=-\,\frac{\pi }{2}[{Y}_{0}(b)\pm {H}_{0}(b)],$$
(13)
$${\mathscr{P}}{\int }_{0}^{\infty }\,da\frac{{J}_{1}(a)}{a(a\mp b)}=-\,\frac{2+\pi b[{Y}_{1}(b)\pm {H}_{1}(b)]}{2{b}^{2}},$$
(14)
$${\mathscr{P}}{\int }_{0}^{\infty }\,da(\frac{{J}_{2}(a)}{a-b}+\frac{{J}_{2}(a)}{a+b})=-\,\frac{4}{{b}^{2}}-\pi {Y}_{2}(b),$$
(15)

where Hn(b) is the Struve function.

The dynamical equations for any atomic observables \(Q\equiv \langle \hat{Q}\rangle \) tracing over the fields can be expressed in Lindblad forms (See Methods),

$$\dot{Q}(t)=\sum _{\nu \ne \mu }^{N}\,\sum _{\mu =1}^{N}\,i{\rm{Im}}({J}_{\mu ,\nu })[{\sigma }_{\mu }^{\dagger }{\sigma }_{\nu },Q]+ {\mathcal L} [Q],$$
(16)
$$ {\mathcal L} [Q]=\sum _{\nu =1}^{N}\,\sum _{\mu =1}^{N}\,{\rm{Re}}({J}_{\mu ,\nu })[{\sigma }_{\mu }^{\dagger }Q{\sigma }_{\nu }-\frac{1}{2}({\sigma }_{\mu }^{\dagger }{\sigma }_{\nu }Q+Q{\sigma }_{\mu }^{\dagger }{\sigma }_{\nu }\mathrm{)].}$$
(17)

The coherent and dissipative coupling forms of g(ξ) and f(ξ) respectively denote the frequency shifts and decay rates between any pairs of the atoms. They should satisfy the Kramers-Krönig relation, which is required for causality in a physical response function47.

In Fig. 1, we plot 2D RDDI for two orthogonal light polarizations with \(\hat{{\bf{p}}}\Vert {\hat{{\bf{r}}}}_{\mu ,\nu }\) and \(\hat{{\bf{p}}}\perp {\hat{{\bf{r}}}}_{\mu ,\nu }\) respectively. Both dissipative parts at small ξ are similar and approach unity, which are within Dicke’s superradiant regime. For the coherent parts in the same limit, the leading order of the asymptotics is 2 ln(ξ/2)/π, which is non-analytic at ξ = 0 but diverges much slower than 1/ξ3 in 3D RDDI. As shown in Fig. 1, the frequency shift is still in the order of Γ2D in as short as ξ/(2π) = 0.01, where \(|g(\xi )| \sim 2.1\) and 1.5 respectively in Fig. 1(a,b). This shows a prevailing effect of 2D RDDI on the radiations at such small scale of ξ, in huge contrast to the 3D case where divergent frequency shift forbids any atomic excitations. This promises a short-range and strongly interacting regime in the 2D RDDI, similar to the 1D case as its coherent parts sinξ ≈ 0. We note that g(ξ) goes to −∞ for both parallel and orthogonal dipoles in Fig. 1, in contrast to 3D case where the corresponding collective frequency shift Ωμ,ν(ξ) (with an explicit form in Methods) goes to ∓∞ respectively. This can be attributed to the prefactor of \(\hat{{\bf{p}}}\)⋅\({\hat{{\bf{r}}}}_{\mu ,\nu }\) in Ωμ,ν(ξ) where parallel and orthogonal dipoles change signs of the interaction energy as ξ → 0.

Figure 1
figure 1

Resonant dipole-dipole interactions in a 2D reservoir. The pairwise couplings strengths of collective decay rates [f(ξ), solid] and frequency shifts [g(ξ), dash-dotted], normalized by Γ2D, are shown for light polarizations (a) parallel and (b) perpendicular to atomic separations \({\hat{{\bf{r}}}}_{\mu ,\nu }\).

Full size image

For longer \(\xi \gg 1\), f(ξ) → 1/ξ3/2 (\(\mathrm{1/}\sqrt{\xi }\)) for light polarization parallel (orthogonal) to atomic separations, in contrast to the asymptotic form of RDDI in free space, which is 1/ξ2 (1/ξ). This longer-range dependence of \(\mathrm{1/}\sqrt{\xi }\) is evident in Fig. 1(b), which can be seen as a crossover from 3D to 1D RDDI that eventually lead to infinite-range couplings. This length scaling in this particular polarization configuration can be reinterpreted by ξ−(d−1)/2 where d represents the dimension of reservoir from which RDDI emerge. As a consequence, the f(ξ) in the case of \(\hat{{\bf{p}}}\perp {\hat{{\bf{r}}}}_{\mu ,\nu }\) weakens less rapidly over distances, which can still maintain a significant strength of \(|f(\xi )|/f\mathrm{(0)}\, \sim \,\mathrm{50 \% }\) at \(\xi \gtrsim 10\) or equivalently \({r}_{\mu ,\nu }/{\lambda }_{L}\gtrsim 1.6\). This will make a significant effect on super- and subradiant properties, which are unique from the results in 1D and 3D reservoirs.

Collective super- and subradiant couplings

In the following, we investigate the collective decay constants in a 2D lattice with 2D RDDI, and use ξ hereafter to denote the dimensionless scale of lattice period ds ≡ rμ,μ+1. We consider single photon interacting with an equidistant atomic array, and on absorption the atoms can be excited to the symmetric state,

$$|{\rm{\Psi }}\rangle =\frac{1}{\sqrt{N}}\sum _{\mu =1}^{N}\,{e}^{i{{\bf{k}}}_{L}\cdot {{\bf{r}}}_{\mu }}\,{\sigma }_{\mu }^{\dagger }\mathrm{|0}\rangle ,$$
(18)

where \({e}^{i{{\bf{k}}}_{L}\cdot {{\bf{r}}}_{\mu }}\) is the traveling phase carried by the photon and |0〉 denotes the ground state for all atoms. From the pairwise couplings under the symmetric state, that is \(\langle {\rm{\Psi }}|{J}_{\mu ,\nu }{\sigma }_{\mu }^{\dagger }{\sigma }_{\nu }|{\rm{\Psi }}\rangle \), we obtain the cooperative decay constants and associated frequency shift respectively,

$${{\rm{\Gamma }}}_{N}=\frac{1}{N}\sum _{\nu =1}^{N}\,\sum _{\mu =1}^{N}\,{e}^{-i{{\bf{k}}}_{L}\cdot ({{\bf{r}}}_{\mu }-{{\bf{r}}}_{\nu })}{\rm{Re}}[{J}_{\mu ,\nu }],$$
(19)
$${{\rm{\Delta }}}_{N}=\frac{1}{N}\sum _{\nu \ne \mu }^{N}\,\sum _{\mu =1}^{N}\,{e}^{-i{{\bf{k}}}_{L}\cdot ({{\bf{r}}}_{\mu }-{{\bf{r}}}_{\nu })}\,{\rm{Im}}[{J}_{\mu ,\nu }],$$
(20)

from which the radiation intensity of spontaneously emitted photon can be described by a simple form of exp(−ΓNt + iΔNt). We note that the above sums feature conjugate summands when exchanging μ and ν, and thus two exponentials will combine to a cosine function. In Fig. 2(a), we show the superradiant properties of the symmetric state in a 2D Nx × Nz array with kL along \(\hat{z}\). In Dicke’s limit where \(\xi \ll 1\), we expect of similar ΓN from 2D or 3D RDDI in the same lattice configurations. ΓN saturates quite fast as \({N}_{x}\gg {N}_{z}\), showing an independence of the number of atoms as Nx increases in the direction perpendicular to the light excitation. On the contrary, two contrasting dependences of Nz can be located at \({N}_{x}\ll 10\) and \({N}_{x}\gtrsim 10\), which are \( \sim {N}_{z}^{0.65}\) and \({N}_{z}^{0.97}\) respectively for Nx = 2 and 30. This shows a suppressed scaling in a needle-like 2D lattice compared to the square structure. Similar distinguishing features are also present in ΔN, where the needle-like structure allows significant red shifts, whereas for \({N}_{z} < {N}_{x}\), blue shifts emerge instead. In the region of Nz < Nx, we have a relatively broad and less varying dependence of lattice structures.

Figure 2
figure 2

Super- and subradiant decay constants and frequency shifts. (a) Superradiant decay constants ΓN and frequency shifts ΔN are shown in the upper and lower panels respectively at ξ = 1 with \(\hat{{\bf{p}}}\Vert \hat{x}\), where the lattice period is ds = λL/(2π). Specific numeric values are labeled on every other contour lines for clarity. (b) Some selective ξ's show subradiant decay behaviors as a dependence of Nx for Nz = 1 (solid and dash-dotted), 2 (dotted), and 3 (dashed), in the same configuration of \(\hat{{\bf{p}}}\). ξ = 5, 6, 9, 40 correspond to ds/λL = 0.8, 0.95, 1.43, 6.37 respectively.

Full size image

More interesting decay behavior of 2D RDDI results from the oscillatory negative couplings in the case of \(\hat{{\bf{p}}}\perp {\hat{{\bf{r}}}}_{\mu ,\nu }\) in Fig. 1(b). As shown in Fig. 2(b), the subradiant decay can be supported at some selective ξ's in an optically-thin lattice structure. This even sustains in longer distance, for example of ξ = 40 in the plot. If we put this 2D lattice mediating 3D RDDI in free space, the subradiance under the symmetric states becomes less significant (\( \sim \,\mathrm{15 \% }\) more of the ΓN at \({N}_{x}\, \sim \,40\)), and thus 2D RDDI comparing the 3D case specifically show a notable long-range effect, resembling the infinite-range sinusoidal forms of 1D RDDI.

Phase-imprinted subradiant states

Next, we further study the subradiance from 2D RDDI, which can be enabled by imprinting linearly increasing phases on the atomic arrays20,21 or via a side excitation with a π phase shift in the sub-ensembles18. This spatially varying phase can be imprinted by applying pulsed gradient magnetic or electric fields, or directly using light which carries orbital angular momentum in atomic ring structures26,27. We first construct a complete Hilbert space of single excitation, which reads

$$|{{\rm{\Phi }}}_{m}\rangle =\frac{1}{\sqrt{N}}\sum _{\mu =1}^{N}\,{e}^{i{{\bf{k}}}_{L}\cdot {{\bf{r}}}_{\mu }}{e}^{i2m\pi (\mu -\mathrm{1)/}N}\,{\sigma }_{\mu }^{\dagger }\mathrm{|0}\rangle ,$$
(21)

where m ∈ [1, N]. The above set should be orthonormal, where the inner products of them is \(\langle {{\rm{\Phi }}}_{m}|{{\rm{\Phi }}}_{n}\rangle ={N}^{-1}\,{\sum }_{\mu \mathrm{=1}}^{N}\) \({e}^{\frac{i2\pi (\mu -\mathrm{1)(}m-n)}{N}}\) = δm,n, satisfied by De Moivre formula. The Hilbert space of equation (21) includes both super- and subradiant states, which has also been applied in forward- and backward-propagating eigenstates to reveal the emergent universal dynamics in a 1D nanophotonic system34. Note that in general there are infinite ways to construct singly-excited Hilbert space, and therefore equation (21) is not unique to the setting of N atoms interacting with single photon. The states of equation (21) can be prepared collectively where all atoms are excited uniformly, and allow studies on super- and sub-radiance systematically by varying the imprinted phases. Though these states can be controlled dynamically, their fidelities may suffer from an inefficient phase-imprinting protocol using pulsed lasers or limitation of large gradient magnetic fields20. Nevertheless, the phase imprinting construction allows a controllable way to manipulate these orthonormal states collectively.

Since 2D RDDI involve a long-range functional form, it is not possible to write down the analytical eigenstates in general. Therefore, we numerically derive the eigenbases, and the time evolutions of |Φm(t)〉 can be obtained by solving the Schrödinger equations, ∂|Φ(t)〉/∂t = −J|Φ(t)〉, where the matrix elements of J consists of the couplings \({J}_{\mu ,\nu }^{\ast }\), and \(|{\rm{\Phi }}(t)\rangle ={\sum }_{\mu }{b}_{\mu }(t){\sigma }_{\mu }^{\dagger }\mathrm{|0}\rangle \). By diagonalizing J, we obtain \(\overrightarrow{b}(t)\) = \(S{e}^{\overrightarrow{\lambda }t}{S}^{-1}\overrightarrow{b}(t\mathrm{=0)}\), where S and \(\overrightarrow{\lambda }\) are the eigenbases and eigenvalues respectively. For some initially prepared state |Φ(t = 0)〉 = |Φm〉, we obtain its time evolution as \({A}_{m}(t)\equiv {\sum }_{l}{w}_{l}(m){e}^{{\lambda }_{l}t}\), where wl(m) = \(({h}^{\dagger }S)\cdot ({S}^{-1}h)\) leads to the weightings |wl(m)|2 of |Φm〉 on respective l th eigenmodes, with a column vector h consisting of the imprinted phases \({N}^{-\mathrm{1/2}}{e}^{i{\bf{k}}\cdot {{\bf{r}}}_{\mu }+i2m\pi (\mu -\mathrm{1)/}N}\).

In Fig. 3(a), we show the distributions of the eigenmodes in an ascending order in a 2D square lattice. When \(\xi \lesssim 5\) or the mutual distance is less than the resonant wavelength, the 2D system allows significant super- and subradiant eigen-decay constants, as expected and similar to the results from 3D RDDI in a strongly interacting regime. By contrast, as ξ extends further, 2D RDDI still permit the lowest decay rate below 10−2 Γ2D in Fig. 3(a), indicating of long-range atom-atom correlations. As a comparison, in the same 2D lattice configuration but in a 3D reservoir, the eigenmodes show a level dependence and have reached the noninteracting regime. In Fig. 3(b), as an example, we further show the radiation intensity |Am(t)|2 of two subradiant states assuming they are initially created. Increasing m means larger gradient fields required to prepare these states. For the selective state of m = 5, we see quite a slow subradiant decay with two beating frequencies as time evolves. This originates from three dominating eigenmodes as shown in the weightings of the inset (less obvious for \(|{w}_{10}(m)| \sim \,0.15\)), where two most significant modes of l = 11 and 18 occupy relatively small eigen-decay constants of \( \sim 2\times {10}^{-4}\) and \( \sim 8\times {10}^{-4}{{\rm{\Gamma }}}_{2D}\) respectively. The beating frequencies can be determined by the differences of Im(λl) in the respective modes. Other example of m = 7 state in the lower plot of Fig. 3(b) also shows the subradiant decay behavior with four significant weightings on the eigenmodes instead. The lifetime of the oscillatory waveform can be approximately decided by the most significant mode of l = 39 which has an eigen-decay rate \( \sim 8\times {10}^{-3}{{\rm{\Gamma }}}_{2D}\). In addition to the fast oscillating ripples in the radiation pattern in both Fig. 3(a,b), a slowly-varying envelope also appears and extends to long time scales, implying the dominance of the subradiant modes.

Figure 3
figure 3

Eigen-decay rates and time evolutions of subradiant states for 10 × 10 array. (a) The eigen-decay constants can be obtained by the real parts of the eigenvalues λl at ξ = 1 (solid), 5 (dash-dotted), 10 (dotted), where l denotes the l th eigenvalues from the eigen-analysis. ξ = 1, 5, 10 correspond to ds/λL = 0.16, 0.8, 1.6 respectively. As a comparison, we show the results for the same 2D lattice configuration from 3D RDDI at ξ = 10 (dashed), and the horizontal line guides the eye for a natural decay constant. (b) Time dynamics of some selective subradiant states of m = 5 and 7 at ξ = 1, with most notable state weightings |wl(m)|2 on the eigenmodes in the respective insets. m denotes the strength of imprinting phase gradient. Clear beatings in the radiation can be seen in both plots due to finite Im(λl).

Full size image

Finally, we study a striped 2D lattice structure for two orthogonal light excitations in Fig. 4. For both super- and subradiant states in the example of 2D lattice with \({N}_{z}\gg {N}_{x}\), the decay behavior can be approximately separated into two time scales, an early fast drop and late subradiant decay, which also manifests in a dilute but optically thick 3D cloud17. In Fig. 4(a), the superradiant state of m = 0 becomes exactly the symmetric state of equation (18), where ΓN governs the decay behavior in the beginning when |Am(t)|2 \(\gtrsim \) 0.01. The later oscillatory subradiance indicates of multiple though less occupied subradiant modes. Comparing the lifetime determined when its initial probability drops to e−1 in the early stage, an optically-thick striped lattice in the case of \({{\bf{k}}}_{L}\Vert \hat{z}\) shows an enhanced decay rate by only a factor of \( \sim 4\) over the case of \({{\bf{k}}}_{L}\Vert \hat{x}\). On the other hand for the subradiant states in Fig. 4(b) with a finite phase imprinting, the contrasting reduction factor of the decay rates becomes \( \sim 100\) in the optically-thick configuration. This magnifying factor in the subradiant time scale suggests a potential photon routing relying on the 2D lattice mediating 2D RDDI, where light going through an optically-thick direction delays and almost stops within the time \( \sim 100{{\rm{\Gamma }}}_{2D}^{-1}\). Furthermore, potential chiral implementations using the phase-imprinted many-body states can be feasible in various atomic systems, for example cavity-optomechanical circuits48, 2D coupled ring resonators44, or superconducting qubits and quantum dots in the photonic waveguides49 under an effectively emulated 2D reservoir.

Figure 4
figure 4

Contrasting time evolutions at ξ = 5 (ds/λL = 0.8) for a Nx × Nz = 4 × 20 array. (a) Superradiant state dynamics of m = 0 with kL along \(\hat{z}\) (solid) and \(\hat{x}\) (dotted) respectively, comparing the natural decay of \({e}^{-2{\Gamma }_{2D}t}\) (dashed). (b) Subradiant state evolutions of m = 10 with kL along \(\hat{z}\) (solid) and \(\hat{x}\) (dotted) respectively. The inset shows the long time behavior for the subradiant state in the case of kL along \(\hat{z}\), which shows an early drop with an oscillatory subradiance afterward, similar to (a).

Full size image

Conclusion

In conclusion, we have derived the explicit form of the RDDI from a confined two-dimensional reservoir. We demonstrate distinctive characteristics of 2D RDDI, which allows subradiance under a singly-excited symmetric state more significantly than the 3D case. This indicates long-range atom-atom correlations which are different from the induced RDDI in either 1D or 3D reservoirs. By imprinting spatially dependent phases on the 2D atomic arrays, we propose to prepare single-excitation subradiant states in a potentially deterministic and controllable way. Our results put forward potential applications in manipulating quantum information and preparations of many-body subradiant states in a 2D reservoir.

Methods

General formalism for resonant dipole-dipole interaction in a three-dimensional reservoir

Here we review the general formalism of resonant dipole-dipole interaction (RDDI)4,5 in a free space of three-dimensional (3D) reservoir. The RDDI originates from the common quantized light fields rescattering multiple times in the dissipation process. This collective dipole-dipole interaction in an ensemble of two-level quantum emitters is responsible for cooperative spontaneous emissions, so-called superradiance1,2 and subradiance, and collective frequency shift50,51. Only recently that significantly small collective frequency shift can be observed in some versatile atomic systems, including the embedded atoms in the planar cavity52, a vapor cell53, an ionic system54, and cold atoms7.

The spontaneous decay behavior in a system of N two-level quantum emitters, with |g〉 and |e〉 for the ground and excited states respectively, can be described by a 3D reservoir of quantized bosonic light fields interacting with the medium. With a dipole approximation, the Hamiltonian reads5,

$$H=\sum _{\mu =1}^{N}\,\hslash {\omega }_{e}{\hat{\sigma }}_{\mu }^{\dagger }{\hat{\sigma }}_{\mu }-\sum _{\mu =1}^{N}\,\sum _{q}\,{g}_{q}({e}^{i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }-i{\omega }_{q}t}{\hat{a}}_{q}+{e}^{-i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }+i{\omega }_{q}t}{\hat{a}}_{q}^{\dagger })({\hat{\sigma }}_{\mu }+{\hat{\sigma }}_{\mu }^{\dagger }),$$
(22)

where the atomic raising operator is \({\hat{\sigma }}_{\mu }^{\dagger }\equiv |e{\rangle }_{\mu }\langle g|\) with \({\hat{\sigma }}_{\mu }=({\hat{\sigma }}_{\mu }^{\dagger }{)}^{\dagger }\), and quantized fields \({\hat{a}}_{q}\) should satisfy the bosonic commutation relations \([{\hat{a}}_{q},{\hat{a}}_{q^{\prime} }^{\dagger }]={\delta }_{q,q^{\prime} }\). The coupling constant \({g}_{q}\equiv d/\hslash \sqrt{\hslash {\omega }_{q}\mathrm{/(2}{\varepsilon }_{0}V)}({\overrightarrow{\varepsilon }}_{q}\cdot \hat{d})\) involves a dipole moment d with its unit direction \(\hat{d}\), two possible polarizations of the fields \({\overrightarrow{\varepsilon }}_{q}\) with the modes q, and a quantization volume V. The above Hamiltonian involves the non-rotating wave terms which are necessary for a complete description of the frequency shift (dispersion) of the RDDI in the dissipation. Therefore, the dispersion and absorption of RDDI should satisfy the Kramers-Kronig relation.

Following the derivations in ref.5, we continue to formulate a Heisenberg equation for an atomic operator \(\hat{Q}\), that is \(d\hat{Q}/dt=i[H,\hat{Q}]\) (let \(\hslash =1\)). We obtain

$$\frac{d\hat{Q}}{dt}=i{\omega }_{e}\sum _{\mu }\,[{\hat{\sigma }}_{\mu }^{\dagger }{\hat{\sigma }}_{\mu },\hat{Q}]-i\sum _{\mu }\,\sum _{q}\,{g}_{q}\{{e}^{i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }}[{\hat{\sigma }}_{\mu }+{\hat{\sigma }}_{\mu }^{\dagger },\hat{Q}]{\hat{a}}_{q}(t)-{e}^{-i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }}{\hat{a}}_{q}^{\dagger }(t)[\hat{Q},{\hat{\sigma }}_{\mu }+{\hat{\sigma }}_{\mu }^{\dagger }\mathrm{]\}.}$$
(23)

In the above, \({\hat{a}}_{q}\) can be further substituted by solving \(d{\hat{a}}_{q}/dt=i[H,{\hat{a}}_{q}]\), which is

$${\hat{a}}_{q}(t)={\hat{a}}_{q}\mathrm{(0)}{e}^{-i{\omega }_{q}t}+i\sum _{\mu }\,{g}_{q}{e}^{-i{{\bf{k}}}_{q}\cdot {{\bf{r}}}_{\mu }}{\int }_{0}^{t}\,dt^{\prime} [{\hat{\sigma }}_{\mu }(t^{\prime} )+{\hat{\sigma }}_{\mu }^{\dagger }(t^{\prime} )]{e}^{-i{\omega }_{q}(t-t^{\prime} )}.$$
(24)

With the Born-Markov approximation of \({\omega }_{e}t\gg 1\) and \(t\gg {({r}_{\mu \nu })}_{max}/c\) (rμν ≡ |rμ − rν|), we obtain the dynamical equation of \(Q\equiv {\langle \hat{Q}\rangle }_{0}\) in Lindblad forms by considering the vacuum initial bosonic fields 〈〉0,

$$\dot{Q}(t)=\sum _{\mu \ne \nu }\,i{{\rm{\Omega }}}_{\mu ,\nu }[{\sigma }_{\mu }^{\dagger }{\sigma }_{\nu },Q]+ {\mathcal L} (Q),$$
(25)
$$ {\mathcal L} (Q)=\sum _{\mu ,\nu }\,{\gamma }_{\mu ,\nu }[{\sigma }_{\mu }^{\dagger }Q{\sigma }_{\nu }-\frac{1}{2}({\sigma }_{\mu }^{\dagger }{\sigma }_{\nu }Q+Q{\sigma }_{\mu }^{\dagger }{\sigma }_{\nu })].$$
(26)

The Ωμ,ν and γμ,ν describe the collective frequency shifts and decay rates respectively. These represent the coherent and dissipative parts of the pairwise couplings, Jμ,ν ≡ (γμ,ν + i2Ωμ,ν)/2, which are defined as

$$\begin{array}{rcl}{J}_{\mu ,\nu } & = & \sum _{q}\,|{g}_{q}{|}^{2}{\int }_{0}^{\infty }\,dt^{\prime} {e}^{i{{\bf{k}}}_{q}\cdot ({{\bf{r}}}_{\mu }-{{\bf{r}}}_{\nu })}[{e}^{i({\omega }_{e}-{\omega }_{q})t^{\prime} }+{e}^{-i({\omega }_{e}+{\omega }_{q})t^{\prime} }],\\ & = & \sum _{q}\,|{g}_{q}{|}^{2}{\int }_{0}^{\infty }\,dt^{\prime} {e}^{i{{\bf{k}}}_{q}\cdot ({{\bf{r}}}_{\mu }-{{\bf{r}}}_{\nu })}[\pi \delta ({\omega }_{q}-{\omega }_{e})+\pi \delta ({\omega }_{q}+{\omega }_{e})+i{\mathscr{P}}{({\omega }_{e}-{\omega }_{q})}^{-1}-i{\mathscr{P}}{({\omega }_{q}+{\omega }_{e})}^{-1}],\end{array}$$
(27)

where \({\mathscr{P}}\) is the principal value of the integral.

For a 3D reservoir, we consider continuous limits of modes \({\sum }_{q}\to {\sum }_{{\overrightarrow{\varepsilon }}_{q}}{\int }_{-\infty }^{\infty }\frac{V}{{\mathrm{(2}\pi )}^{3}}{d}^{3}q\) with two possible field polarizations \({\overrightarrow{\varepsilon }}_{q}\). In spherical coordinates, we show the main results of Jμ,ν in free space5,

$$\begin{array}{rcl}{\gamma }_{\mu ,\nu }(\xi ) & \equiv & \oint \,d{{\rm{\Omega }}}_{q}\mathrm{[1}-{(\hat{{\bf{q}}}\cdot \hat{{\bf{p}}})}^{2}]{\int }_{0}^{\infty }\,dq{q}^{2}{\bar{g}}_{q}^{2}\frac{V}{{\mathrm{(2}\pi )}^{3}}[\pi \delta ({\omega }_{q}-{\omega }_{e})+\pi \delta ({\omega }_{q}+{\omega }_{e})],\\ & = & \frac{3{\rm{\Gamma }}}{2}\{[1-{(\hat{{\bf{p}}}\cdot {\hat{r}}_{\mu \nu })}^{2}]\frac{\sin \,\xi }{\xi }+[1-\mathrm{3(}\hat{{\bf{p}}}\cdot {\hat{r}}_{\mu \nu }{)}^{2}](\frac{\cos \,\xi }{{\xi }^{2}}-\frac{\sin \,\xi }{{\xi }^{3}})\},\end{array}$$
(28)
$$\begin{array}{rcl}{{\rm{\Omega }}}_{\mu ,\nu }(\xi ) & \equiv & -\oint \,d{{\rm{\Omega }}}_{q}\mathrm{[1}-{(\hat{{\bf{q}}}\cdot \hat{{\bf{p}}})}^{2}]{\int }_{0}^{\infty }\,dq{q}^{2}{\bar{g}}_{q}^{2}\frac{V}{{\mathrm{(2}\pi )}^{3}}[i{\mathscr{P}}{({\omega }_{q}-{\omega }_{e})}^{-1}+i{\mathscr{P}}{({\omega }_{q}+{\omega }_{e})}^{-1}],\\ & = & \frac{3{\rm{\Gamma }}}{4}\{-\mathrm{[1}-{(\hat{{\bf{p}}}\cdot {\hat{r}}_{\mu \nu })}^{2}]\frac{\cos \,\xi }{\xi }+\mathrm{[1}-\mathrm{3(}\hat{{\bf{p}}}\cdot {\hat{r}}_{\mu \nu }{)}^{2}](\frac{\sin \,\xi }{{\xi }^{2}}+\frac{\cos \,\xi }{{\xi }^{3}})\},\end{array}$$
(29)

where dΩq denotes an integration of a solid angle, \({\bar{g}}_{q}^{2}\) ≡ \({(d/\hslash )}^{2}[\hslash {\omega }_{q}\mathrm{/(2}{\varepsilon }_{0}V)]\), \(\hat{{\bf{p}}}\) parallels the excitation field polarization, the natural decay constant \({\rm{\Gamma }}={d}^{2}{\omega }_{e}^{3}\mathrm{/(3}\pi \hslash {\varepsilon }_{0}{c}^{3})\), and dimensionless ξ ≡ kL|rμ − rν| with kL = ωe/c. As ξ → 0, Dicke’s regime is reached where γμ,ν → Γ, while Ωμ,ν goes to infinity. This divergence shows the inapplicability of quantum optical treatment in RDDI or in other words, it simply forbids any atomic excitations by external fields.