Subradiant Decay in 2D and 3D Atomic Arrays

1. Introduction

Dicke superradiance corresponds to the enhancement of the radiation rate by N emitters due to coherence effects [1,2], as opposed to subradiance for which it is inhibited [3,4]. Recently, the subradiant suppression of radiative decay has spurred strong interest due to its potential applications for quantum memories [5,6], excitation transfer [7,8] and topological photonics [9]. Many-atom subradiance has been observed almost ten years ago in disordered clouds of cold atoms, with long lifetimes evidencing the phenomenon [10].

More recently, ordered sub-wavelength arrays of atoms have emerged as a promising configuration to achieve subradiance, with distinct advantages over disordered ensembles for applications such as photon storage or photonic gates [11,12,13,14,15,16]. In these systems, the emitters are regularly arranged with a spatial period below the atomic transition wavelength, so the dipolar interactions lead to strong cooperative optical properties.

In this context, infinite lattices offer interesting insights on the scattering properties of these arrays. The single-excitation eigenmodes of the scattering problem are translationally invariant and obey Bloch’s theorem. This allows one to define the quasi-momentum $\mathbf{k}$ and calculate explicitly the collective frequency shift and decay rate of each of these modes [15]. Differently, in finite systems these collective features can only be accessed by numerical diagonalization, leading to severe limitations on the system size which can be simulated. Furthermore, most studies on subradiance in ordered systems have focused on one-dimensional (1D) atomic chains [6,7,8,17,18,19,20,21] and two-dimensional (2D) arrays [9,13,14,15,22], with investigations on the three-dimensional configurations lacking [23]. We note that while the band structure of 2D and 3D arrays have been studied through the density of states [24,25,26], the connection between this phenomenon and subradiance has not been yet clarified.

In this work, we investigate the cooperative decay rates of 2D and 3D atomic array, using a technique where the discrete sums over the emitters is turned into an angular integral. It provides the cooperative rate for the quasi-modes in infinite arrays [27]. But it is also particular convenient to study the rates of large size systems, for which the technique provides a good approximation, and numerical diagonalization are out of reach. Used successfully for 1D systems [27], we show the collective decay rates of the associated generalized Dicke state (also called generalized Bloch state [28]) can also be derived explicitly for 2D and 3D arrays. The work thus opens new perspectives for the study of cooperative scattering in large, regular atomic systems.

2. Microscopic Modeling of Atomic Arrays

2.1. Coupled Dipole Model

Let us consider N atoms with transition frequency ${\omega }_0=c{k}_0$ between the ground state $|g_j\rangle$ and the excited state $|e_j\rangle$ ($j=1,\dots ,N$), linewidth ${\Gamma }_0=k_0^3{d}_{ge}^2/(3\pi {ϵ}_0\hslash )$, where ℏ is the Planck constant, ${ϵ}_0$ the vacuum permittivity and $d_{ge}$ the electric-dipole transition matrix element. The atoms are at fixed positions ${\mathbf{r}}_j$ and have transition dipole moments ${\mathbf{d}}_j$, with ${\mathbf{d}}_j=d_{ge}{\widehat{\mathbf{d}}}_j$. In the Born-Markov approximation, one can trace out the quantized light fields and obtain the field-mediated dipole-dipole couplings between the atoms described by an effective Hamiltonian [29]:

$$H=-\hslash \sum _{j,m=1}^N{\widehat{\mathbf{d}}}_j^{*}·\mathbf{G}({\mathbf{r}}_j-{\mathbf{r}}_m)·{\widehat{\mathbf{d}}}_m{\sigma }_j^{†}{\sigma }_m,$$

(1)

where ${\sigma }_j=|g_j\rangle \langle e_j|$ is the lowering operator for atom j. The dyadic Green’s tensor $\mathbf{G}$ for the electromagnetic field in the 3D free space in vacuum is given by

$$\mathbf{G}\left(\mathbf{r}\right)={\displaystyle \frac{3{\Gamma }_0}{2}}{\displaystyle \frac{e^{i{k}_{0}r}}{k_{0}r}}\left\{\mathbb{I}-\widehat{\mathbf{n}}\widehat{\mathbf{n}}+\left(\mathbb{I}-3\widehat{\mathbf{n}}\widehat{\mathbf{n}}\right)\left({\displaystyle \frac{i}{k_{0}r}}-{\displaystyle \frac{1}{k_0^2{r}^2}}\right)\right\},$$

(2)

where $\mathbb{I}$ is the identity tensor $3\times 3$, $r=\left|\mathbf{r}\right|$ and $\widehat{\mathbf{n}}=\mathbf{r}/r$. The dyadic Green’s tensor (2) can be obtained from the scalar Green’s function $G\left(r\right)=exp\left(i{k}_{0}r\right)/\left(k_{0}r\right)$ [30]

