Cooperative Decay of an Ensemble of Atoms in a One-Dimensional Chain with a Single Excitation

Abstract

We propose a new expression of the cooperative decay rate of a one-dimensional chain of N two-level atoms in the single-excitation configuration. From it, the interference nature of superradiance and subradiance arises naturally, without the need to solve the eigenvalue problem of the atom–atom interaction Green function. The cooperative decay rate can be interpreted as the imaginary part of the expectation value of the effective non-Hermitian Hamiltonian of the system, evaluated over a generalized Dicke state of N atoms in the single-excitation manifold. Whereas the subradiant decay rate is zero for an infinite chain, it decreases as $1/N$ for a finite chain. A simple approximated expression for the cooperative decay rate is obtained as a function of the lattice constant d and the atomic number N. The results are obtained first for the scalar model and then extended to the vectorial light model, assuming all the dipoles aligned.

1. Introduction

Cooperative spontaneous emission by N excited two-level atoms has been extensively studied since the seminal works by Dicke in 1954 [1] and Lehmberg in 1970 [2]. Whereas superradiance, i.e., enhanced spontaneous emission due to constructive interference between the emitters, has been well understood [3,4], subradiance, i.e., inhibited emission due to destructive interference between the emitters, is more elusive and difficult to observe [5,6,7]. Nevertheless, the number of studies on subradiance has seen a large increase in the last few years, as it offers the opportunity of storing photons in emitter ensembles for times longer than the single emitter lifetime [8,9,10,11,12,13,14,15].

In disordered systems, the cooperative decay must be studied numerically, usually by solving the dynamics of an initially excited ensemble [16]. Eventually, in the linear regime of weak excitation, useful information can be obtained from the spectrum of the eigenvalues of the system [17], determined numerically by diagonalizing the finite matrix associated with the Green operator describing the coupling between the emitters. Things are apparently more simple when atoms form ordered arrays, which are easier to treat theoretically. This is also the reason why propagation of excitations in lattices with different dimensionalities has been generally investigated more than in disordered systems. In fact, the existing literature on the interaction of coherent light with atomic lattices is very rich [18,19,20,21,22]. Most of it presents a ‘solid-state physics’ point of view, investigating the transport of particular photonic (but not necessarily) modes through ordered samples, often taken as infinite but finite in a few exceptions. For instance, infinite and finite chains of two-level atoms have been considered in ref. [12]. Conversely, cooperative effects as superradiance and subradiance are typically interpreted as quantum optics phenomena, and are seen associated with cooperativity in the emission or light scattering of ensemble of atoms [3,8]. The bridge between these two complementary descriptions is not yet completely established.

The aim of this paper is to add some further insight on the study of a one-dimensional array of atoms in the single-excitation configuration, i.e., with only one atom excited among N. Starting from the effective non-Hermitian Hamiltonian, which includes an imaginary part describing the cooperative spontaneous decay and a real part describing the cooperative energy shift [2,23], we focus on the cooperative decay only, without solving the eigenvalue spectrum problem, but adopting a different approach: we introduce a collective decay function ${\mathsf{\Gamma }}_k$, parametrized by a continuous label k which, in the case of an infinite chain, represents the Fourier transform of the system. For a finite chain, the collective decay ${\mathsf{\Gamma }}_k$ is able to catch the main characteristics of the superradiant and subradiation emission by the excited atoms in the chain, enlightening the symmetry properties of the N-atoms state associated with such cooperative phenomena. In particular, starting from the solution for an infinite chain, obtained as a limit case of the general solution, we obtain an approximated analytical expression for the cooperative decay from a finite chain of N atoms when N is large. The analysis is carried out, initially assuming the scalar model of light and neglecting the vectorial nature of the dipoles. The scalar model is particularly attractive because, since the polarization direction does not play any role, it is able to catch the main features of cooperativity just considering the relative phases of the emitters, taken at a fixed distance along the chain. So, the control parameters are, in this case, only the lattice constant d and the atom number N. Then, we extend the results to the vectorial light model for a set of N equally oriented dipoles. We see that the scalar light model still provides a good approximation of the vectorial light model, either when the directions of the dipoles are randomized or in the limit of large lattice constant.

2. Scalar Model

We consider N two-level atoms with the same atomic transition frequency ${\omega }_0=c{k}_0$, linewidth $\mathsf{\Gamma }$ and dipole $\mu$. The atoms are prepared in a single-excitation state; $|g_j\rangle$ and $|e_j\rangle$ are the ground and excited states, respectively, of the j-th atom, $j=1,\dots ,N$, which is placed at position ${\mathbf{r}}_j$. We consider here the single-excitation effective Hamiltonian in the scalar approximation, whereas the exact vectorial model will be considered in the next section. If we assume that only one photon is present, when tracing over the radiation degrees of freedom, the dynamics of the atomic system can be described by the non-Hermitian Hamiltonian [24,25]

