Optical waveguiding by atomic entanglement in multilevel atom arrays

Spin Model for Multilevel Atoms.

Intuitively, the interaction with light allows for processes of photon-mediated emission and reabsorption, where the excitation of one atom decays and another is excited. We describe such dynamics by means of a spin model (46⇓⇓⇓⇓⇓⇓–53) where we integrate out the photons and find an all-atomic density matrix that depends only on the internal degrees of freedom of the atoms. Specifically, the evolution of the atomic density matrix ρ^A obeys ρ^̇A=−(i/ℏ)[H,ρ^A]+L[ρ^A], where H is the Hamiltonian, and L[ρ^A] is the Lindblad operator. For atoms with hyperfine structure and resonance frequency ω0, these operators readH=ℏ∑i,j=1N∑q,q′=−11Jijqq′Σ^iq†Σ^jq′,[3a]L[ρ]=∑i,j=1N∑q,q′=−11Γijqq′22Σ^jq′ρΣ^iq†−Σ^iq†Σ^jq′ρ−ρΣ^iq†Σ^jq′,[3b]where we have defined the polarization-dependent spin-exchange and decay rates asJijqq′=−μ0ω02ℏ|℘|2 e^q⋅Re G(ri,rj,ω0)⋅e^q′*,[4a]Γijqq′=2μ0ω02ℏ|℘|2 e^q⋅Im G(ri,rj,ω0)⋅e^q′*,[4b]with ℘=⟨Fg‖er^‖Fe⟩ being the reduced matrix element associated with the transition. In the above equations, e^±1=∓12(x^±iŷ), and e^0=z^ are spherical basis vectors. Physically, the spin-exchange and decay rates are proportional to the classical field amplitude projected along polarization q at position ri, due to a classical oscillating dipole of polarization q′ at rj. Naturally, they are given in terms of the free-space electromagnetic Green’s tensor, G(ri,rj,ω0)≡G(rij,ω0), with rij=ri−rj, which readsG(r,ω0)=eik0r4πk02r3(k02r2+ik0r−1)1+  +(−k02r2−3ik0r+3)r⊗rr2,[5]where r≡|r| and k0=2π/λ0=ω0/c is the wave number corresponding to the atomic transition energy. Note that the Green’s function Gαβ is a tensor quantity ({α,β}={x,y,z}), as both the electromagnetic field and the atomic transition have specific polarizations. These equations are derived within the Markovian approximation, which is highly justified for our system under consideration (13, 54, 55).

The dynamics under the master equation can analogously be described in the quantum jump formalism (56), where the last 2 terms in the parentheses of Eq. 3b are combined with H to form a non-Hermitian Hamiltonian that characterizes the deterministic evolution. Within this formalism, the effective non-Hermitian Hamiltonian is readily given byHeff=ℏ∑i,j=1N∑q,q′=−11Jijqq′−iΓijqq′2Σ^iq†Σ^jq′.[6]While this equation is quite general, we now highlight several points, which are key to understanding collective decay:

• • First, for a single, isolated atom, we obtain the total decay rate from excited sublevel me as Γ0=∑qΓiiqqCme+q,q2, which is equal for all me. The branching fraction into a specific ground state mg=me+q is simply given by Cme+q,q2.

• • Second, the form of Heff illustrates that, in general, lowering one atom angular-momentum projection by q does not necessarily imply that another atom’s projection is raised by q (q′≠q), as q refers to the dipole matrix element and not the global polarization of the emitted photon (Fig. 1C). Thus, a system initialized in a 2-level subspace does not generally remain in that space under coherent dynamics. In particular, even if the initially involved excited state decays to a unique ground state, a different excited state can become populated later in time (see atoms i and k in Fig. 1C). We note that one way to suppress such undesired excitation pathways (e.g., the σ+ transition of atom k) is by applying large Zeeman shifts. While this might indeed constitute a practical way to maintain the desired effects of subradiance to some extent given hyperfine structure, our focus here is on identifying a true many-body mechanism for subradiance, which fundamentally persists in the presence of multiple pathways.