$$\mathbf{G}\left(\mathbf{r}\right)={\displaystyle \frac{3{\Gamma }_0}{2}}\left[\mathbb{I}+{\displaystyle \frac{1}{k_0^2}}\nabla \nabla \right]G\left(r\right).$$

(3)

If we assume that all transition dipoles point in the same direction, ${\widehat{\mathbf{d}}}_j=\widehat{\mathbf{d}}$, such as in presence of a strong external field, the interaction term between atoms j and m read

$$G_{jm}={\widehat{\mathbf{d}}}_j^{*}·\mathbf{G}\left({\mathbf{r}}_{jm}\right)·{\widehat{\mathbf{d}}}_m={\displaystyle \frac{3\Gamma }{2}}\left[1+{\displaystyle \frac{1}{k_0^2}}{(\widehat{\mathbf{d}}·{\nabla }_{{\mathbf{r}}_{jm}})}^2\right]{\displaystyle \frac{e^{i{k}_0{r}_{jm}}}{k_0{r}_{jm}}},$$

(4)

with ${\mathbf{r}}_{jm}={\mathbf{r}}_j-{\mathbf{r}}_m$ and $r_{jm}=\left|{\mathbf{r}}_{jm}\right|$. Note that $G_{jm}$ has both a real and imaginary part, which describe the photon exchange between the atoms and their collective emission out of the system, respectively. Throughout this work we focus our attention on the imaginary part of $G_{jm}$, defining

$$\begin{array}{ccc}\hfill {\Gamma }_{jm}& =& \mathrm{Im}\left[G_{jm}\right]={\displaystyle \frac{3{\Gamma }_0}{2}}\left[1+{\displaystyle \frac{1}{k_0^2}}{(\widehat{\mathbf{d}}·{\nabla }_{{\mathbf{r}}_{jm}})}^2\right]{\displaystyle \frac{sin\left(k_0{r}_{jm}\right)}{k_0{r}_{jm}}}.\hfill \end{array}$$

(5)

Introducing the angular average

$$\begin{array}{c}\hfill {〈f({\theta }_0,{\varphi }_0)〉}_{\Omega }={\displaystyle \frac{1}{4\pi }}{\int }_0^{2\pi }d{\varphi }_0{\int }_0^{\pi }sin{\theta }_{0}f({\theta }_0,{\varphi }_0)d{\theta }_0,\end{array}$$

and using the property ${\langle e^{-ixcos\theta }\rangle }_{\Omega }=(sinx)/x$, we can rewrite ${\Gamma }_{jm}$ as an angular average of the photon exchange term between the two atoms

$$\begin{array}{ccc}\hfill {\Gamma }_{jm}& =& {\displaystyle \frac{3{\Gamma }_0}{2}}{〈\left[1-{(\widehat{\mathbf{d}}·{\widehat{\mathbf{k}}}_0)}^2\right]e^{-i{\mathbf{k}}_0·({\mathbf{r}}_j-{\mathbf{r}}_m)}〉}_{\Omega },\hfill \end{array}$$

(6)

where ${\mathbf{k}}_0=k_0(sin{\theta }_{0}cos{\varphi }_0,sin{\theta }_{0}sin{\varphi }_0,cos{\theta }_0)$ and ${\widehat{\mathbf{k}}}_0={\mathbf{k}}_0/k_0$. The remarkable property of Eq. (6) is that the contribution of the two atoms can now be factorized, which will allow us to calculate explicitly the collective decay rates, as we will later see.

2.2. Generalized Dicke State

We hereafter consider regular atomic arrays with lattice step d, with atoms at positions ${\mathbf{r}}_{j_x,j_y,j_z}=d[(j_x-1){\widehat{\mathbf{e}}}_x+(j_y-1){\widehat{\mathbf{e}}}_y+(j_z-1){\widehat{\mathbf{e}}}_z]$ with $j_x=1,\dots N_x$, $j_y=1,\dots ,N_y$ and $j_z=1,\dots ,N_z$, so $N=N_x{N}_y{N}_z$ is the total number of atoms. We then introduce the generalized Dicke states [27], or generalized Bloch states [28]:

$$|\mathbf{k}\rangle ={\displaystyle \frac{1}{\sqrt{N}}}\sum _{j=1}^N{e}^{i\mathbf{k}·{\mathbf{r}}_j}|j\rangle .$$

(7)