$$\begin{array}{ccc}\hfill \widehat{H}& =& -i{\displaystyle \frac{\hslash }{2}}\sum _{j,m}{G}_{jm}{\widehat{\sigma }}_j^{†}{\widehat{\sigma }}_m,\hfill \end{array}$$

(1)

where ${\widehat{\sigma }}_j=|g_j\rangle \langle e_j|$ and ${\widehat{\sigma }}_j^{†}=|e_j\rangle \langle g_j|$ are the lowering and raising operators, and $G_{jm}$ is the scalar Green function,

$$G_{jm}=\left\{\begin{array}{cc}{\mathsf{\Gamma }}_{jm}-i{\Omega }_{jm}\hfill & \mathrm{if}j\ne m,\hfill \\ \mathsf{\Gamma }\hfill & \mathrm{if}j=m,\hfill \end{array}\right.$$

(2)

and

$${\mathsf{\Gamma }}_{jm}=\mathsf{\Gamma }{\displaystyle \frac{sin\left(k_0{r}_{jm}\right)}{k_0{r}_{jm}}},{\Omega }_{jm}=\mathsf{\Gamma }{\displaystyle \frac{cos\left(k_0{r}_{jm}\right)}{k_0{r}_{jm}}},$$

(3)

where $r_{jm}=|{\mathbf{r}}_j-{\mathbf{r}}_m|$. $\widehat{H}$ contains both real and imaginary parts, which take into account that the excitation is not conserved, since it can leave the system by emission. We focus our attention on the decay term ${\mathsf{\Gamma }}_{jm}$. It can be obtained as the angular average of the radiation field propagating between the two atomic positions ${\mathbf{r}}_j$ and ${\mathbf{r}}_m$, with wave-vector $\mathbf{k}=k_0(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta )$,

$${\mathsf{\Gamma }}_{jm}={\displaystyle \frac{\mathsf{\Gamma }}{2}}{〈e^{-i\mathbf{k}·({\mathbf{r}}_j-{\mathbf{r}}_m)}+\mathrm{c}.\mathrm{c}.〉}_{\Omega }$$

(4)

where the angular average is defined as

$${〈f(\theta ,\varphi )〉}_{\Omega }={\displaystyle \frac{1}{4\pi }}{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }sin\theta f(\theta ,\varphi )d\theta .$$

Equation (4) provides a simple interpretation of ${\mathsf{\Gamma }}_{jm}$ as the coupling between the jth atom and the mth atom, mediated by the photon shared between the two atoms and averaged over all the vacuum modes [26,27,28]. The interaction also includes the energy shift ${\Omega }_{jm}$, which we will not consider here, since it is not relevant for our study. The important technical point is that Equation (4) allows us to factorize ${\mathsf{\Gamma }}_{jm}$ in the product of two terms, before averaging them over the total solid angle.

Things become particularly simple if we consider N atoms placed along a linear chain with lattice constant d, i.e., with positions ${\mathbf{r}}_j=d(j-1){\widehat{\mathbf{e}}}_z$, with $j=1,\dots ,N$. Then, we can write

$$e^{-i\mathbf{k}·{\mathbf{r}}_j}=E_j=e^{-i{k}_{0}d(j-1)cos\theta }$$

(5)

and

$${\mathsf{\Gamma }}_{jm}={\displaystyle \frac{\mathsf{\Gamma }}{2}}{〈E_j{E}_m^{\ast }+\mathrm{c}.\mathrm{c}.〉}_{\Omega }.$$

(6)

2.1. Collective Decay Rate

We introduce the following collective variable,

$${\mathsf{\Gamma }}_k={\displaystyle \frac{1}{N}}\sum _{j=1}^N\sum _{m=1}^N{\mathsf{\Gamma }}_{jm}{e}^{ikd(j-m)}$$

(7)

depending on the continuous index $k\in [0,2\pi /d)$. Equation (7) represents the Fourier transform of ${\mathsf{\Gamma }}_{jm}$ only in the case of the infinite chain, such that the system is periodic and k discrete; nevertheless, the continuous variable k is still a good label for the modes when N is sufficiently large [12]. We refer to ${\mathsf{\Gamma }}_k$ as the ‘continuous spectrum’ of the decay rate, outlining that we are not referring to the discrete spectrum of the eigenvalues of the system, described by the non-Hermitian Hamiltonian (1), but we introduce a quantity related to the phase of the collective state of the N atoms emitting a single photon. A more precise interpretation of it will be provided at the end of the section. For the moment, let us just proceed in its evaluation. By using Equation (6) in Equation (7), we can write

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k& =& {\displaystyle \frac{\mathsf{\Gamma }}{N}}\langle |F_k\left(\theta \right){{|}^2\rangle }_{\Omega }\hfill \end{array}$$

(8)

where

$$\begin{array}{ccc}\hfill |F_k{\left(\theta \right)|}^2& =& {\left|\sum _{j=1}^N{e}^{i(k-k_{0}cos\theta )d(j-1)}\right|}^2={\displaystyle \frac{{sin}^2[(k-k_{0}cos\theta )dN/2]}{{sin}^2[(k-k_{0}cos\theta )d/2]}}\hfill \end{array}$$

(9)

and

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k& =& {\displaystyle \frac{\mathsf{\Gamma }}{4\pi N}}{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }sin\theta {\left|F_k\left(\theta \right)\right|}^{2}d\theta ={\displaystyle \frac{\mathsf{\Gamma }}{k_{0}dN}}{\int }_{(k-k_0)d/2}^{(k+k_0)d/2}{\displaystyle \frac{{sin}^2\left(Nt\right)}{{sin}^{2}t}}dt\hfill \end{array}$$

(10)

where we changed the integration variable from $\theta$ to $t=(k-k_{0}cos\theta )d/2$. For large N, we can approximate in the integral of Equation (10),

$${\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],$$

(11)

where $\mathrm{sinc}\left(x\right)=sinx/x$, so that

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k& =& ={\displaystyle \frac{\mathsf{\Gamma }N}{k_{0}d}}\sum _{m=-\infty }^{+\infty }{\int }_{(k-k_0)d/2}^{(k+k_0)d/2}{\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]dt.\hfill \end{array}$$

(12)

2.2. Infinite Chain

In the limit $N\to \infty$,

$${\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]\to {\displaystyle \frac{\pi }{N}}\delta (t-m\pi )$$

(13)

where $\delta \left(x\right)$ is the Dirac delta function. Thus,

$${\mathsf{\Gamma }}_k={\displaystyle \frac{\mathsf{\Gamma }\pi }{k_{0}d}}\sum _{m=-\infty }^{+\infty }\Pi [2m\pi -k_{0}d<kd<2m\pi +k_{0}d]$$

(14)

where $\Pi (a<x<b)$ is the rectangular function, equal to 1 for $a<x<b$ and 0 elsewhere. Equation (14) is the solution for an infinite chain, where k is the true index of the modes. The fact that ${\mathsf{\Gamma }}_k$ is a Fourier transform in the continuous variable k is a consequence of Bloch’s theorem. Although the result (14) is not new (see, for instance, ref. [12], obtained analytically using the vectorial model, but only for $k_{0}d<\pi$), it will be useful for our study for two reasons:

(a)

We aim to compare the scalar model with the vectorial model and see what are the most relevant differences; the scalar model is more appealing, since it neglects the vectorial nature of the dipoles composing the chain, focusing only on their phases. Also, the scalar model is not completely unrealistic, since we will see that it can be recovered from the vectorial model assuming the dipoles randomly oriented.

(b)

We are interested in describing a finite chain, with a finite number N of atoms, where the infinite chain represents a limit case of it. Many previous studies of infinite and finite chains (usually based on the numerical evaluation of the complex eigenvalues of the matrix $G_{jm}$ in Equation (1)) have a ‘solid-state physics’ approach [12], investigating the guided propagation of photonic modes along the chain. On the contrary, we are interested to the point of view of the atoms, where superradiance and subradiance arise from constructive and destructive interferences of N emitters, respectively. In the paper by Bettles et al. [11], these two points of view are equally well discussed.

From the result (14) we see that, when $k_{0}d<\pi$ (i.e., $d/{\lambda }_0<0.5$), ${\mathsf{\Gamma }}_k=0$ for $k_0<k<2\pi /d-k_0$, i.e., full subradiance (see Figure 1a), whereas for $0<k<k_0$ and $2\pi /d-k_0<k<2\pi /d$, ${\mathsf{\Gamma }}_k=\mathsf{\Gamma }(\pi /k_{0}d)>\mathsf{\Gamma }$, i.e., enhanced radiance [12,14]. For larger lattice constant, slight subradiance and enhanced radiance alternate every half-wavelength, i.e., at increasing interval $\pi$ of $k_{0}d$ (see Figure 1b–d): it results then that, when $m\pi <k_{0}d<(m+1)\pi$ in the interval $(m+1)\pi <(k+k_0)d<(m+2)\pi$, ${\mathsf{\Gamma }}_k=\mathsf{\Gamma }(m+1)\pi /k_{0}d>\mathsf{\Gamma }$ for m odd and ${\mathsf{\Gamma }}_k=\mathsf{\Gamma }m\pi /k_0<\mathsf{\Gamma }$ for m even. In the external intervals $(k+k_0)d<\pi$ and $(k+k_0)d>(m+2)\pi$, ${\mathsf{\Gamma }}_k=\mathsf{\Gamma }m\pi /k_{0}d<\mathsf{\Gamma }$ for m odd and ${\mathsf{\Gamma }}_k=\mathsf{\Gamma }(m+1)\pi /k_{0}d>\mathsf{\Gamma }$ for m even. Hence, for each subsequent $\pi$ interval of $k_{0}d$, two regions of the spectrum can be identified: one where ${\mathsf{\Gamma }}_k$ is less than $\mathsf{\Gamma }$, and the other one where ${\mathsf{\Gamma }}_k$ is larger than $\mathsf{\Gamma }$. The difference between these two values is $\Delta \mathsf{\Gamma }=\pi /k_{0}d$. The only case of full subradiance, with zero decay rate, occurs for $k_{0}d<\pi$.

2.3. Finite Chain

The case of a finite chain requires the evaluation of the integral in Equation (12). Figure 2 shows ${\mathsf{\Gamma }}_k/\mathsf{\Gamma }$ vs. $k_{0}d$ for $k=0$ (blue line) and $kd=\pi$ (red line), for $N=10$ (Figure 2a) and $N=100$ (Figure 2b). The dashed line in Figure 2a is the analytic solution (14) for infinite chain. We see that, for an infinite chain and $k_{0}d<2\pi$, ${\mathsf{\Gamma }}_{k=0}=\mathsf{\Gamma }(\pi /k_{0}d)$, whereas for $\pi <k_{0}d<2\pi$, ${\mathsf{\Gamma }}_{k=\pi /d}=\mathsf{\Gamma }(2\pi /k_{0}d)$. We now have the elements to discuss superradiance and subradiance in a finite chain in detail. Let us consider the blue line of Figure 2a,b, corresponding to $k=0$. For $k_{0}d\ll 1/N$, ${\mathsf{\Gamma }}_{k=0}\sim \mathsf{\Gamma }N$, which corresponds to the Dicke superradiance of N atoms confined in a region smaller than the optical wavelength, such that they emit all in phase. Increasing $k_{0}d$, the distance between the atoms increases and they become less correlated: the cooperative emission rate decreases until, at $k_{0}d=\pi$, it reaches the uncorrelated value $\mathsf{\Gamma }$ (see Figure 2a). Then, beyond $k_{0}d=\pi$, the interference between the emitters becomes destructive and ${\mathsf{\Gamma }}_{k=0}$ drops below $\mathsf{\Gamma }$ (‘slight’ subradiance), until $k_{0}d$ reaches the value $2\pi$. At this value of the lattice constant, ${\mathsf{\Gamma }}_{k=0}$ has a rapid jump from the minimum $\mathsf{\Gamma }(\pi /k_{0}d)$ to the maximum $\mathsf{\Gamma }(3\pi /k_{0}d)$ (see Figure 2b, blue line). Then, the process repeats itself, with a jump of $\mathsf{\Gamma }(m+1)/m$ each time $k_{0}d$ crosses the value $m\left(2\pi \right)$.

Conversely, subradiance is described by the red line of Figure 2a,b, with $kd=\pi$: when $k_{0}d\ll 1/N$, the nearest neighbor atoms emit with opposite phases, and spontaneous emission is inhibited. Surprisingly, this destructive interference is also preserved for larger values of the lattice constant, up to $k_{0}d=\pi$. Then, beyond $k_{0}d=\pi$, we observe an alternation of enhanced and inhibited emission, mirroring the behavior of ${\mathsf{\Gamma }}_{k=0}$. Although full subradiance, ${\mathsf{\Gamma }}_k=0$, has a transparent interpretation as modes guided without radiative losses beyond the light line $k=k_0$ [12], the interpretation ‘from the point of view of the atoms’, in terms of destructive interference between the emitters in a finite chain, is less intuitive, as discussed, for instance, in ref. [11]. To this aim, we need to consider the superradiant and subradiant states of the system, as originally defined by Dicke [1]. Before treating them in the next subsection, it is useful to obtain further analytic expressions of the collective decay rate ${\mathsf{\Gamma }}_k$ for the finite chain. An approximated analytic expression can be obtained replacing the function ${\mathrm{sinc}}^2\left(x\right)$ in Equation (12) with the Lorentzian function $1/(1+x^2)$, which has the same peak value of unity and the same normalization value of $\pi$. Then, the integral in Equation (12) yields

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k& =& {\displaystyle \frac{\mathsf{\Gamma }}{k_{0}d}}\sum _{m=-\infty }^{+\infty }\left[arctan\left(b_m\right)-arctan\left(a_n\right)\right]\hfill \end{array}$$

(15)

where $a_m=[(k-k_0)d/2-m\pi ]N$ and $b_m=[(k+k_0)d/2-m\pi ]N$. Equation (15) approximates the behavior obtained from the exact result of Equation (12) well, as shown in Figure 3 for $k=\pi /d$ and $N=10$.

This expression allows to evaluate the subradiant decay rate in the limit of large N, using the identity $arctan\left(z\right)=\pm \pi /2-arctan(1/z)$, where the positive sign is for $z>0$ and the negative sign for $z<0$. Then, for $k_{0}d<\pi$ and the band interval $k_0<k<2\pi /d-k_0$, the main contribution in the sum in Equation (15) comes for the terms $m=0$ and $m=1$ which, for large N and far from the interval edges, yields

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k& \approx & {\displaystyle \frac{4\mathsf{\Gamma }}{N}}\left[{\displaystyle \frac{1}{{\left(kd\right)}^2-{\left(k_{0}d\right)}^2}}+{\displaystyle \frac{1}{{(kd-2\pi )}^2-{\left(k_{0}d\right)}^2}}\right].\hfill \end{array}$$

(16)

We observe the dependence on $1/N$, typical of subradiance [6]. In particular, at the center of the band, $k=\pi /d$, and for $k_{0}d=\pi /2$, ${\mathsf{\Gamma }}_k\approx (32/3{\pi }^{2}N)\mathsf{\Gamma }$. Conversely, superradiance is obtained from Equation (15) for $k=0$ and $m=0$,

$${\mathsf{\Gamma }}_{k=0}={\displaystyle \frac{2\mathsf{\Gamma }}{k_{0}d}}arctan\left({\displaystyle \frac{k_{0}dN}{2}}\right).$$

(17)

In the limit $k_{0}d\ll 1/N$, ${\mathsf{\Gamma }}_{k=0}\approx \mathsf{\Gamma }N[1-{\left(k_{0}dN\right)}^2/12+\dots ]$, as can be observed in Figure 2a.

2.4. Symmetric and Anti-Symmetric States

It is easy to prove that

$${\mathsf{\Gamma }}_k=-{\displaystyle \frac{2}{\hslash }}\mathrm{Im}\langle k|\widehat{H}|k\rangle ,$$

(18)

where $\widehat{H}$ is the non-Hermitian Hamiltonian (1) and

$$|k\rangle ={\displaystyle \frac{1}{\sqrt{N}}}\sum _{j=1}^N{e}^{ikd(j-1)}|j\rangle ,$$

(19)

where $k\in [0,2\pi /d)$ and $|j\rangle =|g_1,\dots ,e_j,\dots ,g_N\rangle$. For $k=0$, the state (19) corresponds to the Dicke state [1], whereas for $k=k_0$ is the timed Dicke state introduced by Scully and coworkers [29], corresponding to the entangled state of N atoms, where a single photon with momentum $\hslash k_0$ along the axis chain has been absorbed. For $kd=\pi$ the state (19) is subradiant, since nearest-neighbor atoms have opposite phases. This picture is in agreement with the analysis of ref. [11], where the nearest-neighbor phase difference ${\varphi }_{i+1}^{\ell }-{\varphi }_i^{\ell }$ for $i=1,\dots ,N$ is determined for the N eigenmodes of the system, with index $\ell =1,\dots ,N$ ordered from the largest to the smaller decay rate. In our case, the nearest-neighbor phase difference for the state $|k\rangle$ is $kd$, which provides a simple interpretation of enhancement or inhibition of spontaneous emission by N in-phase or out-of-phase atoms in the chain. It is interesting to observe that, for an infinite chain and in the scalar model, the imaginary part of the expectation value $\langle k\left|\widehat{H}\right|k\rangle$ is constant. This is not true for a finite chain, but it is a very good approximation for a sufficiently large N (see Figure 1a). Moreover, the fact that ${\mathsf{\Gamma }}_k$ is independent on k for the infinite chain is a consequence of the scalar model, for which the spontaneous emission in the vacuum modes is isotropic. This is not true for the vectorial model, as we will see in the next section.

It is interesting to note that the states $|k\rangle$ form a complete basis of the single-excitation manifold. In fact,

$${\displaystyle \frac{d}{2\pi }}{\int }_0^{2\pi /d}dk|k\rangle \langle k|=\sum _{j=1}^N|j\rangle \langle j|=1.$$

(20)

As expected, the states $|k\rangle$ are not orthogonal for a finite chain, since

$$\langle k^{\prime }|k\rangle ={\displaystyle \frac{sin[(k-k^{\prime })dN/2]}{sin[(k-k^{\prime })d/2]}}{e}^{i(k-k^{\prime })d(N-1)/2},$$

(21)

but they become so for an infinite chain, $\langle k^{\prime }|k\rangle \to \delta (k-k^{\prime })$ for $N\to \infty$.

2.5. Spectrum and Eigenvalues

Finally, it is interesting to see how ${\mathsf{\Gamma }}_k$ approximates the eigenvalues spectrum of a finite chain. Figure 4 shows ${\mathsf{\Gamma }}_k$, as calculated from Equation (12), and the N eigenvalues ${\lambda }_i$ of the $N\times N$ matrix ${\mathsf{\Gamma }}_{jm}$, ordered from the largest to the smallest, as function of $kd=\pi (i-1/2)/N$, with $i=1,\dots ,N$ for $k_{0}d=\pi /2$, $N=10$ (a) and $N=50$ (b). We see that ${\mathsf{\Gamma }}_k$ rather satisfactorily reproduces the features of the eigenvalues of the matrix ${\mathsf{\Gamma }}_{jm}$ and describes the transition from superradiance to subradiance going from $kd=0$ to $kd=\pi$, a transition that becomes sharper by increasing N.

3. Vectorial Model

We now extend the previous analysis to the vectorial model, taking into account the polarization of the electromagnetic field. The non-Hermitian Hamiltonian is now

$$\begin{array}{ccc}\hfill \widehat{H}& =& -i{\displaystyle \frac{\hslash }{2}}\sum _{\alpha ,\beta }\sum _{j,j^{\prime }}{G}_{\alpha ,\beta }({\mathbf{r}}_j-{\mathbf{r}}_{j^{\prime }}){\widehat{\sigma }}_{j,\alpha }^{†}{\widehat{\sigma }}_{j^{\prime },\beta }.\hfill \end{array}$$

(22)

where $\alpha ,\beta =(x,y,z)$. Here, ${\widehat{\sigma }}_{j,x}=({\widehat{\sigma }}_j^{m_J=1}+{\widehat{\sigma }}_j^{m_J=-1})/2$, ${\widehat{\sigma }}_{j,y}=({\widehat{\sigma }}_j^{m_J=1}-{\widehat{\sigma }}_j^{m_J=-1})/2i$ and ${\widehat{\sigma }}_{j,z}={\widehat{\sigma }}_j^{m_J=0}$, where ${\widehat{\sigma }}_j^{m_J}=|g_j\rangle \langle e_j^{m_J}|$ is the lowering operator between the ground state $|g_j\rangle$ and the three excited states $|e_j^{m_J}\rangle$ of the jth atom with quantum numbers $J=1$ and $m_J=(-1,0,1)$. The vectorial Green function in Equation (22) is [2,13]

$$G_{\alpha ,\beta }\left(\mathbf{r}\right)={\displaystyle \frac{3\mathsf{\Gamma }}{2}}{\displaystyle \frac{e^{i{k}_{0}r}}{i{k}_{0}r}}\left[{\delta }_{\alpha ,\beta }-{\widehat{n}}_{\alpha }{\widehat{n}}_{\beta }+\left({\delta }_{\alpha ,\beta }-3{\widehat{n}}_{\alpha }{\widehat{n}}_{\beta }\right)\left({\displaystyle \frac{i}{k_{0}r}}-{\displaystyle \frac{1}{k_0^2{r}^2}}\right)\right]$$

(23)

with $r=\left|\mathbf{r}\right|$ and ${\widehat{n}}_{\alpha }$ being the components of the unit vector $\widehat{\mathbf{n}}=\mathbf{r}/r$. We consider the linear chain with lattice constant d, i.e., ${\mathbf{r}}_j=d(j-1){\widehat{\mathbf{e}}}_z$, with $j=1,\dots ,N$, and all the dipoles aligned with an angle $\delta$ with respect to the chain’s axis, so that ${\widehat{n}}_{\alpha }={\widehat{n}}_{\beta }=cos\delta$ and

$$G_{\alpha ,\alpha }\left(\mathbf{r}\right)={\displaystyle \frac{3\mathsf{\Gamma }}{2}}{\displaystyle \frac{e^{i{k}_{0}r}}{i{k}_{0}r}}\left[{sin}^2\delta +(1-3{cos}^2\delta )\left({\displaystyle \frac{i}{k_{0}r}}-{\displaystyle \frac{1}{k_0^2{r}^2}}\right)\right].$$

(24)

Notice that if the dipoles are randomly oriented, ${\langle G_{\alpha ,\beta }\left(\mathbf{r}\right)\rangle }_{\Omega }=0$ for $\alpha \ne \beta$ and in Equation (24) ${\langle {cos}^2\delta \rangle }_{\Omega }=1/3$, so that ${\langle G_{\alpha ,\alpha }\left(\mathbf{r}\right)\rangle }_{\Omega }=e^{i{k}_{0}r}/\left(i{k}_{0}r\right)$, i.e., we recover the scalar model. This, in general, is not true in disordered systems, where the short-range terms $1/r^2$ and $1/r^3$ play a role at large densities [30].

The decay rate for the vectorial model is given by the real part of $G_{\alpha ,\alpha }$,

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}^{\left(\delta \right)}\left(r_{jm}\right)& =& {\displaystyle \frac{3\mathsf{\Gamma }}{2}}\left[{sin}^2\delta j_0\left(k_0{r}_{jm}\right)+(3{cos}^2\delta -1){\displaystyle \frac{j_1\left(k_0{r}_{jm}\right)}{k_0{r}_{jm}}}\right]\hfill \end{array}$$

