Theory of Subradiant States of a One-Dimensional Two-Level Atom Chain

TL;DR

This paper theoretically analyzes subradiant states in a one-dimensional two-level atom chain, revealing their fermionic nature and connecting multiply-excited states to the Tonks-Girardeau limit, while also identifying new correlated states.

Contribution

It provides a theoretical framework for singly- and multiply-excited subradiant states and links their behavior to the Tonks-Girardeau model, introducing new correlated states.

Findings

Subradiant states exhibit fermionic behavior.

Multiply-excited subradiant states relate to the Tonks-Girardeau limit.

Discovery of a new family of correlated states.

Abstract

Recently, the subradiant states of one-dimensional two-level atom chains coupled to light modes were found to have decay rates obeying a universal scaling, and an unexpected fermionic character of the multiply-excited subradiant states was discovered. In this Letter, we theoretically obtain the singly-excited subradiant states, and by eliminating the superradiant modes, we demonstrate a relation between the multiply-excited subradiant states and the Tonks-Girardeau limit of the Lieb-Liniger model which explains the fermionic behavior. In addition, we identify a new family of states with correlations different from the fermionic ansatz.

Equations

$$H_{\text{eff}}=-\frac{i}{2}{\Gamma }_{1D}\sum _{m,n=1}^N{e}^{i{k}_{1D}|z_m-z_n|}{\sigma }_m^{†}{\sigma }_n,$$

$$H_{\text{eff}}|k⟩={\omega }_k|k⟩-i\frac{{\Gamma }_{1D}}{2}(g_k|k_{1D}⟩-h_k|-k_{1D}⟩),$$