The single-excitation states are $|j\rangle =|g_1,\dots ,e_j,\dots ,g_N\rangle$ and $\mathbf{k}\in \mathcal{S}$, where $\mathcal{S}$ is the volume of first Brillouin zone in the reciprocal lattice, i.e., $|k_x|,|k_y|,|k_z|\le \pi /d$. It includes the symmetric Dicke state [1] for $\mathbf{k}=0$, as well as the timed-Dicke state introduced in Ref. [31,32] for a driving with wavevector $\mathbf{k}={\mathbf{k}}_0$. The family $\left\{\right|\mathbf{k}\rangle \}$ satisfies the completeness relation

$${\left({\displaystyle \frac{d}{2\pi }}\right)}^3{\int }_{\mathcal{S}}{d}^{3}k|\mathbf{k}\rangle \langle \mathbf{k}|=1.$$

(8)

For a finite lattice these states are not orthogonal since

$$\langle {\mathbf{k}}^{\prime }|\mathbf{k}\rangle ={\displaystyle \frac{sin[(k_x-k_x^{\prime })d{N}_x/2]}{sin[(k_x-k_x^{\prime })d/2]}}{\displaystyle \frac{sin[(k_y-k_y^{\prime })d{N}_y/2]}{sin[(k_y-k_y^{\prime })d/2]}}{\displaystyle \frac{sin[(k_z-k_z^{\prime })d{N}_z/2]}{sin[(k_z-k_z^{\prime })d/2]}}{e}^{i(\mathbf{k}-{\mathbf{k}}^{\prime })·{\mathbf{r}}_{N_x,N_y,N_z}}.$$

(9)

nevertheless, in the limit of an infinite lattice the asymptotic property $\langle {\mathbf{k}}^{\prime }|\mathbf{k}\rangle \to {\delta }^{\left(3\right)}(\mathbf{k}-{\mathbf{k}}^{\prime })$ turns the set $\left\{\right|\mathbf{k}\rangle \}$ into an orthonormal basis for the single-excitation manifold.

2.3. Collective Decay Rate

The collective decay rate of mode $\mathbf{k}$ is defined as [27]

$$\Gamma \left(\mathbf{k}\right)=-{\displaystyle \frac{1}{\hslash }}\mathrm{Im}\langle \mathbf{k}\left|H\right|\mathbf{k}\rangle ={\displaystyle \frac{1}{N}}\sum _{j=1}^N\sum _{m=1}^N{\Gamma }_{jm}{e}^{i\mathbf{k}·{\mathbf{r}}_{jm}}.$$

(10)

Eq. (10) corresponds to the Fourier transform of ${\Gamma }_{jm}$ only if the lattice is infinite, when modes (7) are the Bloch states and set $\left\{\mathbf{k}\right\}$ is discrete. Nevertheless, treating $\mathbf{k}$ as a continuous variable in Eq. (10) is a good approximation for the mode when N is sufficiently large [15]. We hereafter refer to $\Gamma \left(\mathbf{k}\right)$ as the “continuous spectrum” of the decay rate. We highlight that we are not referring to the discrete spectrum of the eigenvalues of the system, described by the non-Hermitian Hamiltonian (1), but rather to the imaginary part of the expectation value of the effective Hamiltonian on the collective state of Eq. (7), describing the interference among the (single-photon) fields emitted by N atoms. In particular, in the infinite N limit, $\Gamma \left(\mathbf{k}\right)$ is precisely the decay rate of mode $\mathbf{k}$, yet in the finite-N case, the exact rate must be obtained by diagonalizing the full matrix (1).

Inserting the cross decay term (6) into (10) leads to

$$\begin{array}{ccc}\hfill \Gamma \left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\Gamma }_0}{2N}}\sum _{j=1}^N\sum _{m=1}^N{〈\left[1-{(\widehat{\mathbf{d}}·{\widehat{\mathbf{k}}}_0)}^2\right]e^{i(\mathbf{k}-{\mathbf{k}}_0)·{\mathbf{r}}_{jm}}〉}_{\Omega }\hfill \\ & =& {\displaystyle \frac{3{\Gamma }_0}{2N}}{〈\left[1-{(\widehat{\mathbf{d}}·{\widehat{\mathbf{k}}}_0)}^2\right]{\left|F\left(\mathbf{k}\right)\right|}^2〉}_{\Omega },\hfill \end{array}$$

(11)

where we have introduced the structure factor

$$F\left(\mathbf{k}\right)=\sum _{j=1}^N{e}^{i(\mathbf{k}-{\mathbf{k}}_0)·{\mathbf{r}}_j}.$$

(12)

Using this explicit equation, let us now investigate the decay rates for a 2D square lattice and a 3D cubic lattice.

3. Two-Dimensional Square Lattice

3.1. General Expression of the Decay Rate

For a 2D atom array in the plane $(x,y)$, with ${\mathbf{r}}_j=d[(j_x-1){\widehat{\mathbf{e}}}_x+(j_y-1){\widehat{\mathbf{e}}}_y]$ with $j_x=1,\dots N_x$ and $j_y=1,\dots ,N_y$, the rate (11) reads