• • Finally, we consider the situation where an excited atom has an allowed decay channel into a ground state that is unoccupied by any atoms (illustrated by atom k in Fig. 1C following absorption of a σ+ photon). In that case, there is no interference; the excited state decays into that ground-state level with the full strength of a single, isolated atom. Presumably, any mechanism for strong subradiance must then prevent this process from occurring.

It will also be useful to calculate the emitted quantum field associated with any given atomic state. It can be shown that the positive frequency component readsE^+(r)=μ0ω02∑j=1N∑q=−11G(r,rj,ω0)⋅e^q* ℘ Σ^jq.[7]From Eq. 7, we can find the negative-frequency component E^−(r) by taking the Hermitian conjugate of E^+(r). The field also in principle contains a vacuum noise component (25), which will not affect our quantities of interest and is thus not explicitly written.

One-Dimensional Chain.

From now on, we focus on the problem of a 1D chain. We will see that this system both encodes the classical waveguiding effect previously identified for 2-level atoms and allows for additional highly correlated subradiant states due to multilevel structure. Moreover, 1D allows for numerics on relatively large systems (not dominated by boundaries), and we can take advantage of additional conserved quantities that allow us to simplify the problem. For the sake of simplicity, we choose the direction of the chain to align to the quantization axis of the atoms (i.e., along z). This guarantees the conservation of polarization of the emitted photons, as Gαβ=0 along the axis if α≠β (i.e., Gαβ is a diagonal tensor in the polarization indexes). Therefore, the only nonzero interactions are those preserving polarization (i.e., q=q′) and the effective Hamiltonian simplifies toHeff=ℏ∑i,j=1N∑q=−11Jijq−iΓijq2Σ^iq†Σ^jq,[8]where Jijq≡Jijqq and Γijq≡Γiiqq.

Until now, our spin model captures the dynamics of atoms with any kind of hyperfine structure. Hereafter, we focus on what we consider to be the minimal toy model that captures all of the relevant physics: 6-level atoms with Fg=1/2 and Fe=3/2 (Fig. 1B). This specific hyperfine structure displays closed transitions, where an excited state with maximum angular-momentum projection me=±3/2 can decay only to only one ground state, also with maximum angular-momentum projection mg=±1/2 in the ground-state manifold by emitting a circularly polarized photon (q=∓1). It also exhibits excited states with me=±1/2 where decay into 2 different ground states mg=±1/2 is allowed, which involve emission of photons of different polarizations (q=±1 or q=0). To simplify notation, in what follows we label ground states with mg={−1/2,1/2} as {0,1}, respectively, and excited states with me={−3/2,−1/2,1/2,3/2} as {2,3,4,5}, respectively (Fig. 1B).

The conservation of angular-momentum projection along the direction of the chain allows us to diagonalize the Hamiltonian by blocks with well-defined Fz=∑i=1Nmi, where the sum over Zeeman sublevels includes both ground and excited states (Fig. 2). Note that subspaces with equal magnitude of angular-momentum projection but opposite sign (i.e., ±Fz) display the same energy and decay spectra (i.e., the physics is identical if we simultaneously flip mg to −mg and me to −me). In what follows, we study single-excitation subradiant states of various angular momentum manifolds, starting by maximum |Fz| (as we will show, the most “classical” manifold) and ending by Fz=0 (the most complex and where subradiant states enabled by rich entanglement live).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Classification of (single-excitation) basis states within each subspace of conserved angular-momentum projection Fz along the z axis. For |Fz|=Fzmax, the states live within a 2-level subspace, where one (N−1) atoms occupy the excited (ground) states of minimum me(mg). As |Fz| is lowered, one type of basis state consists of replacing some number of the ground-state atoms of minimum mg=−1/2 with atoms (shown in red) with mg=+1/2. These basis states containing “defect” atoms predominantly compose subradiant states, as described in the main text.