$${\delta }_{\xi }=\frac{\xi \pi }{Nd}\times \{\begin{array}{ll}1-i\frac{1}{N}\cot (\frac{k_{1D}}{2}d),& k\approx 0\\ 1+i\frac{1}{N}\tan (\frac{k_{1D}}{2}d),& k\approx -\pi /d\end{array}$$

$${\gamma }_{\xi }={\Gamma }_{1D}\frac{{\pi }^2}{2}\frac{{\xi }^2}{N^3}\times \{\begin{array}{ll}\frac{{\cos }^2(k_{1D}d/2)}{{\sin }^4(k_{1D}d/2)},& k\approx 0\\ \frac{{\sin }^2(k_{1D}d/2)}{{\cos }^4(k_{1D}d/2)},& k\approx -\pi /d.\end{array}$$

$$|{\varphi }_{k_{\xi }}⟩\propto$$

$$=$$

$$H_{3D,\mathrm{eff}}=-i\frac{3{\gamma }_0}{4k_0}{\int }_0^{k_0}\frac{d\tilde{k}}{2\pi }{\rho }_{+}(\tilde{k})\sum _{m,n=1}^N{e}^{i\tilde{k}|z_m-z_n|}{\sigma }_m^{†}{\sigma }_n-\frac{3{\gamma }_0}{4k_0}{\int }_0^{+\infty }\frac{d\tilde{k}}{2\pi }{\rho }_{-}(\tilde{k})\sum _{m,n=1}^N{e}^{-\tilde{k}|z_m-z_n|}{\sigma }_m^{†}{\sigma }_n,$$

$$Q=-\frac{{\Gamma }_{1D}}{8}\sum _{ϵ=\pm }\sum _{p,q}{b}_{ϵk_{1D}}^{†}{b}_{p+q-ϵk_{1D}}^{†}{b}_p{b}_q.$$

$$V_{\mathrm{sub}}$$

$$=\frac{1}{8}{\Gamma }_{1D}\sum _{m=1}^N(b_m^{†}{)}^2(b_m{)}^2.$$

$$\mathcal{H}=\sum_{i=1}^{n_{e}}\bigg{[}\frac{-\partial_{x_{i}}^{2}}{2m_{*}}+V(x_{i})\bigg{]}+2c_{LL}\sum_{i<j=1}^{n_{e}}\delta(x_{i}-x_{j}),$$

$${\omega }_{\xi }=\frac{{\Gamma }_{1D}}{2}\cot (\frac{k_{1D}}{2}d)+{\Gamma }_{1D}\frac{\cos (k_{1D}d/2)}{{\sin }^3(k_{1D}d/2)}(\frac{\xi \pi }{2N}{)}^2;$$

$${\omega }_{\xi }=-\frac{{\Gamma }_{1D}}{2}\tan (\frac{k_{1D}}{2}d)-{\Gamma }_{1D}\frac{\sin (k_{1D}d/2)}{{\cos }^3(k_{1D}d/2)}(\frac{\xi \pi }{2N}{)}^2.$$

$$H_{3D,\mathrm{eff}}=-{\mu }_0{\omega }_0^2\sum _{i,j=1}^N{\mathbf{d}}_i^{\ast }\cdot \mathbf{G}({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\cdot {\mathbf{d}}_j{\sigma }_i^{†}{\sigma }_j,$$

$$\displaystyle\mathbf{G}(\mathbf{r},\omega_{0})=\frac{e^{ik_{0}r}}{4\pi k_{0}^{2}r^{3}}\bigg{[}(k_{0}^{2}r^{2}$$

$$\displaystyle(-k_{0}^{2}r^{2}-3ik_{0}r+3)\frac{\mathbf{r}\mathbf{r}}{r^{3}}\bigg{]}$$

$$\sum _{i=1}^N\to \frac{1}{d}{\int }_0^{Nd}dx,b_i\to \sqrt{d}{b}_x.$$

$$H_{\text{eff}}^I|k_1,k_2⟩=\frac{N{\Gamma }_{1D}}{4}$$

$$-h_{k_2,ϵ}|b_{k_1,ϵk_{1D}}⟩+$$

$$g_{k,ϵ}$$

$$h_{k,ϵ}$$

$$H_{\mathrm{eff}}|k_1,k_2⟩=$$

$$-i\frac{N{\Gamma }_D}{2}$$

$$+c_{k_1,k_2,-}|-k_D,k_1+k_2+k_D⟩$$

$$\displaystyle\qquad+c_{k_{1},k_{2},+}|k_{1}+k_{2}-k_{D},k_{D}\rangle\bigg{)}.$$

Full text

To achieve controllable and deterministic photon-atom interfaces for applications in quantum information processing and quantum sensing, large atom ensembles may be used to enhance the coupling to photons Hammerer et al. (2010). The photons induce both coherent and dissipative atom-atom interactions that yield collective phenomena of super- or sub-radiance Van Loo et al. (2013), wherein a collective excitation of the atom ensemble decays faster or slower than individual atomic excitations. While superradiance has been extensively studied since the seminal work of Dicke Dicke (1954), subradiance of a large ensemble was observed only very recently in cold atom clouds Guerin et al. (2016); Weiss et al. (2018) and metamaterial arrays Jenkins et al. (2017). Comprehensive theoretical tools for the subradiance are still elusive Plankensteiner et al. (2015); Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019) due to the complicated long-range interactions and many-body features of the atomic ensembles Solano et al. (2017); Noh and Angelakis (2016); Olmos et al. (2013). A one-dimensional (1D) chain of equally spaced two-level atoms offers the simplest geometry to gain insight in the collective decay mechanisms, and implementation of such chains coupled to nanofibers Kornovan et al. (2016), 1D waveguides Haakh et al. (2016); Ruostekoski and Javanainen (2016, 2017); Zoubi (2014), and the full vacuum electromagnetic field in 3D free space Sutherland and Robicheaux (2016); Olmos et al. (2013); Bettles et al. (2016); Jen et al. (2016); Asenjo-Garcia et al. (2017) has attracted considerable attention. Super- and subradiance phenomena are in these systems supplemented by further interesting properties and applications such as atomic mirrors Chang et al. (2012), photon Fock state synthesis González-Tudela et al. (2017), enhancement of cooperativity Plankensteiner et al. (2017) and applications in quantum computation Paulisch et al. (2016).

Recently, the subradiant states of such 1D chains of $N$ qubits in 3D free space and coupled to 1D waveguide were numerically found to have a series of seemingly universal properties Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019): In the one-excitation sector where only one of the $N$ atoms is excited, if we sort all eigenstates (to be elaborated) by increasing decay rates with integer labels from $\xi=1$ to $\xi=N$, the most subradiant states ($\xi\ll N$) have decay rates $\gamma_{\xi}\propto\xi^{2}/N^{3}$. In the multi-excitation sectors, the most subradiant states have a fermionic character, e.g., a most subradiant state with two excitations is given by $|F_{1,2}\rangle\propto\sum_{i<j}(c_{1,i}c_{2,j}-c_{1,j}c_{2,i})|e_{i},e_{j}\rangle$, built from subradiant states $|\psi_{1(2)}\rangle=\sum_{i}c_{1(2),i}|e_{i}\rangle$ in the one-excitation sector, where $|e_{i}\rangle$ ($|e_{i},e_{j}\rangle$) represents the state with the $i$th ($i$th and the $j$th) atom excited to $|e\rangle$ while all other atoms are in the ground state $|g\rangle$. The decay rate of $|F_{1,2}\rangle$ is the sum of the decay rates of $|\psi_{1}\rangle$ and $|\psi_{2}\rangle$. The infidelity of the fermionic ansatz $|F_{1,2}\rangle$ to exact numerical results scale as $N^{-2}$ and $N^{-1}$ for two different classes of states Albrecht et al. (2019).

It is intriguing why these properties appear for both the infinite range atom-atom interactions mediated by 1D guided fields Albrecht et al. (2019); Henriet et al. (2019), and the long-range ($\sim 1/r$) or short-range ($\sim 1/r^{3}$) interactions mediated by the 3D free-space field Asenjo-Garcia et al. (2017); Henriet et al. (2019). A thorough theoretical understanding is needed to guide further experimental studies and applications of subradiance. In this Letter, we provide such understanding based on the theoretical treatment of the physics summarized in Fig. 1.

Spin Models- For light-matter interactions where the Markovian approximation is applicable, we can eliminate the light modes to obtain a master equation describing only the atoms Dung et al. (2002). The master equation is equivalent to a Monte Carlo wave function formalism Mølmer et al. (1993), where the atomic state evolves stochastically under quantum jumps and deterministically under $H=H_{0}+H_{\text{eff}}$, where $H_{0}$ is the bare Hamiltonian of the atoms and the non-Hermitian $H_{\text{eff}}$ describes both coherent and dissipative atom-atom interactions mediated by the vacuum field. The right eigenstates of $H$, or equivalently of $H_{\text{eff}}$ Asenjo-Garcia et al. (2017), in each manifold of states with any given number of atomic excitations have decay rates that are twice the imaginary part of the corresponding right eigenvalues. We focus our analysis on the qubit chain coupled to a 1D waveguide, but our treatment provides sufficient insight to also account for the case of coupling to the 3D vacuum field. For an atom chain coupled to a 1D waveguide, we have Chang et al. (2012); Albrecht et al. (2019):

[TABLE]

where $\Gamma_{1D}$ is the decay rate of a single atom coupled to the waveguide Chang et al. (2012), $k_{1D}$ is the wavenumber of the waveguide mode resonant with the atomic transition, and $\sigma_{m}=|g\rangle_{m}\langle e|$ acts on the $m$th atom. We assume the atoms are equidistantly spaced by $d$. For convenience, we shall denote $H_{\mathrm{eff}}=H^{R}_{\text{eff}}-iH^{I}_{\text{eff}}$.

One-Excitation Sector-The eigenstates of $H^{I}_{\text{eff}}$ divide the one-excitation sector into, a two-dimensional superradiant subspace spanned by Bloch states $|{\pm}k_{1D}\rangle=N^{-1/2}\sum_{m=1}^{N}e^{\pm ik_{1D}z_{m}}|e_{m}\rangle$ with eigenvalue $N\Gamma_{1D}/4$; and an $(N{-}2)$-dimensional dark space with eigenvalue 0, see Fig. 1. The dark states acquire weak (subradiant) decay rates because of their admixture of superradiant states induced by the perturbation from $H^{R}_{\text{eff}}$.

While the perturbation view is informative, a more direct approach to the subradiant states applies the following exact result for the Bloch states $|k\rangle$ ($k\neq\pm k_{D}$),

[TABLE]

where $\omega_{k}{=}\frac{\Gamma_{1D}}{4}\sum_{\epsilon=\pm}\cot(\frac{k_{1D}+\epsilon k}{2}d)$, and the “tails” $g_{k}=\frac{e^{i(k-k_{1D})z_{1}}}{1-e^{i(k-k_{1D})d}}$ and $h_{k}=\frac{e^{i(k+k_{1D})z_{N}}}{e^{-i(k+k_{1D})d}-1}$. It follows that a superposition of two degenerate states, $|k\rangle$ and $|{-}k\rangle$, is an eigenstate of $H_{\mathrm{eff}}$ with eigenvalue $\omega_{k}$ and has no tails if $k$ is a solution to the equation $g_{k}h_{-k}=g_{-k}h_{k}$. This equation has only solutions for complex values of $k$. In the regimes $k\approx 0$ or $\pm\pi/d$ (center or edges of the first Brillouin zone), supposing $k_{\xi}=0+\delta_{\xi}$ and $k_{\xi}=-\pi/d+\delta_{\xi}$ respectively, we find to order $N^{-2}$,

[TABLE]

with $\xi=1,2,3\cdots$, $\xi\ll N$. Note that Eq. (3) amounts to an $1/N^{2}$-order imaginary correction to the Bloch wavenumber.

Next, we substitute Eq. (3) into the expression for $\omega_{k}$, which is parabolic near $k\approx 0$ and $\pm\pi/d$, i.e., $\omega_{k}\propto\delta_{\xi}^{2}$. Then the imaginary corrections directly yield the $\xi^{2}/N^{3}$-scaling of the decay rates Albrecht et al. (2019):

[TABLE]

The eigenstates are written as

[TABLE]

where $k_{\xi}^{(0)}=\xi\pi/(Nd)$ or $-\pi/d+\xi\pi/(Nd)$.

Universality-The $\xi^{2}/N^{3}$-scaling has also been numerically found for 1D atom chains coupled to 3D free-space modes Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019), where the effective Hamiltonian $H_{3D,\mathrm{eff}}$ is determined by the Green’s dyadic tensor (see the Supplemental Material). Fourier transformation of the Green’s tensor reveals a hidden similarity between the coupling to the 1D and 3D quantized radiation fields: $H_{3D,\mathrm{eff}}$ can be written as weighted integrals of terms resembling $H_{\mathrm{eff}}$ with real-valued $k_{1D}\in[0,\,k_{0}]$ and imaginary-valued $k_{1D}\in[i0,+i\infty]$:

[TABLE]

where $\gamma_{0}$ is the spontaneous emission rate and $k_{0}$ is the resonant wavenumber. If the atoms are polarized parallel to the chain, $\rho_{\pm}(\tilde{k})=2\pi(1\mp\tilde{k}^{2}/k_{0}^{2})$ and the atom-atom interaction is short-range ($\sim 1/r^{3}$). If the atoms are polarized transverse to the chain, $\rho_{\pm}(\tilde{k})=\pi(1\pm\tilde{k}^{2}/k_{0}^{2})$ and the atom-atom interaction is long-range ($\sim 1/r$).

In combination with the two exact features of our analytical results for $H_{\mathrm{eff}}$:

The leading order solutions of $\delta_{\xi}$ and $|\phi_{k_{\xi}}\rangle$ are independent of the values of $k_{1D}$,see Eq. (5); 2. 2.

The proportionality $\delta_{\xi}\propto\xi$ and the parabolic dispersion relation $\omega_{k}\propto\delta_{\xi}^{2}$ hold to order-$N^{{-}2}$, for all values of $k_{1D}$,

this explains the universality of the $\xi^{2}/N^{3}$-scaling: Feature 1 implies that the leading order eigenstates of $H_{\mathrm{eff}}$, shown in Eq. (5), are shared simultaneously by all terms integrated in $H_{3D,\mathrm{eff}}$, and thus by the full $H_{3D,\mathrm{eff}}$ due to linearity. Feature 2, hence implies that the corresponding decay rates, scaling as $\xi^{2}/N^{3}$, also apply to the subradiant states of $H_{3D,\mathrm{eff}}$. The prerequisite is that $H^{I}_{3D,\mathrm{eff}}$ must have dark states with $k\approx\pm\pi/d$. Hence we require $k_{0}<\pi/d$ which implies that the ensemble is only sub-radiant in the 3D field if the atom-atom distance is less than half the resonant wave length Olmos et al. (2013); Asenjo-Garcia et al. (2017).

Subradiant multiply-excited states- When the number of atomic excitations $n_{e}\ll N$, the leading order Holstein-Primakoff (HP) approximation Holstein and Primakoff (1940) usually applies and one may replace $\sigma^{\dagger}_{m}\sigma_{n}$ of Eq. (1) with the bosonic operators $b_{m}^{\dagger}b_{n}$ and obtain a quadratic bosonic $H_{\mathrm{eff}}$. It works well for the superradiant modes with wavenumber $\pm k_{1D}$. But for the subradiant multiply-excited states, the bosonic creation operators prepare exchange symmetric combinations of subradiant one-excitation states with decay rates scaling as $N^{-1}$ Asenjo-Garcia et al. (2017); Albrecht et al. (2019) which is much larger than the numerically observed $N^{-3}$-scaling Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019). Instead, the numerical results were found to favor fermionic exchange anti-symmetric combinations of the subradiant one-excitation states Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019).