(25)

where $j_0\left(x\right)=sinx/x$ and $j_1\left(x\right)=sinx/x^2-cosx/x$ are the spherical Bessel functions of order $n=0$ and $n=1$. As before, it is possible to write ${\mathsf{\Gamma }}^{\left(\delta \right)}\left(r_{jm}\right)$ as angular average of the radiation field emitted between the two atoms. By using the identities

$$\begin{array}{ccc}\hfill {\langle e^{ixcos\theta }\rangle }_{\Omega }& =& {\displaystyle \frac{1}{4\pi }}{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }sin\theta e^{ixcos\theta }d\theta =j_0\left(x\right)\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\langle {cos}^2\theta e^{ixcos\theta }\rangle }_{\Omega }& =& {\displaystyle \frac{1}{4\pi }}{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }sin\theta {cos}^2\theta e^{ixcos\theta }d\theta =j_0\left(x\right)-2{\displaystyle \frac{j_1\left(x\right)}{x}}.\hfill \end{array}$$

(27)

we can write

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}^{\left(\delta \right)}\left(r_{jm}\right)& =& {\displaystyle \frac{3\mathsf{\Gamma }}{4}}\left\{{\displaystyle \frac{1+{cos}^2\delta }{2}}{\langle E_j{E}_m^{\ast }\rangle }_{\Omega }+{\displaystyle \frac{1-3{cos}^2\delta }{2}}{\langle {cos}^2\theta E_j{E}_m^{\ast }\rangle }_{\Omega }+\mathrm{c}.\mathrm{c}.\right\}\hfill \end{array}$$