Motivated by previous work on 1D and 2D arrays with simple atomic structure (25), here we seek to elucidate the properties of eigenstates ψξ of Heff in the single-excitation manifold. In particular, such states will have complex eigenvalues ωξ=Jξ−iΓξ/2 characterizing the energy shifts and decay rates. Of specific interest is the identification of states where the decay rate approaches zero as N→∞, which implies that the states decouple from radiation fields and correspond to guided modes, and to elucidate the properties of the eigenstates that enable the waveguiding phenomenon. Importantly, while for simple 2-level atoms, the Hilbert space of the single-excitation manifold increases as ∼N and encodes classical linear optics, here the manifold size increases exponentially, raising the possibility for waveguiding through entanglement.

Phase-Separated States.

Decreasing angular momentum increases the number of defects, which tend to localize at the chain edges. In Fz=0, approximately half of the chain is made of defects, and phase separation and domain walls emerge, which we discuss in this section. Phase-separated states occur when the number of defects is so large that half of the chain is in a “defect state” with respect to the other half. The “active” phase involves one atom excited in 2 (or 5) and N/2−2 in 0 (1), which participate in 2-level subradiance: The atoms exchange photons with q=1(−1) and get excited and decay via closed transitions. The inactive or defect phase consists of N/2+1 atoms in 1 (0). It should be noted that the eigenstate is a spatial superposition, where the active part can equally occupy the left and right sides of the chain, with a domain wall in between the phases, as shown in Fig. 5A.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Phase-separated subradiant states in Fz=0. (A) Sketch of the most subradiant eigenstate. Blue (red) circles represent the phase with one excited atom in 2 (5) and N/2−2 atoms in 0 (1). The faded circles represent the “inactive” or defect phase. (B) Evidence of phase separation in a chain of N=14 atoms through correlations between 0 states at atoms i,j. (C) Scaling of decay rate with atom number for phase-separated states (orange line). In red is scaling of the most subradiant state of a chain of N/2−1 atoms in |Fzmax|. The solid line is a guide to the eye showing a scaling of Γ/Γ0∼1/(N/2−1)3. For B and C, d=0.3λ0.

To observe signatures of phase separation, we calculate connected correlation functions between different Zeeman sublevels m,n. These are defined asCmn(i,j)=⟨σ^mmiσ^nnj⟩−⟨σ^mmi⟩⟨σ^nnj⟩,[12]where the expectation value is calculated over the most subradiant eigenstate, and i,j are atomic indexes. Fig. 5B shows the correlations C00(i,j) between 0 states in different locations for the most subradiant state of a chain of N=14 atoms. We observe 0 states clustering around each other and “repelling” atoms in 1 states, displaying anticorrelations. The orange line in Fig. 5C shows the scaling of the decay rate with respect to atom number. Due to the exponentially large size of the Hilbert space, we can exactly diagonalize only up to N=14 atoms. However, based on our intuition of phase separation, we expect the decay rate to closely coincide with that obtained for the most subradiant eigenstate of a chain of N/2−1 atoms in the |Fzmax| subspace. This calculation is plotted in Fig. 5C (red curve) and seems to agree with the scaling of the phase-separated states.

As in our previous analysis of other manifolds, thus far these results do not point to some fundamentally new mechanism for waveguiding: These states can be understood as the chain separating into 2 domains, where 2-level subradiance occurs independently from the other.

Symmetric States.