$$\begin{array}{ccc}\hfill \Gamma (k_x,k_y)& =& {\displaystyle \frac{3{\Gamma }_0}{2N_x{N}_y}}{〈\left[1-{(\widehat{\mathbf{d}}·{\widehat{\mathbf{k}}}_0)}^2\right]{\left|F(k_x,k_y)\right|}^2〉}_{\Omega },\hfill \end{array}$$

(13)

where the structure factor now factorizes for each direction:

$$F(k_x,k_y)=\left(\sum _{j_x=1}^{N_x}{e}^{i(k_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)d(j_x-1)}\right)\left(\sum _{j_y=1}^{N_y}{e}^{i(k_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)d(j_y-1)}\right).$$

(14)

Summing over $j_x$ and $j_y$, we obtain

$$|F(k_x,k_y){|}^2={\displaystyle \frac{{sin}^2[(k_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)d{N}_x/2]}{{sin}^2[(k_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)d/2]}}\times {\displaystyle \frac{{sin}^2[(k_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)d{N}_y/2]}{{sin}^2[(k_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)d/2]}}.$$

(15)

The property

$${\displaystyle \frac{{sin}^2\left(Nt\right)}{{sin}^{2}t}}\approx N^2\sum _{m=-\infty }^{+\infty }{\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]$$

(16)

can then be used, along with the assumption of large $N_x$ and $N_y$, to write the decay rate as

$$\begin{array}{ccc}\hfill \Gamma (k_x,k_y)& =& {\displaystyle \frac{3{\Gamma }_0}{8\pi }}{N}_x{N}_y\sum _{(m_x,m_y)\in {\mathbb{Z}}^2}{\int }_0^{2\pi }d{\varphi }_0{\int }_0^{\pi }d{\theta }_{0}sin{\theta }_0\hfill \\ & \times & \left[1-{({\widehat{d}}_{x}sin{\theta }_{0}cos{\varphi }_0+{\widehat{d}}_{y}sin{\theta }_{0}sin{\varphi }_0+{\widehat{d}}_{z}cos{\theta }_0)}^2\right]\hfill \\ & \times & {\mathrm{sinc}}^2\left[\left(k_x-{\displaystyle \frac{2\pi m_x}{d}}-k_{0}sin{\theta }_{0}cos{\varphi }_0\right){\displaystyle \frac{d{N}_x}{2}}\right]\hfill \\ & \times & {\mathrm{sinc}}^2\left[\left(k_y-{\displaystyle \frac{2\pi m_y}{d}}-k_{0}sin{\theta }_{0}sin{\varphi }_0\right){\displaystyle \frac{d{N}_y}{2}}\right].\hfill \end{array}$$

(17)

where $g_{\alpha }=2\pi m_{\alpha }/d$, $\alpha =x,y$, are the components of the reciprocal lattice vector $\mathbf{g}$. Note that different values of $\mathbf{g}$ correspond to different Brillouin zones. In order to perform the integration over the angles ${\theta }_0$ and ${\varphi }_0$, we introduce the variables $v_x=(k_x-g_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)/2$ and $v_y=(k_y-g_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)/2$, which leads to

$$\begin{array}{ccc}\hfill \Gamma (k_x,k_y)& =& {\displaystyle \frac{3{\Gamma }_0}{\pi k_0^2}}{N}_x{N}_y\sum _{m_x,m_y}\int d{v}_x\int d{v}_y{\displaystyle \frac{{\mathrm{sinc}}^2\left(v_{x}d{N}_x\right){\mathrm{sinc}}^2\left(v_{y}d{N}_y\right)}{\sqrt{1-C^2(v_x,v_y)}}}\hfill \\ & \times & \left\{1-{\left[{\widehat{d}}_x{C}_x\left(v_x\right)+{\widehat{d}}_y{C}_y\left(v_y\right)+{\widehat{d}}_z\sqrt{1-C^2(v_x,v_y)}\right]}^2\right\}.\hfill \end{array}$$

(18)

We have here defined

$$\begin{array}{ccc}\hfill C_x\left(v_x\right)& =& (k_x-g_x-2v_x)/k_0,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill C_y\left(v_y\right)& =& (k_y-g_y-2v_y)/k_0,\hfill \end{array}$$

(20)

so that $C^2(v_x,v_y)=C_x^2\left(v_x\right)+C_y^2\left(v_y\right)<1$. The final expression forthe rate then depends on the lattice size, $(N_x,N_y)$, and we first discuss the case of an infinite 2D system.

3.2. Infinite Square Array