This somewhat surprising result inspires a closer scrutiny of the HP transformation. Including second order corrections due to saturation, the HP transformation reads $\sigma_{m}=(1-b_{m}^{\dagger}b_{m}/2)b_{m}$ so that we can write $H^{I}_{\text{eff}}=H_{\mathrm{SR}}+Q+Q^{\dagger}$ with the quadratic term $H_{\mathrm{SR}}=N\Gamma_{1D}/4\sum_{\epsilon=\pm}b_{\epsilon k_{1D}}^{\dagger}b_{\epsilon k_{1D}}$ and quartic terms

[TABLE]

Here, $b^{\dagger}_{k}=N^{-1/2}\sum_{m}e^{ikz_{m}}b^{\dagger}_{m}$ and the summation over wavenumber is taken over an orthonormal basis $\{|k\rangle\}_{k}$ containing $|{\pm}k_{1D}\rangle$. The quadratic term $H_{SR}$ has a prefactor $N$-times larger than those of $Q$ and $Q^{\dagger}$, but to assess their influence, we should take account of not only the prefactors but also the magnitudes of the operator terms. For $H_{SR}$, the magnitude of $b_{\epsilon k_{1D}}^{\dagger}b_{\epsilon k_{1D}}$ can be estimated by its typical expectation value over the relevant Hilbert space, i.e., the subradiant states.

