Subradiance Generation in a Chain of Two-Level Atoms with a Single Excitation

1. Introduction

Since the seminal paper by Dicke [1], collective spontaneous emission by an ensemble of two-level atoms has been studied by many authors [2,3,4]. In particular subradiance, i.e. inhibited emission due to destructive interference between the emitters, has recently received great attention, both in disordered systems as in a cloud [5,6,7,8], and in ordered systems, as in atomic chains or 2D and 3D lattices [9,10,11,12,13,14]. The theoretical studies on subradiance has been focused mostly in the study to the eigenvalues of the system [15,16] and in particular on the collective decay rate, which for subradiance is less than the single-atom decay $\Gamma$. Few attention has been devoted to the generation of subradiance, from the preparation of the excited atoms or when the atoms are driven by a laser.

In disordered system, subradiance has been investigated numerically and experimentally by switching off a continuous detuned laser driving a thermal cloud and calculating the decay rate of the fluorescence light intensity by reporting it in a semilog plot [5,6]. After the initial fast decay, subradiance manifests itself in a slowly decaying fluorescence intensity with a rate below the single atom decay rate. At first, the subradiant decay is not purely exponential, since several modes decay simultaneously. For longer times, it then ends up with a pure exponential decay (referred as subradiant decay) when only one long-lived mode dominates. A similar approach can be adopted for ordered system, as atoms in a linear chain. Such linear chains have been investigated theoretically by different authors [17,18,19,20,21,22,23,24]. In the single-excitation approximation, subradiance has been studied in the previous literature by calculating the eigenvalues of the system, determined numerically by diagonalizing the finite $N\times N$ matrix associated with the Green operator describing the coupling between the emitters.

Recently, we proposed a different method to study subradiance in a finite linear chain, based on the evaluation of the expectation value of the effective Hamiltonian on a generalized Dicke state, depending on a continuous variable k[14]. This has several advantages: 1) it allows to calculate the collective decay rate $\Gamma \left(k\right)$ as a function of a continuous parameter, mimic the exact Fourier spectrum; 2) The sum over N can be calculated exactly, making N a mere parameter of the system and no more the dimension of the matrix whose eigenvalues are to be evaluated; 3) the study of $\Gamma \left(k\right)$ allows us to identify the subradiant regions of the spectrum. The last point is important, because it can be used to detect how subradiance may be dynamically generated. By solving the coupled-dipole equations for N atoms with given initial conditions or in presence of a driving laser, it is possible to project the solution at a given time on the generalized Dicke state. The result gives the probability distribution of the state as a function of the continuous variable k. Then, by identifying the subradiant regions of the spectrum, it is possible to see if subradiance has been generated. In particular, we find which is the initial state maximally generating subradiance: this is useful in order to understand the symmetry properties of the subradiant state.

The paper is organized as follow. In the first part we review the main results of Ref.[14], defining the generalized Dicke state and adding the calculation of the collective frequency shift, both in the scalar and vectorial models. The second part is devoted to the generation of subradiance by suitable initial atomic excitations or when atoms are excited by a weak laser field. The dynamics of the system is investigated by solving the N coupled-dipole equations and projecting the single-particle state a the time t over the generalized Dicke state. We will see that it is possible to build a maximally subradiant state, which in the limit of an infinite chain gives a vanishing spontaneous decay rate. Finally, we interpret the results in terms of the fluorescence intensity spatial distribution.

2. Modelling Emission from a Chain of Atom with a Single Excitation

Here we define the collective frequency shift and decay rate as the real and imaginary part of the expectation value of the effective Hamiltonian over the generalized Dicke state. Preliminary results have been previously reported in Ref. [14]. The calculations are carried on first for the scalar model and then for the vectorial model.

2.1. Scalar Model

We consider N two-level atoms with the same atomic transition frequency ${\omega }_0=c{k}_0$, linewidth $\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, while the 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 following non-Hermitian Hamiltonian [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, $G_{jm}$ is the scalar Green function,

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

(2)

and

$${\Gamma }_{jm}=\Gamma {\displaystyle \frac{sin\left(k_0{r}_{jm}\right)}{k_0{r}_{jm}}},{\Omega }_{jm}=\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 takes into account that the excitation is not conserved since it can leave the system by emission. The real part 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 )$,