In the limit $N_x,N_y\to \infty$, using the property ${\mathrm{sinc}}^2\left(v_{x,y}d{N}_{x,y}\right)\to (\pi /d{N}_{x,y})\delta \left(v_{x,y}\right)$ in the rate expression (18) leads to

$$\begin{array}{ccc}\hfill \Gamma (k_x,k_y)& =& {\displaystyle \frac{3{\Gamma }_0\pi }{{\left(k_{0}d\right)}^2}}\sum _{m_x,m_y}{\displaystyle \frac{1}{\sqrt{1-C^2(0,0)}}}\left(1-{\left[{\widehat{d}}_x{C}_x\left(0\right)+{\widehat{d}}_y{C}_y\left(0\right)+{\widehat{d}}_z\sqrt{1-C^2(0,0)}\right]}^2\right),\hfill \\ \hfill \end{array}$$

(21)

with the condition $C^2(0,0)<1$, i.e., ${(k_x-g_x)}^2+{(k_y-g_y)}^2<k_0^2$: Each term of the sums in Eq. (21) is different from zero inside a circle of radius $k_0$ and center in $\mathbf{g}$. This allows us to recover the result of Ref. [15] for parallel polarization, with the dipoles oriented in the array plane (${\widehat{d}}_z=0$):

$$\begin{array}{ccc}\hfill {\Gamma }^{\Vert }\left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\Gamma }_0\pi }{{\left(k_{0}d\right)}^2}}\sum _{\mathbf{g}}{\displaystyle \frac{k_0^2-{[\widehat{\mathbf{d}}·(\mathbf{k}-\mathbf{g})]}^2}{\sqrt{k_0^2-{|\mathbf{k}-\mathbf{g}|}^2}}},\hfill \end{array}$$

(22)

as well as for perpendicular polarization, with the dipoles orthogonal to the 2D array plane (${\widehat{d}}_x={\widehat{d}}_y=0$):

$$\begin{array}{ccc}\hfill {\Gamma }^{\perp }\left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\Gamma }_0\pi }{{\left(k_{0}d\right)}^2}}\sum _{\mathbf{g}}{\displaystyle \frac{{|\mathbf{k}-\mathbf{g}|}^2}{\sqrt{k_0^2-{|\mathbf{k}-\mathbf{g}|}^2}}},\hfill \end{array}$$

(23)

with $\mathbf{g}=(g_x,g_y)$ the reciprocal array vector introduced previously. Figure 1 shows $\Gamma (k_x,k_y)/{\Gamma }_0$ as a function of $k_{x}d$ and $k_{y}d$ for $d={\lambda }_0/5$ and $d={\lambda }_0$ (${\lambda }_0=2\pi /k_0$ is the atomic transition wavelength), and for dipoles oriented along the x-axis. The black areas correspond to dark states, with $\Gamma (k_x,k_y)=0$: They have a suppressed emission, akin to subradiant modes, yet reaching $\Gamma (k_x,k_y)=0$ in this infinite-size limit. Note that for $d<{\lambda }_0/2$ only the terms $g_x=0$ and $g_y=0$ contribute to the sums in Eq. (21), while for ${\lambda }_0/2<d<{\lambda }_0$, the terms $g_x=\pm 2\pi /d$ and $m_y=\pm 2\pi /d$ also contribute.

3.3. Finite Square Array

The collective decay rates for a finite square are obtained from integral (17), for finite $N_x$ and $N_y$. A particular case is that of the mode $(k_x,k_y)=(0,0$), which can be addressed, for example, using a laser with wave vector perpendicular to the atomic plane, as done in Ref. [22]. Figure 2 shows $\Gamma (0,0)/{\Gamma }_0$ as a function of the normalized step $k_{0}d$, for a square array of $N\times N=10\times 10$ atoms, with dipole polarization $\widehat{\mathbf{d}}=(0,0,1)$ for panel (a) and $\widehat{\mathbf{d}}=(1,0,0)$ for panel (b).

The dashed red lines are the solution for an infinite lattice, see Eq. (21), and we observe that they provide a good estimate of the decay rates, apart from the strongly superradiant or subradiant ones. In particular, the rate of the symmetric mode goes to $\Gamma (0,0)=N^2{\Gamma }_0$ for a vanishing step $k_{0}d\to 0$ (out of scale in the figure): This corresponds to the superradiant case in the subwavelength regime, which scales as the number of particles [1].

As for subradiance, for perpendicular polarization [panel (a) in Figure 2], the collective decay rate is close to zero for subwavelength array step, $0<k_{0}d<2\pi$ (i.e., $d<{\lambda }_0$). Differently, for parallel polarization [panel (b)] the decay rate is less than the single-atom decay one ${\Gamma }_0$ only for $\pi <k_{0}d<2\pi$, and its value is overall much larger. The origin of this much more efficient suppression of the emission in the perpendicular polarization channel can be found in the pattern of the dipole radiation, which is in the atomic plane in that case so the light remains confined in the system, whereas for parallel polarization the dipole radiation has a substantial component out of the atomic array, thus promoting stronger leaks of the radiation out of the lattice.