Reasonably, one may expect that a subradiant state contains no excitation of superradiant modes, i.e., typically $\langle b_{\epsilon k_{1D}}^{\dagger}b_{\epsilon k_{1D}}\rangle\approx 0$. Thus the magnitude of $H_{SR}$ is suppressed. Meanwhile, Eq. (7) shows that $Q$ annihilates a two-boson dark state, $b^{\dagger}_{p}b^{\dagger}_{q}|\varnothing\rangle$, with respect to $H_{SR}$, and generates a superradiant two-boson state $b^{\dagger}_{\epsilon k_{1D}}b^{\dagger}_{p+q-\epsilon k_{1D}}|\varnothing\rangle$ ($|\varnothing\rangle$ denotes the boson vacuum). This demonstrates that the saturation correction to the HP approximation plays a significant role even in the low excitation regime ($n_{e}\ll N$), in contrast to its role in many other applications.

An effective theory for how $Q$ couples the dark states to superradiant states and hereby determines their subradiant behavior is illustrated in the right panel of Fig. 1. The effect of $Q$ and $Q^{\dagger}$ is distilled by eliminating the superradiant states, in a manner similar to the adiabatic elimination of excited state manifolds to restrict the effective dynamics of quantum systems to their ground state manifold Reiter and Sørensen (2012).