We also numerically find evidence of another type of subradiant state, with significant population in states 3 and 4. These states also show no particular length scales or features in pairwise correlations, suggesting a qualitatively different mechanism for subradiance. To better grasp the underlying physics, we first introduce a “toy model” Hamiltonian, which is inspired by Eq. 9. In particular, while the Hamiltonian can exactly be written in terms of pure spin wave operators for an infinite system, we consider a hypothetical finite system whose Hamiltonian also takes the same form,H=ℏ∑q=−11∑kJk,q−iΓk,q2Ŝk,q†Ŝk,q,[13]where Ŝk,q†=N−1/2∑jeikdjΣ^jq† and Jk,q and Γk,q are given in Eq. 10. In particular, as these quantities correspond to the results for an infinite system, we have that Γk,q=0 when |k|>k0. The wavevector k is now taken to be a discrete index, with k=2πn/N (with n∈[0,N−1]) to ensure periodic boundary conditions.

A perfectly subradiant eigenstate ψk̃, of well-defined quasi-momentum k̃, with zero decay rate fulfills Im{Hψk̃}=0. It should be noted that despite the apparent simplicity of Eq. 13, it is in general challenging to diagonalize, as the involved operators have complicated commutation relations. To proceed, we first find a state that fulfills the less demanding condition of being a zero-eigenvalue eigenstate of all spin operators Ŝk,q if k≠k̃. Such a state is found to beψk̃=N∑j=1Neik̃zj2jD3/2j+β3jD1/2j  +β4jD−1/2j+5jD−3/2j,[14]where β=(C−1/2,1/C1/2,1), N is a normalization constant, andDαj=∑m≠jPm0⊗n0⊗1⊗n1[15]is a Dicke state where all of the atoms except j are in ground states 0 or 1 such that the total angular momentum projection of the Dicke state is α={±1/2,±3/2}. In the above expression, n0(1)=(N−1±2α)/2 is the number of atoms in state 0(1), and ∑mPm denotes the sum over all distinct permutations of the ground states. It should be noted that Dicke states are known to exhibit significant multipartite entanglement (57, 58). More importantly, the entanglement of this state cannot be obtained simply by applying a projection operator on a coherent state, as was the case for 2-level subradiance. Remarkably, the limit where H is composed of only one k vector corresponds to ref. 45, which shows the existence of states that are zeros of a single jump operator. Our states, however, must satisfy many relationships involving all spin wave operators Ŝk≠k̃,q.

For the state in Eq. 14 to be a lossless eigenstate of H it needs to fulfillHψk̃=ℏ∑q=−11Jk̃,q−iΓk̃,q2Ŝk̃,q†Ŝk̃,qψk̃=λk̃ψk̃,[16]with λk̃∈R. In Eq. 16, we have utilized the property that Sk≠k̃,qψk̃=0, to eliminate all but a single spin wave operator from Hψk̃. For |k̃|>k0 (such that Γk̃,q=0), we find that the eigenstate condition Eq. 16 is satisfied provided that the dispersion relations for both polarizations |q|=0,1 coincide at k̃; i.e., Jk̃,q≡Jk̃, with a corresponding eigenvalue λk̃=C−1/2,12ℏJk̃. We subsequently show that an intersection of the dispersion relations at |k̃|>k0 gives rise to waveguiding not only in our toy model, but also in the original physical system.

In free space the dispersion relations Jk,±1 and Jk,0 are in general different (25) and intersect for a given k̃>k0 only for distances d≲0.17λ0. Numerically, however, it is difficult to confirm directly that an intersection in such a case leads to a decay rate approaching zero as N→∞. We attribute this to the fact that for the limited N that we can simulate, the small lattice constant d≪λ0 still results in significant finite-size effects due to the short length of the chain and because the intersection of the dispersion relations occurs close to radiative wavevectors |k̃|∼k0 (that is, for a short system, significant components of the many-body wavefunction still have wavevectors that couple to radiation). Thus, to better confirm our hypotheses, we add an additional short-range term to the q=±1 interaction rates, Jij,q=±1→Jijq=±1+Jijq=±1′, given byJij,q=±1′(r)=−360Δ7dπr4.[17]Here, Δ=Jk=π/d,q=0−Jk=π/d,q=1 is the difference between the original free-space dispersion relations of Eq. 10, evaluated at the edge of the Brillouin zone k=π/d. It can readily be shown that this extra term guarantees that both the dispersion relations and the slopes of the q=0,±1 polarizations coincide at the Brillouin zone edges, |k|=π/d. In the following, we do not focus on the details of how such dispersion engineering is implemented, although we note that it can potentially be realized by introducing some dielectric structure to change the Green’s function itself or by dressing of atomic levels.