(28)

where $E_j$ has been defined in the previous section. The collective decay rate of the dipoles oriented with the angle $\delta$ is

$${\mathsf{\Gamma }}_k^{\left(\delta \right)}={\displaystyle \frac{1}{N}}\sum _{j=1}^N\sum _{m=1}^N{\mathsf{\Gamma }}^{\left(\delta \right)}\left(r_{jm}\right)e^{ikd(j-m)}$$

(29)

with $k\in [0,2\pi /d]$. Then,

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{\left(\delta \right)}& =& {\displaystyle \frac{3\mathsf{\Gamma }}{4N}}\left[(1+{cos}^2\delta )\langle |F_k\left(\theta \right){|}^2{\rangle }_{\Omega }+(1-3{cos}^2\delta )\langle {cos}^2\theta |F_k\left(\theta \right){{|}^2\rangle }_{\Omega }\right]\hfill \end{array}$$

(30)

where $|F_k{\left(\theta \right)|}^2$ is defined in Equation (9). By evaluating the angular average and using the approximation (11),

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{\left(\delta \right)}& =& {\displaystyle \frac{3\mathsf{\Gamma }N}{2k_{0}d}}\sum _{m=-\infty }^{+\infty }{\int }_{(k-k_0)d/2}^{(k+k_0)d/2}\left[{sin}^2\delta +{\displaystyle \frac{1}{2}}(1-3{cos}^2\delta ){\displaystyle \frac{{(kd-2t)}^2-{\left(k_{0}d\right)}^2}{{\left(k_{0}d\right)}^2}}\right]\hfill \\ & \times & {\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]dt.\hfill \end{array}$$