Note that the subset of superradiant states with only a single excitation of the superradiant modes and thus the eigenvalue (decay rate) $N\Gamma_{1D}/4$, has the strongest coupling to the dark/subradiant states. We hence disregard the coupling to other superradiant states and the effective coupling among subradiant states reduces to $V_{\mathrm{sub}}=\frac{4}{N\Gamma_{1D}}P_{DS}Q^{\dagger}P_{SRS}QP_{DS}$, with projection operators $P_{DS(SRS)}$ on the dark and superradiant spaces, respectively. To evaluate this expression we use the operator relation that $(b_{p^{\prime}+q^{\prime}-\epsilon^{\prime}k_{1D}}b_{\epsilon^{\prime}k_{1D}})(b^{\dagger}_{\epsilon k_{1D}}b^{\dagger}_{p+q-\epsilon k_{1D}})=\delta_{\epsilon,\epsilon^{\prime}}\delta_{p^{\prime}+q^{\prime},p+q}$, i.e., no population of the superradiant two-boson modes within the dark/subradiant space. Finally, we obtain

[TABLE]

That is, $V_{\mathrm{sub}}$ induces decay with rate $O(\Gamma_{1D})$ of nominally subradiant states having more than a single HP boson excitation at the same site.

In the absence of $V_{\mathrm{sub}}$, approximate multiply-excited states are created by the operators $b^{\dagger}_{\xi}=\sum_{m}\langle e_{m}|\phi_{k_{\xi}}\rangle b^{\dagger}_{m}$, with $|\phi_{k_{\xi}}\rangle$ the one-excitation eigenstates Eq. (5) of $H_{\mathrm{eff}}$. As $V_{sub}$ cannot be treated as a perturbation, we study the effective Hamiltonian $\mathcal{H}=\frac{1}{2}\sum_{\xi}\gamma_{\xi}b_{\xi}^{\dagger}b_{\xi}+V_{\mathrm{sub}}$, where only the most subradiant states ($\xi\ll N$) are included in the sum. In the Supplemental Material we show that in the continuous limit, $\mathcal{H}$ can be written as the second-quantized form of the Hamiltonian

[TABLE]

where $c_{LL}=d\Gamma_{1D}/8$ and $V(x_{i})$ is a 1D box potential for bosons in the interval $[0,Nd]$ with one-excitation eigenstates $|\phi_{k_{\xi}}\rangle$ given by Eq. (5). The observation behind Eq. (9) is that the $\xi^{2}/N^{3}$-scaling of $\gamma_{\xi}$ takes the same form of a kinetic energy $\gamma_{\xi}=k_{\xi}^{2}/(2m_{*})$ when $k\approx 0$; or the kinetic energy in a gauge field $\gamma_{\xi}=(k_{\xi}+\pi/d)^{2}/(2m_{*})$ when $k\approx-\pi/d$. With the parametrization of the model, the effective mass in the kinetic energy term reads $m_{*}=\xi^{2}\pi^{2}/(N^{2}d^{2}\gamma_{\xi})\propto N$.