Fig. 6A shows the spatial profile of the most subradiant eigenstate of a finite chain of N=16 atoms obtained with dispersion engineering. The populations in levels 2 and 5 are identical, and so are the populations in 3 and 4. In Fig. 6B, we plot the decay rate of the most subradiant state vs. atom number. In particular, we observe a scaling of the decay rate with atom number of Γ/Γ0∼1/N3, identical to that observed for classical waveguiding. Fig. 6C shows the fidelity between the most subradiant eigenstate of the finite chain and the toy model state ψk̃ given by Eq. 14. There is extremely good agreement at k̃=π/d. In particular, the infidelity at k̃=π/d scales as ϵ=1−|⟨ψ|ψk̃⟩|2∼1/N2, which indicates that our toy model accurately captures the physics of the actual system. Finally, Fig. 6D shows the field intensity emitted by the most subradiant eigenstate of the chain with dispersion engineering. The intensity pattern is consistent with waveguiding along the chain, where most of the field is radiated through the edges. We note that these states are robust against classical fluctuations in the atomic positions, as shown in SI Appendix.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Symmetric states in the manifold of zero angular-momentum projection (Fz=0). (A) Spatial profile of the most subradiant eigenstate of a chain of N=16 atoms, where the |q|=1 component of the Green’s function has been altered according to Eq. 17. (B) Scaling of the decay rate as a function of atom number. The green circles correspond to the case where the |q|=1 component of the Green’s function has been altered according to Eq. 17, and the gray circles are obtained for the unaltered case. The solid line is a guide to the eye and scales as Γ/Γ0∼1/N3. (C) Overlap between the state ψk̃ given by Eq. 14 and the most subradiant state of a finite chain obtained by numerical diagonalization, as a function of the wavevector k̃. The gray line is obtained for a chain of N=12 atoms in free space, while the blue line is obtained for N=16 atoms in a medium with dispersion engineering given by Eq. 17. The red shaded area represents the radiative wavevectors. (D) Field intensity (arbitrary units) created by the most subradiant mode in a chain of N=16 atoms, obtained by including dispersion engineering. Red circles denote atomic positions. For all plots, d=0.3λ0.

We emphasize the importance of realizing an intersection of the dispersion relations Jk,±1 and Jk,0, to yield this class of highly entangled, waveguiding states. In particular, in the gray circles and curves in Fig. 6 B and C, respectively, we plot the corresponding results where only the free-space Green’s function is used, without the additional term of Eq. 17. It can be seen that although a moderate reduction of decay rate Γ/Γ0∼0.1 can be achieved, this rate apparently does not decrease with increasing N. Furthermore, the overlap fidelity between the numerically obtained eigenstate and the ansatz state does not approach 1, as the differing coherent interactions associated with q=0,±1 mix in additional contributions to the eigenstate. We note that, in contrast to the dispersion-engineered situation, these are not the most subradiant eigenstates in the Fz=0 manifold (as they are given by the phase-separated states). Thus, we have found these states by exactly diagonalizing the Hamiltonian and filtering out states such that ∑j⟨σ33j⟩>∑j⟨σ22j⟩, and then selecting the most subradiant states that satisfy that condition. In this situation, we need to diagonalize the full Hamiltonian matrix, instead of simply finding the eigenvalue with the smallest decay rate by using iterative sparse-matrix diagonalization algorithms, which limits the size of the chain that we can successfully diagonalize to N=12.