3.4. Large-N Limit for the Finite Square Array

For a large but finite lattice ($1\ll N_x,N_y<\infty$), and considering modes along the x direction ($k_y=0$) and with dipoles perpendicular to the 2D array ($\widehat{\mathbf{d}}=\widehat{z}$), an approximate expression for the subradiant region of the spectrum $\Gamma (k_x,k_y)$ can be derived. Indeed, assuming $k_x>k_0$ and $d<{\lambda }_0/2$ in the integral (18) we obtain, for $\mathbf{g}=0$ and up to terms $O(1/N_{x,y}^2)$ (see Appendix A):

$$\begin{array}{ccc}\hfill \Gamma (k_x,0)& =& {\displaystyle \frac{3\pi {\Gamma }_0}{2{\left(k_{0}d\right)}^3\sqrt{k_{x}d}}}\left\{k_{x}d\sqrt{N_x}{\displaystyle \frac{sin\left[\frac{1}{2}arctan\left(\frac{1}{v_0}\right)\right]}{{(1+v_0^2)}^{1/4}}}-{\displaystyle \frac{4\sqrt{2}}{\sqrt{N_x}}}\sqrt{{\displaystyle \frac{v_0+\sqrt{1+v_0^2}}{1+v_0^2}}}\right\},\hfill \\ \hfill \end{array}$$

(24)

where $v_0=(d{N}_x/4k_x)(k_x^2-k_0^2)$. In the limit $d{N}_x\gg 4k_x/{(k_x^2-k_0)}^2$ (that is, far from the subradiant boundary $|k_x|=k_0$), Eq. (24) can be approximated by

$$\begin{array}{ccc}\hfill \Gamma (k_x,0)& \approx & {\displaystyle \frac{6\pi {\Gamma }_0}{{\left(k_{0}d\right)}^3{N}_x}}{\displaystyle \frac{k_x(2k_0^2-k_x^2)}{{(k_x^2-k_0^2)}^{3/2}}},\hfill \end{array}$$

(25)

which is valid for $k_0<\left|k_x\right|<\sqrt{2}{k}_0$. Hence, in the subradiant region ($|k_x|>k_0$) the rate scales as $\Gamma (k_x,0)\sim 1/N_x$. In the superradiant region $|k_x|<k_0$, for large $N_x$ the expression is close to the solution for an infinite chain (21), which diverges for $|k_x|=k_0$. Indeed, for $|k_x|=k_0$ and large $N_x$, Eq. (24) approximates as $\Gamma (k_0,0)\approx 3\pi {\Gamma }_0\sqrt{N_x/{\left(2k_{0}d\right)}^3}$.

These features are illustrated in Figure 3, where the exact solution of $\Gamma (k_x,0)/{\Gamma }_0$ from Eq. (17) [continuous blue line] is compared with the approximate solution (24) [dashed red line] as a function of $k_{x}d$, in the region $0<k_x<\pi /d$. The dash-dotted black line is the infinite chain solution (21), with the vertical dotted line marking the value $k_x=k_0$.

3.5. Finite Square Array: Radial Mode Distribution

For an infinite square array with perpendicularly oriented dipoles [$\widehat{\mathbf{d}}=(0,0,1)$], the rate $\Gamma (k_x,k_y)$ presents a radial symmetry around $(m_x,m_y)=(0,0)$ – see Eq.(23). Let us exploit this symmetry, now for a finite array, with $N_x=N_y=N\gg 1$. Substituting $k_x=k_{\perp }cos\theta$, $k_y=k_{\perp }sin\theta$, $v_{x}dN=vcos{\theta }^{\prime }$ and $v_{y}dN=vsin{\theta }^{\prime }$ into Eq. (18), we then use the approximation ${\mathrm{sinc}}^2\left(v_{x}dN\right){\mathrm{sinc}}^2\left(v_{y}dN\right)\approx 1/(1+v^4/4)$. For a subwavelength step $d<{\lambda }_0/2$, it leads to

$$\begin{array}{ccc}\hfill \Gamma \left(k_{\perp }\right)& =& {\displaystyle \frac{3{\Gamma }_0}{\pi {\left(k_{0}d\right)}^2}}{\int }_0^{\infty }{\displaystyle \frac{v}{1+v^4/4}}{\int }_0^{2\pi }{\displaystyle \frac{C^2(v,\theta )}{\sqrt{1-C^2(v,\theta )}}}d\theta dv,\hfill \end{array}$$

(26)