(31)

In the limit $N\to \infty$, using (13),

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{\left(\delta \right)}& =& {\displaystyle \frac{3\mathsf{\Gamma }\pi }{2k_{0}d}}\sum _{m=-\infty }^{+\infty }\left\{{sin}^2\delta +{\displaystyle \frac{1}{2}}(1-3{cos}^2\delta ){\displaystyle \frac{{(kd-2\pi m)}^2-{\left(k_{0}d\right)}^2}{{\left(k_{0}d\right)}^2}}\right\}\hfill \\ & \times & \Pi [2m\pi -k_{0}d<kd<2m\pi +k_{0}d],\hfill \end{array}$$

(32)

which is the Fourier transform of the decay rates for an infinite chain. We observe that, for $k_{0}d<\pi$, ${\mathsf{\Gamma }}_k^{\left(\delta \right)}$ is still zero in the interval $k_0<k<2\pi /d-k_0$ (full subradiance). Instead, ${\mathsf{\Gamma }}_k^{\left(\delta \right)}$ is no more uniform in k as in the scalar model. In particular, for $k_{0}d<\pi$ and in the intervals $0<k<k_0$ and $2\pi /d-k_0<k<2\pi /d$,

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{\left(\delta \right)}& =& {\displaystyle \frac{3\mathsf{\Gamma }\pi }{2k_{0}d}}\left\{{sin}^2\delta +{\displaystyle \frac{1}{2}}(1-3{cos}^2\delta )\left({\displaystyle \frac{k^2}{k_0^2}}-1\right)\right\}\hfill \end{array}$$