We recognize Equation (9) as the Lieb-Liniger model Lieb and Liniger (1963) originally proposed for 1D gases of hard-core bosons. As the effective mass $m_{*}$ diverges in the large $N$ limit, the kinetic energy-like part of Eq. (9) vanishes. This implies that Eq.(9) reaches the Tonks-Girardeau regime Tonks (1936); Girardeau (1960) of the Lieb-Liniger model, where the eigenstates of $\mathcal{H}$ can be obtained via a fermion-boson mapping Girardeau (1960); Cazalilla et al. (2011): For a free fermion model described by $\sum_{\xi}\frac{1}{2}\gamma_{\xi}f_{\xi}^{\dagger}f_{\xi}$, where $f^{\dagger}_{\xi}=\sum_{m}\langle e_{m}|\phi_{\xi}\rangle f^{\dagger}_{m}$, we write down its eigenstates (e.g., two-fermion states) $f^{\dagger}_{\xi_{1}}f^{\dagger}_{\xi_{2}}|\varnothing\rangle$, and replace $f^{\dagger}_{m}f^{\dagger}_{n}$ with $\mathrm{sign}(n-m)b^{\dagger}_{m}b^{\dagger}_{n}$, where $\mathrm{sign}(n-m)$ is necessary to ensure the consistency with the fermionic commutation relation. This yields a fermion-like bosonic state $|F_{\xi_{1},\xi_{2}}\rangle=\sum_{m<n}[\phi_{k_{\xi_{1}}}(z_{m})\phi_{k_{\xi_{2}}}(z_{n})-e^{i\varphi}\phi_{k_{\xi_{2}}}(z_{m})\phi_{k_{\xi_{1}}}(z_{n})]b^{\dagger}_{m}b^{\dagger}_{n}|\varnothing\rangle$ ($e^{i\varphi}=1$ is introduced for later convenience). Since $b^{\dagger}_{m}b^{\dagger}_{n}|\varnothing\rangle=|e_{m},e_{n}\rangle$, we recover the fermionic ansatz of the two-excitation sector Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019). As a direct consequence of their representation as non-interacting fermions the decay rates of the most subradiant multiply-excited states are merely the sum of the decay rates of their one-excitation constituents, e.g., $\gamma_{\xi_{1}}+\gamma_{\xi_{2}}$. This explains the numerical observations of Refs. Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019). The above mapping also applies to more excitations Girardeau (1960); Cazalilla et al. (2011).

For a finite atom chain, $m_{*}$ is finite and Eq. (9) deviates from the Tonks-Girardeau limit, hence the fermionic ansatz, for example the two-excitation state $|F_{\xi_{1},\xi_{2}}\rangle$, deviates slightly from the numerical eigenstates denoted by $|\psi^{\mathrm{num}}_{\xi_{1},\xi_{2}}\rangle$. The deviation quantified by the infidelity $1-|\langle F_{\xi_{1},\xi_{2}}|\psi^{\mathrm{num}}_{\xi_{1},\xi_{2}}\rangle|^{2}$, is numerically found to scale as $N^{-2}$ when both components, $|\phi_{k_{\xi_{1}}}\rangle$ and $|\phi_{k_{\xi_{2}}}\rangle$, come from the same branch of the one-excitation subradiant states, i.e., $k_{\xi_{1}},\,k_{\xi_{2}}\approx\pm\pi/d$ (or both $\approx 0$); otherwise, the infidelity scales as $N^{-1}$ (when $k_{\xi_{1}}\approx\pm\pi/d$ and $k_{\xi_{2}}\approx 0$) Albrecht et al. (2019). These behaviors can also be explained from the Lieb-Liniger model of Eq. (9). The fermion-boson mapping is not exact and the phase factor $e^{i\varphi}$ introduced above deviates from unity by a factor in the form of $(k_{\xi_{1}}-k_{\xi_{2}})/(m_{*}c_{LL})$ Lieb and Liniger (1963). Since $k_{\xi_{1}}-k_{\xi_{2}}$ is $O(N^{-1})$ or $\approx\pi/d$ in the two cases considered, while $m_{*}c_{LL}$ scales as $N$, their ratio scales exactly in the same manner as the numerically observed infidelities Albrecht et al. (2019). Larger discrepancies with the fermionic ansatz are detectable when the decay rates increase.

*Universality-*The mapping to Lieb-Liniger model and the Tonks-Girardeau gas can also be extended to the 1D atomic chain coupled to 3D free-space modes described by Eq. (6). Here $H^{I}_{3D,\mathrm{eff}}$ possesses short-lived eigenstates $|\tilde{k}\rangle$ with $\tilde{k}\in[{-}k_{0},k_{0}]$ and different decay rates $\gamma_{\tilde{k}}$. Each of them will contribute to $V_{\mathrm{sub}}$ a term with prefactor $\gamma_{\tilde{k}}/N^{2}$. Hence we have $V_{\mathrm{sub}}\propto\sum_{\tilde{k}}\gamma_{\tilde{k}}/N^{2}=\gamma_{0}/N$, similar to the coefficient in the first line of Eq. (8). Since the $\xi^{2}/N^{3}$-scaling decay rates apply in the one-excitation sector in 3D free-space, the fermionic ansatz also applies here. Since only subradiant states with $k\approx\pm\pi/d$ appear in 3D free-space, the pertaining $N^{-2}$ scaling applies to the infidelities of all states given by the fermionic ansatz. This matches the numerical results Asenjo-Garcia et al. (2017).