where $C(v,\theta )=\sqrt{{(k_{\perp }/k_0)}^2-4(k_{\perp }/k_0)(v/k_{0}dN)cos\theta +{(2v/k_{0}dN)}^2}$.

Figure 4(a) shows the behavior of $\Gamma \left(k_{\perp }\right)$ for $d={\lambda }_0/4$ and different values of N, from a numerical integration of Eq. (26). If we now introduce $v^{\prime }=v/N$ in Eq. (26), it now takes the form

$$\begin{array}{ccc}\hfill \Gamma \left(k_{\perp }\right)& =& {\displaystyle \frac{3{\Gamma }_0{N}^2}{\pi {\left(k_{0}d\right)}^2}}{\int }_0^{\infty }{\displaystyle \frac{v^{\prime }}{1+N^4{v}^{\prime 4}/4}}{\int }_0^{2\pi }{\displaystyle \frac{C^2(v^{\prime },\theta )}{\sqrt{1-C^2(v^{\prime },\theta )}}}d\theta d{v}^{\prime }\hfill \end{array}$$

(27)

where $C(v^{\prime },\theta )=\sqrt{{(k_{\perp }/k_0)}^2-4(k_{\perp }/k_0)(v^{\prime }/k_{0}d)cos\theta +{(2v^{\prime }/k_{0}d)}^2}$ is independent on N. This allows us to deduce that in the large N limit, Eq. (27) scales as $1/N^2$. This analysis is confirmed by the numerical analysis of Eq. (26), shown in Figure 4(b), where the $1/N^2$ trend is observed already for $N\ge 10$.

4. 3D Cubic Atom Array

Let us now consider a 3D cubic array, with atoms at positions ${\mathbf{r}}_j=d[(j_x-1){\widehat{\mathbf{e}}}_x+(j_y-1){\widehat{\mathbf{e}}}_y+(j_z-1){\widehat{\mathbf{e}}}_z]$ with $j_x=1,\dots N_x$, $j_y=1,\dots ,N_y$ and $j_z=1,\dots ,N_z$. The rate of mode $\mathbf{k}(k_x,k_y,k_z)$ (11) now writes

$$\begin{array}{ccc}\hfill \Gamma \left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\Gamma }_0}{2N_x{N}_y{N}_z}}{〈\left[1-{(\widehat{\mathbf{d}}·{\widehat{\mathbf{k}}}_0)}^2\right]{\left|F\left(\mathbf{k}\right)\right|}^2〉}_{\Omega },\hfill \end{array}$$

(28)

where

$$\begin{array}{ccc}\hfill {\left|F\left(\mathbf{k}\right)\right|}^2& =& {\displaystyle \frac{{sin}^2[(k_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)d{N}_x/2]}{{sin}^2[(k_x-k_{0}sin{\theta }_{0}cos{\varphi }_0)d/2]}}\hfill \\ & \times & {\displaystyle \frac{{sin}^2[(k_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)d{N}_y/2]}{{sin}^2[(k_y-k_{0}sin{\theta }_{0}sin{\varphi }_0)d/2]}}\times {\displaystyle \frac{{sin}^2[(k_z-k_{0}cos{\theta }_0)d{N}_z/2]}{{sin}^2[(k_z-k_{0}cos{\theta }_0)d/2]}}.\hfill \end{array}$$

(29)

Assuming a large array, $N_x,N_y,N_z\gg 1$, and using the same transformation as in Eq. (18) for 2D arrays, the rate rewrites

$$\begin{array}{ccc}\hfill \Gamma \left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\Gamma }_0}{2\pi k_0^2}}{N}_x{N}_y{N}_z\sum _{\mathbf{g}}\int d{v}_x\int d{v}_y{\mathrm{sinc}}^2\left(v_{x}d{N}_x\right){\mathrm{sinc}}^2\left(v_{y}d{N}_y\right)\hfill \\ & \times & {\displaystyle \frac{1-{\left[{\widehat{d}}_x{C}_x\left(v_x\right)+{\widehat{d}}_y{C}_y\left(v_y\right)+{\widehat{d}}_z\sqrt{1-C^2(v_x,v_y)}\right]}^2}{1-\sqrt{1-C^2(v_x,v_y)}}}\hfill \\ & \times & \left\{{\mathrm{sinc}}^2\left[\left(k_z-g_z-k_0\sqrt{1-C^2(v_x,v_y)}\right){\displaystyle \frac{d{N}_z}{2}}\right]\right.\hfill \\ & +& \left.{\mathrm{sinc}}^2\left[\left(k_z-g_z+k_0\sqrt{1-C^2(v_x,v_y)}\right){\displaystyle \frac{d{N}_z}{2}}\right]\right\}\hfill \end{array}$$

(30)