$${\Gamma }_{jm}={\displaystyle \frac{\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 ${\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. Equation (4) allows to factorize ${\Gamma }_{jm}$ in the product of two terms, before averaging them over the total solid angle. 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,\cdots ,N$, we can write

$${\Gamma }_{jm}={\displaystyle \frac{\Gamma }{4}}{\int }_0^{\pi }sin\theta \left[e^{-i{k}_{0}d(j-m)cos\theta }+\mathrm{c}.\mathrm{c}.\right]d\theta .$$

(5)

2.2. Generalized Dicke State.

We define the generalized Dicke states [14]

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

(6)

where $|j\rangle =|g_1\cdots e_j\cdots g_N\rangle$ and $k\in (-\pi /d,\pi /d)$. It includes for $k=0$ the Dicke state [1] and for $k=k_0$ the Timed-Dicke state introduced by Scully [26,27]. It satisfy the completeness relation

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

(7)

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},$$

(8)

but they become so for an infinite chain, $\langle k^{\prime }|k\rangle \to \delta (k-k^{\prime })$ for $N\to \infty$. So the states $|k\rangle$ form an over-completed basis the the single-excitation manifold.

2.3. Collective Frequency Shift and Decay Rate.

Taking the expectation value of the effective Hamiltonian (1) over the generalized Dicke state (6) yields

$$-{\displaystyle \frac{2}{\hslash }}\langle k|\widehat{H}|k\rangle ={\Omega }_N\left(k\right)+i{\Gamma }_N\left(k\right)$$

(9)

were

$${\Omega }_N\left(k\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N\sum _{{\displaystyle \genfrac{}{}{0pt}{}{m=1}{m\ne j}}}^N{\Omega }_{jm}{e}^{ikd(j-m)}$$

(10)

is the collective frequency shift and

$${\Gamma }_N\left(k\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N\sum _{m=1}^N{\Gamma }_{jm}{e}^{ikd(j-m)}$$

(11)

is the collective decay rate. By using Equation (5) in Equation (11) we can write

$$\begin{array}{ccc}\hfill {\Gamma }_N\left(k\right)& =& {\displaystyle \frac{\Gamma }{2N}}{\int }_0^{\pi }sin\theta {\left|F_k\left(\theta \right)\right|}^{2}d\theta \hfill \end{array}$$

(12)

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}$$

(13)

and

$$\begin{array}{ccc}\hfill {\Gamma }_N\left(k\right)& =& {\displaystyle \frac{\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}$$

(14)

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 (14),

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

(15)

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

$$\begin{array}{ccc}\hfill {\Gamma }_N\left(x\right)& =& ={\displaystyle \frac{\Gamma N}{a}}\sum _{m=-\infty }^{+\infty }{\int }_{(x-a)/2}^{(x+a)/2}{\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]dt.\hfill \end{array}$$

(16)

where $a=k_{0}d$ and $x=kd$, with $x\in (-\pi ,\pi )$. Hence, we have transformed the double sum in Equation (11) into an integral, where N plays the role of a simple parameter of the system. The collective frequency shift takes the form

$${\Omega }_N\left(x\right)={\displaystyle \frac{\Gamma }{N}}\sum _{j=1}^N\sum _{{\displaystyle \genfrac{}{}{0pt}{}{m=1}{m\ne j}}}^N{\displaystyle \frac{cos\left(a\right[j-m\left|\right)}{a|j-m|}}{e}^{ix(j-m)}.$$

(17)

2.4. Infinite Chain

If the chain is infinite, $N\to \infty$, the solution for the collective decay rate is

$$\begin{array}{ccc}\hfill {\Gamma }_{\infty }\left(x\right)& =& {\displaystyle \frac{\Gamma \pi }{a}}\sum _{m=-\infty }^{+\infty }\Pi [2m\pi -a<x<2m\pi +a]\hfill \end{array}$$

(18)

where $\Pi (a<x<b)$ is the rectangular function, equal to 1 for $a<x<b$ and 0 elsewhere. In the first Brillouin zone, $m=0$, ${\Gamma }_{\infty }\left(x\right)=\Gamma \pi /a$ for $\left|x\right|<a$ and ${\Gamma }_{\infty }\left(x\right)=0$ for $a<\left|x\right|<\pi$. Atomic modes in the region enclosed within the light line $k=\pm k_0$ are generally unguided and radiate into free space. Outside the light line ($\left|k\right|>k_0$), the modes are guided and subradiant, as the electromagnetic field is evanescent in the directions transverse to the chain. For an infinite chain, ${lim}_{N\to \infty }{\Omega }_N\left(x\right)={\Omega }_{\infty }\left(x\right)$ depends only on the index $\ell =j-m$,

$$\begin{array}{ccc}\hfill {\Omega }_{\infty }\left(x\right)& =& \Gamma \sum _{\ell =-\infty }^{\infty }{\displaystyle \frac{cos\left(a\right|\ell \left|\right)}{a|\ell |}}{e}^{ix\ell }={\displaystyle \frac{\Gamma }{2a}}\sum _{\ell =1}^{\infty }{\displaystyle \frac{1}{\ell }}\left[e^{i(a+x)\ell }+e^{i(a-x)\ell }+\mathrm{c}.\mathrm{c}.\right].\hfill \end{array}$$

(19)

Using the expansion

$$ln(1-z)=-\sum _{n=1}^{\infty }{\displaystyle \frac{z^n}{n}}$$

we write

$$\begin{array}{ccc}\hfill {\Omega }_{\infty }\left(x\right)& =& -{\displaystyle \frac{\Gamma }{2a}}\left\{ln[1-e^{i(a+x)}]+ln[1-e^{i(a-x)}]+\mathrm{c}.\mathrm{c}.\right\}\hfill \\ & =& -{\displaystyle \frac{\Gamma }{a}}ln\left[2|cosa-cosx|\right]\hfill \end{array}$$

(20)

The frequency shift has a logarithmic divergence for $x=\pm a$ and three extremes at $x=0,\pm \pi$, with ${\Omega }_{\infty }\left(0\right)=-(2\Gamma /a)ln\left[2\right|sin(a/2)\left|\right]$ and ${\Omega }_{\infty }(\pm \pi )=-(2\Gamma /a)ln\left[2\right|cos(a/2)\left|\right]$, respectively.

2.5. Finite Chain

If the chain is finite, the collective decay rate ${\Gamma }_N\left(x\right)$ can be calculated by using Equation (16). The collective frequency shift of Equation (17) can be written transforming the double sum in a single sum on the index $\ell =j-m$, with a degenerate factor $g_{\ell }=N-\ell$:

$$\begin{array}{ccc}\hfill {\Omega }_N\left(x\right)& =& {\displaystyle \frac{2\Gamma }{N}}\sum _{\ell =1}^{N-1}(N-\ell ){\displaystyle \frac{cos(a\ell )}{a\ell }}cos(x\ell )\hfill \end{array}$$

(21)

3. Vectorial Model

We now extend the previous expressions 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

$$G_{\alpha ,\beta }\left(\mathbf{r}\right)={\displaystyle \frac{3\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,\cdots ,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^{\left(\delta \right)}\left(r_{jm}\right)={\displaystyle \frac{3\Gamma }{2}}{\displaystyle \frac{e^{i{k}_0{r}_{jm}}}{i{k}_0{r}_{jm}}}\left[{sin}^2\delta +(1-3{cos}^2\delta )\left({\displaystyle \frac{i}{k_0{r}_{jm}}}-{\displaystyle \frac{1}{k_0^2{r}_{jm}^2}}\right)\right].$$

(24)

where $r_{jm}=a|j-m|$. The decay rate for the vectorial model is given by the real part of $G^{\left(\delta \right)}\left(r_{jm}\right)$,

$$\begin{array}{ccc}\hfill {\Gamma }^{\left(\delta \right)}\left(r_{jm}\right)& =& {\displaystyle \frac{3\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$. The frequency shift is given by the negative of the imaginary part of $G^{\left(\delta \right)}\left(r_{jm}\right)$,

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

(26)

It is possible to demonstrate that the collective decay rate is [14]

$$\begin{array}{ccc}\hfill {\Gamma }_N^{\left(\delta \right)}\left(x\right)& =& {\displaystyle \frac{3\Gamma N}{2a}}\sum _{m=-\infty }^{+\infty }{\int }_{(x-a)/2}^{(x+a)/2}\left[{sin}^2\delta +{\displaystyle \frac{1}{2}}(1-3{cos}^2\delta ){\displaystyle \frac{{(x-2t)}^2-a^2}{a^2}}\right]\hfill \\ & \times & {\mathrm{sinc}}^2\left[\left(t-m\pi \right)N\right]dt\hfill \end{array}$$

(27)

while the collective frequency shift of the dipoles oriented with the angle $\delta$ is

$$\begin{array}{ccc}\hfill {\Omega }_N^{\left(\delta \right)}\left(x\right)& =& {\displaystyle \frac{2}{N}}\sum _{\ell =1}^{N-1}(N-\ell ){\Omega }^{\left(\delta \right)}(a\ell )cos(x\ell )\hfill \end{array}$$

(28)

If the chain is infinite, $N\to \infty$,

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

(29)

and

$$\begin{array}{ccc}\hfill {\Omega }_{\infty }^{\left(\delta \right)}\left(x\right)& =& 2\sum _{\ell =1}^{\infty }{\Omega }^{\left(\delta \right)}(a\ell )cos(x\ell )\hfill \\ & =& {\displaystyle \frac{3\Gamma }{2a^3}}\mathrm{Re}\left\{-a^2{sin}^2\delta \left[ln\left(1-e^{i(x+a)}\right)+ln\left(1-e^{i(x-a)}\right)\right]+(3{cos}^2\delta -1)\right.\hfill \\ & \times & \left.\left[-ia{\mathrm{Li}}_2\left(e^{i(x+a)}\right)+ia{\mathrm{Li}}_2\left(e^{i(x-a)}\right)+{\mathrm{Li}}_3\left(e^{i(x+a)}\right)+{\mathrm{Li}}_3\left(e^{i(x-a)}\right)\right]\right\}\hfill \end{array}$$

(30)

where ${\mathrm{Li}}_{\nu }\left(z\right)={\sum }_{\ell =1}^{\infty }{z}^{\ell }/{\ell }^{\nu }$ is the PolyLog function. Figure 1 and Figure 2 show $\Gamma \left(x\right)/\Gamma$ and $\Omega \left(x\right)/\Gamma$ vs x for $a=\pi /2$, for an infinite chain (dashed lines) and for a finite chain with $N=10$ (continuous lines). Black lines are for the scalar model, red and blue lines are for the vectorial model, with $\delta =0$ and $\delta =\pi /2$, respectively. For an infinite chain the collective decay rate is zero for $\left|x\right|>a$ (i.e. for $\left|k\right|>k_0$), both for the scalar and the vectorial model. Notice that for $\delta =0$ the phase shift ${\Omega }_{\infty }^{(\delta =0)}\left(x\right)$ is not diverging at $x=\pm a$.

4. Dynamics

Having characterized the properties of the collective decay rate and frequency shift, identifying the subradiant zone $\left|x\right|>a$ where the collective decay rate is less than the single-atom decay rate $\Gamma$, we are interested now to study how subradiance can be generated by properly exciting the atoms. First we define the probability distribution of the system to be in a given generalized Dicke state $|k\rangle$, expressed in term of the single-particle basis. Then we will study the time evolution of this distribution, leading to subradiance. The behavior in the vectorial model is not much different from that of the scalar model, so we will limit our discussion to the scalar model only.

4.1. The Probability Density $P\left(k\right)$

Let us assume that the state of the system is $|\Psi \rangle =\alpha |g_1\cdots g_N\rangle +|{\Psi }^{\prime }\rangle$, where

$$|{\Psi }^{\prime }\rangle =\sum _{j=1}^N{\beta }_j|j\rangle$$

(31)

describes the state in the single-excitation manifold. By projecting on the generalized Dicke basis (6),

$$|\Psi \rangle ={\displaystyle \frac{d}{2\pi }}{\int }_{-\pi /d}^{\pi /d}dk|k\rangle \langle k|{\Psi }^{\prime }\rangle$$

(32)

where

$$\langle k|{\Psi }^{\prime }\rangle ={\displaystyle \frac{1}{\sqrt{N}}}\sum _{j=1}^N{e}^{-ikd(j-1)}{\beta }_j={\displaystyle \frac{1}{\sqrt{N}}}{A}_N\left(k\right).$$

(33)

Hence, the probability density to be in a state $|k\rangle$ is

$$P\left(k\right)={\displaystyle \frac{dN}{2\pi }}{\displaystyle \frac{|\langle k|{\Psi }^{\prime }\rangle {|}^2}{\langle {\Psi }^{\prime }|{\Psi }^{\prime }\rangle }}={\displaystyle \frac{d}{2\pi }}{\displaystyle \frac{{\left|{\sum }_{j=1}^N{e}^{-ikd(j-1)}{\beta }_j\left(t\right)\right|}^2}{{\sum }_{j=1}^N{\left|{\beta }_j\right|}^2}}={\displaystyle \frac{|A_N{\left(k\right)|}^2}{{\int }_{-\pi /d}^{\pi /d}{\left|A_N\left(k\right)\right|}^{2}dk}}$$

(34)

with

$${\int }_{-\pi /d}^{\pi /d}P\left(k\right)dk=1.$$

(35)

Equation (34) expresses the probability density in terms of the dipole amplitudes of the single atoms, whose time evolution is described in the following section.

4.2. Dynamics of the probability density $P\left(k\right)$

Let consider the time evolution of the atomic system in the presence of an external driving field in the scalar model. In the linear regime, the probability amplitudes ${\beta }_j\left(t\right)$ evolve with the following coupled-dipole equations,

$${\displaystyle \frac{d{\beta }_j}{dt}}=i{\Delta }_0{\beta }_j-i{\displaystyle \frac{{\Omega }_0}{2}}{e}^{ia(j-1)}-{\displaystyle \frac{\Gamma }{2}}\sum _{m=1}^N{G}_{jm}{\beta }_m$$

(36)

where $G_{jm}$ is defined in Equation (2), $a=k_{0}d$, ${\Omega }_0$ and ${\Delta }_0$ are the Rabi frequency and the laser-atom detuning of the driving laser field. From Equation (36), it is possible to obtain the equation for the temporal evolution of the probability amplitude $A_N(x,t)$ defined in Equation (33) (where $x=kd$) (see Appendix A):

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial A_N(x,t)}{\partial t}}& =& \left(i{\Delta }_0-{\displaystyle \frac{\Gamma }{2}}\right)A_N(x,t)-i{\displaystyle \frac{{\Omega }_0}{2}}{\displaystyle \frac{sin\left[\right(x-a)N/2]}{sin\left[\right(x-a)/2]}}{e}^{-i(x-a)(N-1)/2}\hfill \\ & -& {\displaystyle \frac{\Gamma }{a}}\sum _{\ell =1}^{N-1}\left[sin(a\ell )-icos(a\ell )\right]{\displaystyle \frac{cos(x\ell )}{\ell }}{A}_{N-\ell }(x,t)\hfill \end{array}$$

(37)

For an infinite chain,

$$\underset{N\to \infty }{lim}{\displaystyle \frac{sin\left[\right(x-a)N/2]}{sin\left[\right(x-a)/2]}}=2\pi \delta (x-a)$$

and

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial A_{\infty }(x,t)}{\partial t}}& =& \left[i{\Delta }_0+i{\displaystyle \frac{{\Omega }_{\infty }\left(x\right)}{2}}-{\displaystyle \frac{{\Gamma }_{\infty }\left(x\right)}{2}}\right]A_{\infty }(x,t)-i\pi {\Omega }_0\delta (x-a)\hfill \end{array}$$

(38)

where ${\Gamma }_{\infty }$ and ${\Omega }_{\infty }$ are defined in Equations (18) and (20). The solution of Equation (38) is

$$\begin{array}{ccc}\hfill A_{\infty }(x,t)& =& A_{\infty }(x,0)e^{(i{\Delta }_0+i{\Omega }_{\infty }\left(x\right)/2-{\Gamma }_{\infty }\left(x\right)/2)t}\hfill \\ & +& {\displaystyle \frac{2\pi {\Omega }_0\delta (x-a)}{2{\Delta }_0+{\Omega }_{\infty }+i{\Gamma }_{\infty }}}\left(1-e^{(i{\Delta }_0+i{\Omega }_{\infty }\left(x\right)/2-{\Gamma }_{\infty }\left(x\right)/2)t}\right).\hfill \end{array}$$

(39)

5. Generation of Subradiance

Based on the previous expressions, we now discuss how subradiance can be generated studying the temporal evolution from some initial conditions of the probability amplitudes ${\beta }_j$. The case of the excitation by an incident laser field will be discussed in the next section.

5.1. Single Excited Atom

As a first example, we consider a chain of $N=100$ atoms with $a=\pi /2$, no driving laser, ${\Omega }_0=0$ and ${\Delta }_0=0$, and a single initial atom excited in the middle of the chain, with ${\beta }_{N/2}=1$ and all the others ${\beta }_j$ equal to zero. From Equation (34), the initial probability to be in the state $|x\rangle$ is $P(x,0)=1/\left(2\pi \right)$ (where $x=kd$) i.e. it is uniform. Figure 3 shows $P(x,t)$ vs x at different times, from $t=0$ until $t=10/\Gamma$: after few time units, the probability becomes zero for $\left|x\right|<a$ and different from zero in the subradiant interval $a<\left|x\right|<\pi$.

For an infinite chain, $|A_{\infty }{(x,t)|}^2=exp[-{\Gamma }_{\infty }\left(x\right)t]$ where, for $a<\pi$, ${\Gamma }_{\infty }\left(x\right)=\Gamma (\pi /a)$ for $\left|x\right|<a$ and ${\Gamma }_{\infty }\left(x\right)=0$ for $a<\left|x\right|<\pi$. Then,

$$P(x,t)={\displaystyle \frac{e^{-{\Gamma }_{\infty }\left(x\right)t}}{2\left[\pi -a+a{e}^{-(\Gamma \pi /a)t}\right]}}$$

(40)

In the limit $t\to \infty$, $P(x,t)\to 0$ for $\left|x\right|<a$ and $P(x,t)\to {\displaystyle \frac{1}{2(\pi -a)}}$ for $a<\left|x\right|<\pi$. Figure 4 shows $P(x,t)$ at $t=10/\Gamma$ and $a=\pi /2$ (blue continuous line), obtained for a finite chain with $N=100$, together the value $P(x,\infty )=1/\left[2\right(\pi -a\left)\right]$ obtained for infinite chain (dashed line).

Hence, a single excited atom generates a subradiance state with a probability $P\left(x\right)$ which for an infinite chain is uniform in the subradiance spectral region $a<\left|x\right|<\pi$. It is interesting to see the distribution of the dipole amplitudes for the case of Figure 4: Figure 5 shows $|{\beta }_j|$ vs j at $t=10/\Gamma$, where the inset shows the average excitation probability $\langle |\beta {|}^2\rangle$ vs time. The initial excitation for the atoms with $j=N/2$ spread among the adjacent atoms, to a final distribution generating subradiance. Also if not visible in the Figure 5, it is expected that the destructive interference among the atoms inhibit the spontaneous decay of the excitation, as it will be discussed in the next section.

5.2. The Most Subradiant State

Since an initial uniform probability $P\left(x\right)$ generates asymptotically a subradiant distribution which is uniform for an infinite chain, as observed in the previous case, we are now interested to obtain the values of ${\beta }_j$ which are generating such subradiant distribution. We assume

$$A_N\left(x\right)=e^{-ix(N/2-1)}\left\{\begin{array}{ccc}1\hfill & \mathrm{if}a<\hfill & \left|x\right|<\pi \hfill \\ 0\hfill & \mathrm{if}\hfill & \left|x\right|<a.\hfill \end{array}\right.$$

(41)

which describes a subradiant state, with zero probability distribution in the superradiant region $\left|x\right|<a$ and uniform distribution in the subradiant region $a<\left|x\right|<\pi$. Calculating the single atom probability amplitude ${\beta }_j$ we obtain

$${\beta }_j={\displaystyle \frac{1}{2\pi }}{\int }_{-\pi }^{\pi }{e}^{ix(j-1)}{A}_N\left(x\right)dx=\left\{\begin{array}{ccc}1-{\displaystyle \frac{a}{\pi }}\hfill & \mathrm{if}\hfill & j={\displaystyle \frac{N}{2}}\hfill \\ \hfill & \\ -{\displaystyle \frac{sin\left[a\right(j-N/2\left)\right]}{\pi (j-N/2)}}\hfill & \mathrm{if}\hfill & j\ne {\displaystyle \frac{N}{2}}\hfill \end{array}\right.$$

(42)

In Figure 6 we plot $|{\beta }_j|$ vs j for $N=100$ and $a=1$, at the initial time $t=0$ (blue circles) and at $\Gamma t=10$ (red squares), obtained solving numerically Equation (36) with the initial condition (42). The inset shows the average excitation probability $\langle |\beta {|}^2\rangle$ vs time, which remains approximately constant. The ${\beta }_j$ as defined in Equation (42) reproduce the distribution amplitude (41) only in the limit $N\to \infty$ (see Appendix B), since the states $|k\rangle$ are not hortogonal. The distribution shown in Figure 6 has the same features as that shown in Figure 5 obtained from a single initially excited atom. To better appreciate the similarities with the case of Section 5.1, Figure 7 shows the probability density distribution $P(x,t)$ at $t=0$ (black-dashed line) and at $\Gamma t=10$ (blue continuous line): for an infinite chain, $P\left(x\right)=1/\left[2\right(\pi -a\left)\right]$ for $a<\left|x\right|<\pi$ and zero for $\left|x\right|<a$ , as obtained from Equation (41) (red line in Figure 7). Notice the similarities between Figure 5 for the single initially excited atom and Figure 7 for the initial state (42). We conclude that the state described by Equation (42) represents the ’most subradiant’ state for a finite chain of N atoms, with a purely subradiant spectrum in the limit of an infinite chain.

5.3. Atoms Driven by a Laser.

We now study the subradiance generation when the atoms are excited by an external laser field and then switched off. As an example, we consider a chain of $N=100$ atoms with $a=\pi /2$, weakly driven by a detuned laser, with ${\Omega }_0=0.1\Gamma$ and ${\Delta }_0=10\Gamma$. The atoms, initially unexcited (i.e. ${\beta }_j\left(0\right)=0$ for $j=1,\cdots ,N$) are driven by the laser up to $t_0=50/\Gamma$, after which the laser is switched off. Figure 8 shows $P(x,t)$ vs x and at different times t after the laser is switched off.

The distributions $P(x,t)$ at the laser switch-off time $t=50/\Gamma$ (blue continuous line) and at $t=100/\Gamma$ (red line), are shown in Figure 9: we see that the driving laser on brings the atoms close to a Timed-Dicke state, $|k_0\rangle$[26] (dashed line $x=a$ in Figure 9), with a width inversely proportional to the chain’s length $dN$. The inset of Figure 9 shows that, after the initial fast decay, subradiance manifests itself in a slowly decaying. At first, the subradiant decay is not purely exponential, since several modes decay simultaneously. For longer times, it then ends up with a pure exponential decay when only one long-lived mode dominates. However, the precise evaluation of the decay rate can be problematic due to the general non-exponential decay. In our approach, we determine the precise distribution of the subradiant modes: as it can be observed from the red line of Figure 9), at later times after the laser switch-off the distribution $P(x,t)$ is mostly in the subradiant region, $x>a$. The distribution is broad, so there is not a single subradiant mode dominating. As a consequence, the decay is not purely exponential.

In the case of infinite chain, Equation (39) gives

$$|A_{\infty }{(x,\infty )|}^2={\displaystyle \frac{2N{\left(\pi {\Omega }_0\right)}^2}{{(2{\Delta }_0+{\Omega }_{\infty })}^2+{\Gamma }_{\infty }^2}}\delta (x-a),$$

(43)

where we used the relation

$$\underset{N\to \infty }{lim}{\displaystyle \frac{{sin}^2[(x-a)N/2]}{{sin}^2[(x-a)/2]}}=\pi N\delta (x-a).$$

Hence, for a driven infinite chain the asymptotic spectrum is $P(x,\infty )\propto \delta (x-a)$ and no subradiance occurs. To observe subradiance, we need a finite chain, as seen in Figure 8 and Figure 9: in the detuned case, the finite width, proportional to $1/N$, of the driving term (second term in r.h.s. of Equation (37) and blue line in Figure 9) is responsible for the subradiant components of the spectrum until the laser is on, which subsequently evolve without reaching a steady-state value. From the above analysis, the possibility to have access to the full spectral distribution of the subradiant modes is clearly more advantageous than observing the time decaying excitation to obtain the subradiant decay rate, as done for instance in Ref.[5,6].

We now consider again the same chain of $N=100$ atoms with $a=\pi /2$, but driven by a resonant laser, with ${\Omega }_0=0.1$ and ${\Delta }_0=0$. The result is rather surprising: Figure 10 shows the average excitation, $\langle |\beta {|}^2\rangle$ vs $\Gamma t$, where the drive field is switched off at $\Gamma t=50$ (vertical dashed line). The average excitation grows almost linearly when the laser is turned on, typical for a diffusive regime [28,29]. After the laser is off, the excitation decays very slowly, showing that the excitation remains trapped in the atomic chain. We can understand this peculiar behavior by observing the probability density $P(x,t)$ in Figure 11 at the laser switch-off time, $t=50/\Gamma$ (red line) and at $t=100/\Gamma$ (blue line). Contrarily to the detuned case, at resonance the subradiant region of the spectrum is already populated when the laser is on (red line of Figure 11). After the laser has been switched off, at $\Gamma t=100$, $P(x,t)$ is remained almost the same, with only the radiating components, for $x<a$, decayed (blue line of Figure 11). Since now the spectrum is completely in the subradiant region $x>a$, the decay rate is almost zero and the atoms remain excited for a sufficiently long time.

6. Radiated Intensity

The following important question arises: may the probability $P(x,t)$ be determined measuring the scattered intensity at a certain angle $\theta$ with respect to the chain’s axis? We know that the scattered field appears as a sum of spherical waves radiated by the atomic dipoles,

$$E_s(\mathbf{r},t)=-{\displaystyle \frac{\hslash \Gamma }{e\mu }}\sum _{j=1}^N{\displaystyle \frac{exp\left[i{k}_0\right|\mathbf{r}-{\mathbf{r}}_j\left|\right]}{k_0|\mathbf{r}-{\mathbf{r}}_j|}}{\beta }_j\left(t\right)$$

(44)

In the far-field limit, one has $k_0|\mathbf{r}-{\mathbf{r}}_j|\approx k_{0}r-\mathbf{k}·{\mathbf{r}}_j$, where $\mathbf{k}=k_0(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta )$ and ${\mathbf{r}}_j=d(j-1){\widehat{\mathbf{e}}}_z$, so the field of Equation (44) radiated in a direction $\mathbf{k}$ reads

$$E_s(\theta ,t)\approx -{\displaystyle \frac{\hslash \Gamma }{e\mu }}{\displaystyle \frac{e^{i{k}_{0}r}}{r}}\sum _{j=1}^N{e}^{-i{k}_{0}dcos\theta (j-1)}{\beta }_j\left(t\right)$$

(45)

and the scattered intensity is

$$I_s(\theta ,t)=I_1{\left|\sum _{j=1}^N{e}^{-i{k}_{0}dcos\theta (j-1)}{\beta }_j\left(t\right)\right|}^2\propto P(k_z,t)$$

(46)

where $k_z=k_{0}cos\theta$. Hence, the atoms radiate out of the chain’s axis for $|k_z|<k_0$. The subradiant region $|k_z|>k_0$ is not accessible by the scattered field, since it would be $cos\theta >1$ and the electromagnetic field is evanescent in the directions transverse to the chain, since $k_{\perp }=i{k}_0\sqrt{{cos}^2\theta -1}=i\xi$. Very few photons are emitted outside the chain’s axis direction (anyone in the case of an infinite chain). However, from the radiated intensity it is possible to see if the atomic state is subradiant or not, observing if the atoms are emitting in a direction out of the axis’ chain. This can be seen in Figure 12 and Figure 13, where we plot the field intensity $I_s\propto {\left|E_s\right|}^2$ (where $E_s$ is determined by Equation (44)) in the plane $x=5d$, emitted by a chain of $N=50$ atoms with $kd=1$, centered at $x=y=0$, for two different atomic distributions. Figure 12 shows a case of uniform excitation, with ${\beta }_j=1/\sqrt{N}$ such that $P\left(k_z\right)=N{\mathrm{sinc}}^2(k_{z}dN/2)$ and the probability distribution is peaked around $k_z=0$: in this case the atoms emit out of the chain’s axis. Figure 13 shows the emission from the most subradiant state, with ${\beta }_j$ given by Equation (42), such that $P\left(k_z\right)\approx 0$ for $k_z<k_0$: the field is evanescent and does not propagate out of the chain, remaining confined along the chain; since the chain is finite, most of the energy is radiated out at the ends of the chain [11].

7. Conclusions

In conclusion, we have discussed analytically and numerically how subradiance can emerge from the evolution of the dynamics of N two-level atoms in the single-excitation configuration along a linear chain. In the first part, we have characterized the spectrum of the decay rates and frequency shifts of the system, identifying the regions of the spectrum where spontaneous emission is enhanced or inhibited, up to a complete suppression in the case of an infinite chain. Although these results are already known in literature [11,13,14], few attentions have been dedicated to the study of the subradiance generation. We proceeded first by obtaining a relation between the spectrum of emission and the single-particle amplitudes, whose evolution can be determined by solving the coupled-dipole equations. Then we have studied how different initial excitations evolve toward a subradiant state. A single-excited atom leads to an almost uniform population of subradiant modes. This has suggested the idea that the atomic configuration leading to this uniform population can be calculated directly, obtaining what we named the ’most subradiant state’. Then, we investigated how subradiance may be induced by a driving laser, which excite the atoms and successively is switched off, such that the long-lived subradiant modes survive for long times. Finally, we found the relation between spontaneous emitted intensity and subradiance. Subradiance is characterized by a suppression of the emission in the direction transverse to the chain axis. This analysis may be useful to envisage strategies to detect subradiance in ordered systems by measuring the radiation out of the lattice. The results obtained here for a linear chain can be extended to 2D a 3D lattices.