Conclusion and Discussion- In this Letter, we have developed a theory to explain the $\xi^{2}/N^{3}$-scaling of subradiant decay rates and the fermionic behavior of multiply-excited subradiant states identified in numerical calculations on 1D atom chains coupled to both 1D and 3D radiation reservoirs Asenjo-Garcia et al. (2017); Albrecht et al. (2019); Henriet et al. (2019). We find that the universal $\xi^{2}/N^{3}$-scaling results from a parabolic dispersion relation of the atomic excitation, and imaginary corrections of the Bloch quasi-momentum eigenstates of the non-Hermitian Hamiltonian. For multiply excited systems quartic corrections to the Holstein-Primakoff (HP) expansion of the effective spin Hamiltonian for the atom chain dominate the coupling of the sub-radiant states, and lead to a formulation equivalent to the Lieb-Liniger model of a 1D bosonic quantum gas in the Tonks-Girardeau regime Tonks (1936); Girardeau (1960). The fermionic ansatz solution of that problem explains the decay rates and the properties of the solutions found in Ref. Asenjo-Garcia et al. (2017); Albrecht et al. (2019). There is a high current interest and many potential applications of subradiance Chang et al. (2012); González-Tudela et al. (2017); Plankensteiner et al. (2017); Paulisch et al. (2016) and the analytical findings presented here may inspire further study of subradiance in light-matter interactions of more complex geometries, e.g., systems with higher dimensional atom arrays Facchinetti et al. (2016), chiral waveguides that break the parity symmetry Ramos et al. (2014) and setups with topological effects Perczel et al. (2017); Ozawa et al. .

Let us conclude by discussing a remaining theoretical issue. We recall our effective separation of the Hamiltonian into an interaction term, $V_{\mathrm{sub}}$, based on $H^{I}_{\mathrm{eff}}$ and an expansion on subradiant eigenmodes $b_{\xi}$ for which $H^{R}_{\mathrm{eff}}$ contributes the decay rates $\gamma_{\xi}$. Like the numerical calculations, a more rigorous analytical approach should incorporate $H^{R}_{\mathrm{eff}}$ and $H^{I}_{\mathrm{eff}}$ on an equal footing. The fact that our separate treatment applies may be understood from the perturbation view. The leading order approximation of the subradiant states are the dark states of $H^{I}_{\mathrm{eff}}$. They are also approximate eigenstates of $H^{R}_{\mathrm{eff}}$ when restricted to the most subradiant states. It means that the fermionic ansatz, as the leading order approximation, is shared by both $H^{I}_{\mathrm{eff}}$ and $H_{\mathrm{eff}}$. Therefore analyzing $V_{\mathrm{sub}}$ from the simpler $H^{I}_{\mathrm{eff}}$ is sufficient to capture the salient fermionic behavior. This is also verified by a direct construction of the fermionic ansatz without using the HP transformation, for both $H^{I}_{\mathrm{eff}}$ and $H_{\mathrm{eff}}$ (see the Supplemental Material). Interestingly, we find that the fermionic ansatz does not exhaust all the most subradiant eigenstates. For a medium-size ensemble of $N=20$ atoms, we obtained numerical eigenstates of $H_{\mathrm{eff}}$ with very low fermionic state fidelity. The subradiant states of this different character have well defined “center of mass” wavenumber, and well defined spatial separation, as illustrated in Fig. 2(b). Further discussion of these states is beyond the scope of this manuscript, but may be of interest for future work possibly together with the interesting prospects for studying quantum fluctuations Paredes et al. (2004); Jacqmin et al. (2011); Budde and Mølmer (2004) in the Tonk-Girardeau gas theory by detection of the excited state correlations among atoms in a subradiant chain.

Acknowledgments-This work was supported by the Villum Foundation and by the European Unions Horizon 2020 research and innovation program (Grant No. 712721, NanOQTech).

I Supplemental Material

I.1 A. Energy Levels of the Subradiant States

We derived the decay rates of the subradiant states in the one-excitation sector. As a byproduct, their energy levels are given as following: for the subradiant states with $k_{\xi}\approx\xi\pi/(Nd)\;(\xi\ll N)$, we have

[TABLE]