where $C^2(v_x,v_y)<1$. Here, $\mathbf{g}=(2\pi /d)\mathbf{m}$ is the reciprocal lattice vector, with $\mathbf{m}=(m_x,m_y,m_z)\in {\mathbb{Z}}^3$.

4.1. Infinite 3D Array

In the infinite array limit, $N_x,N_y,N_z\to \infty$, one can use the property ${\mathrm{sinc}}^2\left(v_{\alpha }d{N}_{\alpha }\right)\to (\pi /d{N}_{\alpha })\delta \left(v_{\alpha }\right)$ to simplify rate (30) as

$$\begin{array}{ccc}\hfill \Gamma \left(\mathbf{k}\right)& =& {\displaystyle \frac{3{\pi }^2{\Gamma }_0}{2\pi k_0^2{d}^3}}\sum _{\mathbf{g}}\left\{1-{\left[\widehat{\mathbf{d}}·\left(\mathbf{k}-\mathbf{g}\right)\right]}^2\right\}\delta \left[\left|\mathbf{k}-\mathbf{g}\right|-k_0\right].\hfill \end{array}$$

(31)

Each term of the sum in Eq. (31) is thus an infinitely-thin shell of center $\mathbf{g}$ and radius $k=k_0$.

We note that in Ref. [15] subradiant modes in infinite 1D and 2D arrays have been described as "guided", since the associated radiation fields are evanescent in the directions transverse to the array. In contrast, in 3D arrays the radiation modes are extended over space, in all directions. Thus, differently from the 1D and 2D cases, the subradiant modes in finite 3D arrays cannot be identified using the results of the infinite arrays, as it has been possible in 1D and 2D arrays. This has lead to focusing on the numerical diagonalization of the scattering matrix to determine the eigenvalues in the finite-N case [15]. In our approach, subradiant modes in a finite 3D array have a finite width around the surface $k=k_0$, as we will discuss in the next section.

4.2. Decay Rate for Mode Along an Axis of the Finite Cubic Array

Considering a cubic array with a large size, $N_x,N_y,N_z\gg 1$, an approximate solution for the decay rate $\Gamma (k_x,k_y,k_z)$ can be derived from Eq. (30). In the subradiant region of the spectrum, $k>k_0$ with $\mathbf{g}=0$ and considering for instance modes with polarization along the z axis and along the $k_x$ one ($\mathbf{k}=(k_x,0,0)$), we derive, up to terms $O(1/N_{x,y,z}^2)$ (see Appendix B):

$$\begin{array}{ccc}\hfill \Gamma (k_x,0,0)& =& {\displaystyle \frac{3\pi {\Gamma }_0{N}_x}{2{\left(k_{0}d\right)}^2}}{\mathrm{sinc}}^2\left({\displaystyle \frac{d{N}_x(k_x-k_0)}{2}}\right).\hfill \end{array}$$

(32)

Note that for $k_x=k_0$, using that $N=N_x{N}_y{N}_z$ is the total number of atoms in the cubic array with edges $L_y=N_{y}d$ and $L_z=d{N}_z$, the rate rewrites $\Gamma (k_x,0,0)=2b_0{\Gamma }_0$ where $b_0=(3\pi /4)N/\left(k_0^2{L}_y{L}_z\right)$ is the resonant optical thickness. In contrast, $\Gamma (k_x,0,0)$ decreases as $1/N_x$ when $k_x$ moves away from $k_0$, both inside and outside the shell, toward the dark zone. In Figure 5 the rate $\Gamma (k_x,0,0)/{\Gamma }_0$ as a function of $k_{x}d$ shows the very good agreement between the exact expression (28) and the approximate one (32).

In conclusion, for an infinite 3D array subradiance appears almost everywhere, except when $|\mathbf{k}-\mathbf{g}|=k_0$. For a finite cubic array, subradiant states are close to those of the infinite array, but with the sharp surfaces of the spectrum turned into smooth ones, with tails in the dark regions whose widths decrease as the inverse of the atom number N.

5. Conclusion

We have thus characterized the decay rates of 2D and 3D atomic arrays, identifying families of Bloch states with strong subradiance properties. In particular, we have shown that using an integral representation to deal with the sum over the many emitters allows one to obtain explicit expressions for the rates, but also to study finite-size effects in different dimensions. With the recent realization of a subradiant 2D atomic array mirror [22], the development of techniques to study systems too large to be studied directly by numerical diagonalization is a crucial task.

While the present work addressed the single-excitation regime, the possibility to use two-level atoms to manipulate multiple-photons states calls for the extension of such methods to the many-body regime [33,34,35]. In particular, combining the long lifetimes of subradiant modes with the photon-photon interactions achieved through dipolar or Rydberg interactions could allow for the creation of tunable atomic platforms for photon manipulation.