(33)

Equation (33) is in agreement with the results of ref. [12]. In Figure 5, we compare the results of the vectorial and the scalar models for an infinite chain, for the same values of $k_{0}d$ as in Figure 1 and for $\delta =\pi /2$. We observe that increasing $k_{0}d$ the differences between the scalar and the vectorial models become less important, and the spectrum becomes more flat. Also, the behavior of ${\mathsf{\Gamma }}_k^{\left(\delta \right)}$ vs. $k_{0}d$ for $k=0$ and $k=\pi /d$, shown in Figure 6, is similar to that obtained from the scalar model (see Figure 2). We observe that ${\mathsf{\Gamma }}_k^{\left(\delta \right)}$ reduces to the scalar model’s value, ${\mathsf{\Gamma }}_k$, for the special angle $\delta =arccos(1/\sqrt{3})$.

Finally, as performed for the scalar model, we can obtain an approximated expression of ${\mathsf{\Gamma }}_k^{\left(\delta \right)}$ for a finite lattice, substituting the function ${\mathrm{sinc}}^2\left(x\right)$ with the Lorentzian function $1/(1+x^2)$ in Equation (31) and solving the integral. A long but straightforward calculation yields:

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_k^{\left(\delta \right)}(x,a)& =& {\displaystyle \frac{3\mathsf{\Gamma }}{2a}}\sum _{m=-\infty }^{+\infty }\left\{\left[{sin}^2\delta +{\displaystyle \frac{1}{2a^2}}(1-3{cos}^2\delta )[{(x-2\pi m)}^2-a^2]\right]\right.\hfill \\ & \times & (arctan{b}_m-arctan{a}_m)\hfill \\ & +& {\displaystyle \frac{2}{Na}}(1-3{cos}^2\delta )-{\displaystyle \frac{1}{N{a}^2}}(1-3{cos}^2\delta )(x-2\pi m)ln{\displaystyle \frac{1+b_m^2}{1+a_m^2}}\hfill \\ & -& \left.{\displaystyle \frac{2}{N^2{a}^2}}(1-3{cos}^2\delta )(arctan{b}_m-arctan{a}_m)\right\}\hfill \end{array}$$