for those with $k_{\xi}\approx-\pi/d+\xi\pi/(Nd)$ ($\xi\ll N$), we have

[TABLE]

Apart from the constant part, the above expressions show that the subradiant states have Lamb shifts proportional to $\xi^{2}/N^{2}$. The band is parabolic and flat around the extreme point $k=0$ or $k=\pm\pi/d$.

I.2 B. Hamiltonian of 1D Atom Chain in 3D Free-Space

The effective atom-atom coupling Hamiltonian is expressed as

[TABLE]

where $\mathbf{d}_{i}$ and $\mathbf{r}_{i}$ is the dipole moment and the position of the ith atom, respectively, and $\mu_{0}$ is the vacuum permeability. For the case of 3D free-space, the dyadic Green’s tensor $\mathbf{G}(\mathbf{r}_{i},\mathbf{r}_{j},\omega)=\mathbf{G}(\mathbf{r}_{i}-\mathbf{r}_{j},\omega)$ is given as

[TABLE]

where $k_{0}=\omega_{0}/c$. This expression can be transformed to the wave number presentation and yields $H_{3D,\mathrm{eff}}$ presented as Eq. (6) of the main text.

I.3 C. Transformation to Continuous Limit

Equation (9) of the main text is written in a discrete notation. The continuous expression can be obtained from the discrete notation by the mapping

[TABLE]

The bosonic commutation relation changes from $[b_{i},b_{j}^{\dagger}]=\delta_{i,j}$ to $[b_{x},b_{y}^{\dagger}]=\delta(x-y)$.

I.4 D. Direct construction

To see why the fermionic ansatz has the $N^{-3}$ decay rate, we introduce the two-excitation state $|k_{1},k_{2}\rangle=\sum_{m<n}e^{ik_{1}z_{m}+ik_{2}z_{n}}|e_{m},e_{n}\rangle$, and evaluate

[TABLE]

where

[TABLE]

and $c_{k_{1},k_{2},\epsilon}{=}g_{k_{1},\epsilon}^{*}+g_{k_{2},\epsilon}$ and $|b_{k,k^{\prime}}\rangle{=}|k,k^{\prime}\rangle+|k^{\prime},k\rangle$. All these state amplitudes scale as $N^{-1}$ rather than the desired $N^{-2}$-scaling, which is required to obtain the $1/N^{3}$-scaling of the decay rates.

As in the one-excitation sector, we may proceed and construct superpositions of four degenerate states $|k_{1},k_{2}\rangle$, $|{-}k_{1},k_{2}\rangle$, $|k_{1},{-}k_{2}\rangle$ and $|{-}k_{1},{-}k_{2}\rangle$, to $|\phi_{k_{1}},\phi_{k_{2}}\rangle=\sum_{m<n}\phi_{k_{1}}(z_{m})\phi_{k_{2}}(z_{n})|e_{m},e_{n}\rangle$, where $\phi_{k}(z_{m})=\langle e_{m}|\phi_{k}\rangle$. Then in the expression of $H^{I}_{\text{eff}}|\phi_{k_{1}},\phi_{k_{2}}\rangle$, the state amplitudes on $|b_{k_{1(2),\epsilon k_{1D}}}\rangle$, but not on $|b_{k_{1}+k_{2}-\epsilon k_{1D},\epsilon k_{1D}}\rangle$, are successfully reduced to the $N^{-2}$-scaling.

To also suppress the latter, we form the superposition with the permuted state $|\phi_{k_{2}},\phi_{k_{1}}\rangle$. The suitable superposition turns out to be “fermionic”, i.e., $|F_{k_{1},k_{2}}\rangle\propto|\phi_{k_{1}},\phi_{k_{2}}\rangle-|\phi_{k_{2}},\phi_{k_{1}}\rangle$. Different from what we have seen in the one-excitation sector, here the “tails” cannot be erased completely by the superposition. It means that the Fermionic ansatz is only the leading order solution.

With reference to our concluding remarks in the main text on the full treatment of $H_{\mathrm{eff}}^{R}$ and $H_{\mathrm{eff}}^{I}$, the same construction can be applied for the full Hamiltonian $H_{\mathrm{eff}}$, where we obtain

[TABLE]

While Eq. (17) features “tails” that are neither symmetric nor anti-symmetric, its main features, and hence the applicability of the fermionic ansatz, are captured by $H^{I}_{\mathrm{eff}}$ given in Eq. (15). As we show in the main text, the “tails” $|b_{k_{1}+k_{2}-\epsilon k_{1D},\epsilon k_{1D}}\rangle$ of Eq. (15) arise via the second order terms in the Holstein-Primakoff (HP) transformation, which in turn establish the Tonks-Girardeau limit of the Lieb-Liniger model, and hence the fermionic solutions.