(34)

where we defined $x=kd$, $a=k_{0}d$, $a_m=(x-2\pi m-a)N/2$ and $b_m=(x-2\pi m+a)N/2$. A numerical analysis shows that the approximated expression (34) reproduces satisfactorily the exact solution (31). We have verified the $1/N$ scaling of subradiance by evaluating numerically the exact solution (31) for $kd=\pi$ and $k_{0}d=\pi /2$, and two different polarization directions, $\delta =\pi /2$ and $\delta =0$, shown in Figure 7. The linear fit confirms the subradiant law ${\mathsf{\Gamma }}_{\mathrm{sub}}\sim 1/N$, independently from the polarization direction.

4. Conclusions

We have presented a different approach to the study of the cooperative decay in a one-dimensional chain of N atoms in the single-excitation configuration. We have defined a collective function ${\mathsf{\Gamma }}_k$, containing the phase information among the emitters. In particular, it reduces to the Fourier transform of the decay rates in the case of an infinite chain. For a finite chain, it gives an intuitive interpretation of superradiance and subradiance as constructive or destructive interference of the emission by N atoms in the chain. More specifically, ${\mathsf{\Gamma }}_k$ is proportional to the imaginary part of the expectation value of the non-Hermitian Hamiltonian operator, evaluated on a generalized Dicke state of N atoms with nearest-neighbor phase difference $kd$, where d is the lattice constant: superradiance and subradiance occurs for $k=0$ and $k=\pi /d$ (but only for $d/{\lambda }_0<1/2$), respectively. Our analytic results are in agreement with the numerical results of [11,21]. More importantly, this approach can be complementary to the study of the discrete eigenvalue problem, generally numerically limited by the size of N. We have derived an explicit expression of ${\mathsf{\Gamma }}_k$, both for the scalar and vectorial light models. From it, subradiance shows a dependence as $1/N$, similar to that predicted and observed in disordered systems [6,7] and in a 2D lattice [10]. This result seems to be at odds with the analysis of ref. [12], where a dependence as $1/N^3$ is found for the most subradiant eigenvalues. However, that scaling refers to discrete eigenvalues, whose behavior depends in general on the microscopic details, such as, for instance, the polarization of the atoms. On the contrary, ${\mathsf{\Gamma }}_k$ is a mean-field quantity, obtained as an expectation value of the non-Hermitian Hamiltonian operator, evaluated over a collective state with a well-defined phase relationship between nearest-neighbor atoms. In this sense, ${\mathsf{\Gamma }}_k$ could be associated with measurement of the average subradiant decay in a finite linear chain, rather than the less accessible eigenvalues of the non-Hermitian Hamiltonian. It will be interesting to extend this approach of the study of the cooperative decay from one-dimensional chains to two- and three-dimensional finite arrays, eventually also in the presence of disorder.