Effects of finite trapping on the decay, recoil, and decoherence of dark states of quantum emitter arrays

We consider $N$ atoms of equal mass $m$. Each atom $i$ has two internal states, $\ket{g_{i}}$ and $\ket{e_{i}}$, with a transition frequency $\omega_{0}$ and a spontaneous decay rate $\gamma_{0}$. The operators $\sigma_{i}^{+}$ and $\sigma_{i}^{-}$ are the raising and lowering operators of the internal states of atom $i$: $\sigma_{i}^{+}=\ket{e_{i}}\bra{g_{i}}$ and $\sigma_{i}^{-}=\ket{g_{i}}\bra{e_{i}}$.

Each atom is trapped in a harmonic well. The location of the center of the trap of atom $i$ is denoted by $\textbf{R}_{i}$, while the location of the atom $i$ relative to its trap center is $\textbf{r}_{i}$. For the relative position of two traps, we use the notation $\textbf{R}_{ij}=\textbf{R}_{i}-\textbf{R}_{j}$ and the relative displacement as $\textbf{r}_{ij}=\textbf{r}_{i}-\textbf{r}_{j}$. We also use the absolute position of an atom which is denoted by $\tilde{\textbf{r}}_{i}=\textbf{R}_{i}+\textbf{r}_{i}$, with corresponding difference of two atoms $\tilde{\textbf{r}}_{ij}=\tilde{\textbf{r}}_{i}-\tilde{\textbf{r}}_{j}$. A schematic of the system is shown in Fig. 1. For simplicity, the traps are of equal angular frequency $\omega_{t}$ and evenly separated by a distance $d=|\textbf{R}_{i+1}-\textbf{R}_{i}|$. The eigenstates of the trap for atom $i$ in each direction are denoted by $\ket{n_{i}}$ with energy $n_{i}\hbar\omega_{t}$ (where the conventional 1/2 is dropped for convenience).

The density matrix of the system, $\rho$, can be expanded in the tensor product basis of the vibrational and internal states for each atom. We denote such a generic basis state by $\ket{A,V}$ where $A$ and $V$ are the collection of the internal and vibrational quantum numbers of the system, respectively: $\ket{A}=\ket{\alpha_{1},\alpha_{2},\dots}$, and $\ket{V}=\ket{n_{1},n_{2},\dots}$, where at most one $\alpha_{j}$ can be $e_{j}$ while the others in $g_{j}$. In this way, there are $N+1$ allowed internal states. Each $n_{j}$ can be $0,1,2,\dots$. Hence, the density matrix can be written as

For numerical feasibility, we can truncate the vibrational states to $0,1,2,\dots,n_{max}$, where we choose a phonon cut-off, $n_{max}$, to insure convergence of the numerical simulation. This way, we have $N_{vib}=n_{max}+1$ vibrational states per atom. The vibrational Hilbert space will, therefore, contain $(N_{vib})^{N}$ states. The total number of states in the system is $N_{tot}=(N_{vib})^{N}(N+1)$.

The evolution of the system’s density matrix, $\rho$, is governed by a master equation

where $H_{\text{eff}}$ is an effective Hamiltonian that contains a trap term and a non-Hermitian interaction term, and $\mathcal{R}$ is a population recycling superoperator Rubies-Bigorda et al. (2024a) defined below. The trap term is given by

where $a_{i}$ is the lowering operator of the vibrational states of the trap of atom $i$. The interaction term comes from the coupling with the electromagnetic field (which is traced out) and is given by

where $g$ is a complex valued Green’s function proportional to the electric field propagator and is given in Eqs. (7) and (8) below. The imaginary part of this Green’s function gives a non-Hermitian interaction responsible for the collective decay of the excited state. On the other hand, the real part gives a Hermitian interaction responsible for the collective Lamb shifts of the energy of the states and also the induced forces which can shift the positions of the atoms before the decay of the excited state Scully (2009).

The final term in Eq. (2) is the recycling superoperator, $\mathcal{R}$, which accounts for the decayed population and the recoil energy during the photon emission step. The recycling superoperator acts on the density matrix in a relatively complicated way. The simplest description is in the position representation of the density matrix, where $\rho$ is expanded as

where $r$ and $r^{\prime}$ are the collections of the positions of the atoms. The recycling superoperator is given by

where $\text{Im}[g]$ is the imaginary part of the Green’s function, and $\tilde{\textbf{r}}^{\prime}_{ij}=\tilde{\textbf{r}}_{i}-\tilde{\textbf{r}}^{\prime}_{j}$. The multiplication of the density matrix by the position dependent Green’s function in the above equation is responsible for and explains the mechanism of the imparting of momentum and energy to the atoms during the photon emission step. We use a spectral method employing the eigenstates of a truncated position operator $(a+a^{\dagger})$ as the basis to calculate this superoperator action. This results in a faster convergence with the number of basis states as done in the appendix of Ref. Robicheaux (2025).

The Green’s function in Eq. (4) depends on the underlying geometry of the system. For atoms coupled to a 1D waveguide (assuming no other decay channel), the Green’s function has a simple form and leads to coupling given by Asenjo-Garcia et al. (2017b):

where $k_{0}=2\pi/\lambda_{0}=\omega_{0}/c$ is the wave number of light, $\lambda_{0}$ is the wavelength of light, and $c$ is the speed of light. For a waveguide along the x-axis, $\textbf{r}=x\hat{x}$ and $|\textbf{r}|=|x|$ denotes the distance between two atoms in 1D. This function is simple because it is bounded and periodic.

On the other hand, the Green’s function for free space is more complicated and given in Eq. (8):

where $\hat{q}$ is the dipole orientation, $r=|\mathbf{r}|$ is the norm of the vector $\mathbf{r}$, $\hat{r}=\mathbf{r}/r$ is the unit vector along $\mathbf{r}$, and $h_{l}^{(1)}(x)$ are the outgoing spherical Hankel functions of angular momentum $l$; $h_{0}^{(1)}(x)=e^{ix}/[ix]$ and $h_{2}^{(1)}(x)=(-3i/x^{3}-3/x^{2}+i/x)e^{ix}$. Unlike the 1D case, this Green’s function is non-periodic. Its magnitude decays like $1/r$ at large $r$, and its real part diverges like $1/r^{3}$ at small $r$. In Eq. (4), this leads to infinities in the self-interaction terms $i=j$, and we remove the real part of the Green’s function at $i=j$ to avoid the infinities.

The distinction between $\textbf{R}_{i}$ and $\textbf{r}_{i}$ (see Fig. 1) is important for the description of the motional aspect of the system but also has conceptual consequences for the following reason. The atoms interact through the electromagnetic field. As described in Eq. (4) and in Eqs. (7 and 8), the interaction strength depends on the phase difference of light due to the difference of absolute positions: $k_{0}|\tilde{\textbf{r}}_{ij}|$, where $k_{0}$ is the wave number of light. This has contribution from the relative positions of the traps, $k_{0}\textbf{R}_{ij}$, which is a fixed value. The other contribution comes from $k_{0}\textbf{r}_{ij}$, which is a quantum operator that describes the spatial degree of freedom of an atom inside its trap.

In the theoretical limit of infinitely trapped or infinitely massive atoms, this operator can be taken to be taken to be fixed, $k_{0}\textbf{r}_{ij}=0$. This approximation is often used in theoretical studies of collective decay. For finite trapping strength, the position operator

where $\eta$ is the Lamb-Dicke (LD) parameter given by $\eta=k_{0}\sqrt{\frac{\hbar}{2m\omega_{t}}}$. This LD parameter serves two purposes. First, it quantifies the coupling between the internal and vibrational states. The coupling $g(\tilde{\textbf{r}}_{ij})$ in Eq. (4) can be Taylor expanded in the LD parameter to give terms where the vibrational operators, $a_{i}$ and $a_{i}^{\dagger}$, are multiplied by the internal operators $\sigma_{i}^{+}$ and $\sigma_{i}^{-}$. Second, the LD parameter quantifies the spatial spread of vibrational state of an individual atom. For example, the standard deviation of position divided by $\lambda_{0}$ is equal to ${\eta}/{2\pi}$ for the ground vibrational state. Typically, the LD regime is characterized by a spread that is much smaller than the wavelength, corresponding to $\eta\ll 2\pi$.

The master equation in Eq. (2) is computationally expensive to solve. In the case of starting in a pure state of the single excitation subspace with no driving laser, the master equation can be solved more efficiently in two parts: 1) a Schrödinger equation for the electronic excited state that does not depend on the electronic ground state dynamics. 2) a recycling superoperator that accumulates the decayed population in the electronic ground state.

In this paper, we always start the system from a pure state $\ket{\psi(t=0)}=\ket{DS,0}$, where $\ket{DS}$ is an electronic eignestate and $\ket{0}$ is the vibrational ground state of all atoms. In this case, the Schrödinger equation is a pure state equation:

The magnitude of this state, $\innerproduct{\psi(t)}{\psi(t)}$ is equal to the probability of being in the excited state, $p(t)$. This probability decays over time due to the non-hermiticity of the effective Hamiltonian. The rate of change of probability is given by

with a corresponding decay rate defined as $\gamma_{d}(t)=-\frac{dp(t)}{dt}/p(t)$. The decayed population $1-p(t)$ goes into the electronic ground state $\ket{g}\bra{g}$ where it evolves under the trap term.

The full system time dependent density matrix can be reconstructed as $\rho(t)=\rho_{g}(t)+\rho_{e}(t)$, where $\rho_{g}(t)$ is the decayed density matrix, and $\rho_{e}(t)=\ket{\psi(t)}\bra{\psi(t)}$ is the pure excited state density matrix. However, keeping track of $\rho_{g}(t)$ is cumbersome, as it is usually a mixed state that evolves under the trap Hamiltonian.

Nevertheless, having the full $\rho(t)$ is rarely needed. We show how various quantities of interest such as the accumulated energy, infidelity, and entanglement entropy can be calculated from the evolved electronic excited state $\ket{\psi(t)}$ alone.

The accumulated energy at time t, $E_{\text{R}}(t)\equiv\langle H_{\text{trap}}\rangle(t)$, can be calculated by adding the energy of the electronic excited state and the electronic ground state $E_{e}(t)+E_{g}(t)$. $E_{e}(t)$ is equal to $\bra{\psi(t)}H_{\text{trap}}\ket{\psi(t)}$, and $E_{g}(t)$ is accumulated over time from the recycling superoperator:

For highly subradiant configurations, storing a photon for a long time is of interest Asenjo-Garcia et al. (2017b). In particular, we wish to store the photon in a state of type $\ket{DS,0}$, where the electronic part is in a highly subradiant mode, $\ket{DS}$, and all atoms are in the vibrational ground state, $\ket{0}$. In the early dynamics of such a state ($t\ll 1/\gamma_{d}$), the decayed population, $\text{Tr}(\rho_{g}(t))$, is small and the system is mostly in the excited state.

However, due to the coupling between the electronic and vibrational states, the pure electronic eigenmode, $\ket{DS,0}$, generally turns into a superposition that includes other excited states that have different electronic or vibrational excitations. This leads to decoherence and infidelity in the stored excitation. The infidelity due to the mixing of the internally excited states is a measure of the closeness between the excited state at time $t$ and at time $0$ and is given by Nielsen and Chuang (2010)

where $\rho^{\text{proj}}(t)=\rho_{e}(t)/\text{Tr}(\rho_{e}(t))$. The reason we project the density matrix is to isolate the infidelity due to the mixing of the excited states from the infidelity due to the decay of the excited state. In the context of conditional measurements, if there is no photon detected before time $t$, the state of the system is $\rho^{\text{proj}}(t)$ and only the infidelity due to the mixing of the excited states is relevant.

In general, the electronic and vibrational parts of the state $\rho_{e}(t)$ become entangled before the state decays. The entanglement entropy is calculated from the reduced density matrix of the internal states of the atoms $\rho_{A}(t)=\text{Tr}_{V}\left[\rho^{\text{proj}}(t)\right]$, where $\text{Tr}_{V}$ is the trace over the vibrational states of the atoms. The entanglement entropy is given by Nielsen and Chuang (2010)

One goal of this paper is to identify the conditions under which the entanglement entropy and the infidelity are minimized.

We begin with the simplest example of highly subradiant states: two atoms in a 1D waveguide. In the single excitation subspace, the symmetrized combination can lead to subradiant or superradiant states when they are separated by a wavelength. The symmetrized electronic states are:

The effective Hamiltonian for this system is

In the single excitation symmetrized states, the effective Hamiltonian acting on $\ket{\pm}$ becomes Olmos and Lesanovsky (2025):

where $\phi=k_{0}d$ is the phase difference due to the separation of the traps.

For fixed point-like atoms located at the center of the traps, the operator $k_{0}\textbf{r}_{12}$ is zero, and the decay rate is

which leads to the constructive or destructive interference depending on $\phi$. For fixed point-like atoms, perfect dark states can be achieved by positioning the atoms in the $\ket{-}$ state separated by a multiple of $\lambda_{0}$ (or by an odd multiple of half-wavelength for the $\ket{+}$ state). We investigate the $\ket{-}$ state with integer $\lambda_{0}$ spacing as the dark state, and we use the effective Hamiltonian for the $\ket{-}$ state.

If the atoms are in a non-zero spread ground state of the trap, $\ket{-,0}$, the decay rate at integer wavelength separation becomes

where the term $\exp(-\eta^{2})$ accounts for the non-zero spatial spread of the wavefunction. This term arises from the convolution integral of the position probability distribution with the imaginary part of the Green’s function, as in Eq. (11). It prevents the perfect cancellation in Eq. (20), resulting in a lower bound for the decay rate $\sim\eta^{2}\gamma_{0}$, for small $\eta$ Rubies-Bigorda et al. (2024a).

For identically trapped atoms in the dark state configuration, the effective Hamiltonian, $H^{-}_{\text{eff}}$, in Eq. (18) can be decomposed into two components: (1) a center-of-mass Hamiltonian, $H_{\text{cm}}=\hbar\omega_{t}a_{cm}^{\dagger}a_{cm}$, representing a pure harmonic oscillator with no perturbation, and (2) a relative motion Hamiltonian, $H_{\text{rel}}$, which is a harmonic oscillator perturbed by dipole-dipole interactions:

where the non-indexed lowering operator, $a$, defined as $a=(a_{1}-a_{2})/\sqrt{2}$, acts on the relative displacement states. Here, $k_{0}\textbf{r}_{12}$ can be expressed as $\sqrt{2}\eta(a+a^{\dagger})$.

The relative Hamiltonian in Eq. (21) can be further expanded, assuming small $\eta$ and neglecting terms of order $\eta^{3}$ and higher, as

The terms in this expansion represent the trap of strength $\omega_{t}$ (first term), a linear potential of strength $\gamma_{0}\eta$ (second term), which pushes the wave packet away from the trap center, and non-uniform decay of strength $\gamma_{0}\eta^{2}$ (third term), which deforms the wave packet because the wave function’s tails decay faster than its center. Throughout this work, we demonstrate that the relative strength of these terms determines the behavior of decay, recoil, and decoherence in the dark states. It is useful to consider three trap regimes: a strong trap regime where the trap dominates the induced forces, $\omega_{t}\gg\gamma_{0}\eta$, a weak trap regime where the induced forces dominate, $\omega_{t}\ll\gamma_{0}\eta$, and an intermediate trap regime where the trap is only slightly stronger than the induced forces, $\omega_{t}\gtrapprox\gamma_{0}\eta$.

The effects of the induced forces and the non-uniform decay become manifest in weak traps. To isolate the effects of the different terms in (Eq. 22), we simulate three cases: 1) The full Hamiltonian as in Eq. (21) with the induced forces present. 2) The Hamiltonian without the induced real (coherent) potential term: $H_{\text{rel}}\rightarrow\hbar\omega_{t}a^{\dagger}a-\frac{\text{i}\hbar\gamma_{0}}{2}\{1-\cos\left(k_{0}\textbf{r}_{12}\right)\}$. 3) The Hamiltonian with the traps center shifted to cancel the linear potential term: $H_{\text{rel}}\rightarrow H_{\text{rel}}-{\hbar\gamma_{0}\eta}(a+a^{\dagger})/{\sqrt{2}}$. In other words, we cancel the linear part of the induced potential by adding a linear potential to the Hamiltonian.

The results of the simulations are in Figs. 2, 3, and 4. In these simulations, the mass of the atoms is $m=2.21\times 10^{-25}$ kg, and the individual atom decay rate is $\gamma_{0}=10^{8}$ s^-1. The internal transition has $\lambda_{0}=1000$ nm, and the trap frequency is $\omega_{t}=0.01\gamma_{0}=10^{6}$ s^-1 resulting in $\eta\approx 0.1$. The system is initialized in the dark state configuration $\ket{-}\ket{00}$ with $\lambda_{0}$ separation.

In Fig. 2, we plot the time dependent excited state population for the three cases. Initially, all cases exhibit an exponential decay with a rate of approximately $\gamma_{0}\eta^{2}$, but then deviate due to the change of the vibrational wavefunction. For cases where forces are absent or traps are shifted (cases 2 and 3), the decay rate slightly decreases. However, in the presence of forces (case 1), the decay rate greatly increases for $t>50/\gamma_{0}$.

The time dependence of the decay rate can be explained by the evolution of the spatial wavefunction. Figure 3 plots the mean position and width of the wavefunction over time. In the absence of forces (case 2) or with shifted traps (case 3), the wavepacket width decreases slightly due to faster decay of the wings, while the center of the wave packet remains nearly unchanged. This leads to a decrease in the standard deviation of position, $\sigma_{x_{12}}$, and thereby a decrease in the decay rate of Fig. 2. With forces present (case 1), the relative wavefunction shifts rightward, which is due to the repulsive interaction between the atoms’ individual wavepackets. The shift of the position of the wavefunction changes the relative phase $\phi$ from the initial value $2\pi$, resulting in a faster decay as per Eq. (19).

The motion of the wavepacket in Fig. 3 not only affects the decay rate but also results in the atoms gaining kinetic and potential energy before they decay. This contributes to the excited state energy term, $E_{e}$, discussed in Sec. II.5. This term builds up as the atoms are accelerated away from the trap center. This excited state energy is then transferred to the ground state energy term, $E_{g}$, as the atoms decay. When the atoms decay, they acquire additional recoil energy due to the emission of the photon. These two effects add up to the total energy of the atoms, $E_{\text{R}}$. As a consequence of these two contributions, the energy gained by the atoms exhibits qualitatively different behavior based on the trap strength.

In the strong trap regime, $\omega_{t}\gg\gamma_{0}\eta$ (or if the induced forces are absent or negligible), the induced forces are small, resulting in the excited state gaining minimal energy, $E_{e}(t)\approx 0$ for all $t$. In this case, the decay is exponential and the final energy is mostly from the recoil due to the emission of the photon. The total recoil energy of the atoms is inversely proportional to the decay rate, as noted in Ref. Suresh and Robicheaux (2021) leading to

where $E_{\text{r}}=\frac{\hbar^{2}k_{0}^{2}}{2m}$ is the one-atom recoil energy coming from momentum conservation during a photon emission from the state $\ket{e}$. The recoil energy in Eq. (23) after photon emission is derived from the recycling term, assuming a stationary initial excited state. For the dark state of two atoms in a 1D waveguide, the decay rate is approximately $\gamma_{0}\eta^{2}$, giving a recoil energy of

On average, a vibrational quantum is generated in this subradiant case in the strong trap regime. In Fig. 4, the motional energy gained by the atoms is shown for the three cases. When the forces are absent (case 2) or when the traps are shifted (case 3), the final energy is close to $\hbar\omega_{t}$, with difference ($E_{R}(\infty)-\hbar\omega_{t}$) coming from the narrowing of the wavepacket due to the non-uniform decay. However, with forces present (case 1), the final energy significantly exceeds $\hbar\omega_{t}$, mostly coming from the wavepacket displacement.

The energy gained by the atoms may be deposited either in the center-of-mass motion, $H_{\text{cm}}$, of the two atoms or in their relative motion, $H_{\text{rel}}$. For symmetric traps, a small fraction of the recoil energy, $E_{\text{r}}/2=\eta^{2}\hbar\omega_{t}/2$, is deposited in the center-of-mass mode.

The energy gained by the atoms and their decay behavior are manifest in the spectrum of the emitted photon Loudon (2000). The spectrum can be reconstructed following the ‘input-output’ formalism as in Refs. Gardiner and Collett (1985); Caneva et al. (2015); Sheremet et al. (2023). We calculated the spectrum of the emitted photon for the two atoms in a 1D waveguide in cases 1 and 2. In both cases the spectrum is red shifted on average by the final energy deposited in the atoms, $E_{\text{R}}(\infty)/\hbar$, and its width corresponds to the average decay rate of the subradiant state.

For a system of two atoms, the dark state $\ket{-}$ is an eigenstate of the effective Hamiltonian. Consequently, this state remains decoupled from the bright state $\ket{+}$ during the decay process, irrespective of atomic motion. This leads to the entanglement entropy $S(t)$, as defined in Eq. (15), remaining identically zero. However, at least two states are needed for a qubit. Therefore, we study a system of three atoms which can support two dark states.

In the single-excitation subspace, three atoms in a 1D waveguide, separated by $\lambda_{0}$, form two dark states and one bright state Asenjo-Garcia et al. (2017b); Paulisch et al. (2016); Rubies-Bigorda et al. (2024b). The dark states can be defined as:

and the bright state is:

If the atoms are fixed at integer wavelength separation, the dark states remain non-decaying while the bright state decays at a rate of $3\gamma_{0}$ Asenjo-Garcia et al. (2017b). When the atoms have a non-zero spread in a vibrational ground state, the dark states decay at rate of $\simeq\eta^{2}\gamma_{0}$, and the bright state decays at a rate of $\simeq(3-2\eta^{2})\gamma_{0}$ in the Lamb-Dicke regime. Using the two dark states as qubit states allows high fidelity for a duration less than the lifetime ($<1/(\eta^{2}\gamma_{0})$). The lifetime is long because these subradiant states are almost exact eigenstates of the system and decay at the same, slow rate Paulisch et al. (2016).

While the previous discussion focused on decay dynamics, practical traps introduce coupling between the two dark states due to atomic motion. For example, starting from $\ket{DS_{1},0}$, we observe population transfer to the electronic states $\ket{DS_{2}}$ and $\ket{BS}$, as shown in Fig. 5. In this simulation, we chose favorable parameters to preserve fidelity: $\eta=0.01$, $\omega_{t}=0.1\gamma_{0}$, and $d=\lambda_{0}$. The choice of $\eta=0.01$ corresponds to an atom 10 times heavier than a Cesium atom, while maintaining the same $\omega_{0}$ and $\gamma_{0}$ as before.

Unlike the complete decay dynamics in Fig. 2 for two atoms, Fig. 5 only shows the early dynamics where the decayed population, $1-p(t)$, is under $1.4\%$. For reference, if the atoms are confined to the ground state of the trap, the lifetime $1/\gamma_{d}$ is approximately $1/(\eta^{2}\gamma_{0})=10000/\gamma_{0}$, while we only simulate for $t_{f}=100/\gamma_{0}$. In this time window, the state $\ket{DS_{2}}$ (summed over all vibrational states) populates up to approximately $1.95\%$ and oscillates with a period close to the trap period $2\pi/\omega_{t}\approx 62.8/\gamma_{0}$. The bright state $\ket{BS}$ also populates to a small fraction $\approx 1\times 10^{-5}$. The simple sinusoidal behavior results because the dipole-dipole interaction couples $\ket{DS_{1},0}$ with $\ket{DS_{2}}$ in the $n=1$ subspace, and these states have an energy splitting $\simeq\hbar\omega_{t}$, see App. A.

This mixing between the states is only possible owing to the coupling between the internal and the vibrational states. In the process of the population transfer from the initial state, there is entanglement entropy ($S(t)$ as defined in Eq. (15)) generated. The entropy is not shown in the figure but is approximately proportional to the population in $\ket{DS_{2}}$. For example there is a peak in the entropy at $t\approx 30/\gamma_{0}$ which corresponds to the time when the population in $\ket{DS_{2}}$ is at its maximum. This peak is at $S(t)\approx 0.1$ which is small but could affect applications. Like the entropy, the populating $\ket{DS_{2}}$ causes infidelity, $I(t)$, that is proportional to the population in $\ket{DS_{2}}$ and also peaks at $t\approx 30/\gamma_{0}$ at a value approximately $0.01$.

The oscillation amplitude for the entropy or the infidelity curves can serve as a figure of merit to compare the effects of the trap strength and the induced forces across different system parameters, $\eta$ and $\omega_{t}$. In the intermediate or strong trap regime, we found that the amplitude of the oscillations depends on the ratio of the force term to the trap term, similar to the condition for exponential decay for two atoms. Specifically, the infidelity amplitude scales as $\sim(\gamma_{0}\eta/\omega_{t})^{2}$. This scaling, derived in Appendix A, arises because the force term, $\gamma_{0}\eta$, functions as a coupling in an effective two-level system, while the trap term, $\omega_{t}$, serves as a detuning. In the weak trap regime, such as in the case for Fig. 2, the mixing between the internal states is largest, and the entropy, $S(t)$ rapidly reaches a maximum value of $S_{\text{max}}\approx\text{ln}(2)$, which describes a completely mixed effective two-level system. To achieve an order $10^{-4}$ infidelity (or three nines of fidelity in the language of quantum information) for a case with 3 Cs atoms with the transition $6^{2}S_{1/2}\leftrightarrow 6^{2}P_{3/2}$, a relatively strong trap of frequency $\omega_{t}\approx 2\pi\times 3\,$ MHz $\approx\gamma_{0}/2$ is needed.

In this section, we present results for highly subradiant states of atoms interacting in free space, without a waveguide. For the free space case, a highly subradiant configuration of two atoms only occurs at nearly zero separation. For example, for two atoms polarized along the $z$ axis with a separation along the $x$ axis, the atoms need to be placed in the $\ket{-}$ state at $d\approx 0.1\lambda_{0}$ to achieve a $10$-fold decrease from the individual decay rate, $\gamma_{d}=\gamma_{0}/10$ Olmos and Lesanovsky (2025). This is impractical due to the rapid increase of induced dipole forces at such close distances. For instance, assuming spread $\eta=0.02$, a trap frequency $\omega_{t}\gtrapprox\gamma_{0}$ is necessary to prevent the atoms from colliding due to the attractive dipole-dipole forces. An alternative approach to suppressing the decay rate, without going to very small separations, is to arrange multiple atoms in a regular array with separation $d<0.5\lambda_{0}$ Asenjo-Garcia et al. (2017a), which is the main focus of this section.

For a large number, $N$, of point-like atoms arranged in a circular array in free space, the decay rate of the most subradiant eigenstate can be made very small by increasing $N$, provided that the separation between atoms is less than $\lambda_{0}/2$ Asenjo-Garcia et al. (2017a). In the infinite limit $N\rightarrow\infty$, the decay rates of many eigenstates approach zero. These eigenstates are electronic (spin) waves with different wavenumbers $k$. When $k$ exceeds the free-space wavenumber $k_{0}$, the radiation must be evanescent perpendicular to the array, resulting in guided modes around the ring that are completely dark. For finite even $N$, the most subradiant eigenmode is the one with the largest $k$ which is equal to $\pi/d$. This state which we call $\ket{DS_{min}}$ is an equal superposition of all single excitation states with alternating signs. The decay rate of this state, $\gamma_{min}$, decreases exponentially as $\exp(-N/N_{0})\gamma_{0}$, where $N_{0}$ is a constant that depends on the interatom distance, $d$. The resulting highly subradiant states of the ring are promising for the storage of photons and quantum information.

If the atoms are spread in the ground state of a trapping potential, as shown in Fig. 1, the decay rate of the most subradiant state has two primary contributions: (1) the finite size of the array, which scales like $\sim\exp(-N/N_{0})\gamma_{0}$, and (2) the finite spread of the atoms’ wavefunction, which contributes a factor of $\sim\eta^{2}\gamma_{0}$. At sufficiently large $N$, $\gamma_{min}$ reaches an asymptotic value dominated by the atomic spread, $\gamma_{spread}$, as illustrated in Fig. 6 and derived in Appendix B. However, as $d$ increases, more atoms are required to reach this spread-dominated regime. In the case in Fig. 6, the atoms are taken to have position spread in the plane of the atoms (x-y plane), but no spread in the perpendicular direction. Also, their polarization is perpendicular to the plane of the atoms ($\hat{q}=\hat{z}$). In this case, $\gamma_{spread}\approx 0.8\eta^{2}\gamma_{0}$. If the polarization was circular around the z-axis ($\hat{q}=(\hat{x}+\text{i}\hat{y})/\sqrt{2}$) instead, the limiting decay rate, $\gamma_{spread}$, becomes $0.6\eta^{2}\gamma_{0}$.

These minimum decay rates were calculated by diagonalizing the non-Hermitian effective Hamiltonian (restricted to $n_{max}=0$ so all atoms are in the ground vibrational state). An alternative way is to utilize the ring symmetry Antman Ron et al. (2024) which results in the spin waves always being the eigenmodes of the system independent of $\eta$ and $d$. In this way, $\gamma_{min}$ is equal to $2\,\text{Im}\{\bra{DS_{min}}H_{\text{eff}}\ket{DS_{min}}\}/\hbar$. A third way of calculation is to use statistical methods. In the calculation of $\gamma_{min}$, we found that having vibrational ground states is approximately equivalent to point-like atoms that have a Gaussian position distribution equal to that of the vibrational ground state. Computationally, we sample a position configuration from the Gaussian distribution, which gives a particular $H_{\text{eff}}$, and calculate $2\,\text{Im}\{\bra{DS_{min}}H_{\text{eff}}\ket{DS_{min}}/\hbar$. We repeat this process many times and average the results.

Similar to the case of three atoms, forces between atoms in the ring induce entanglement between their internal and vibrational degrees of freedom. Figure 7 shows the mixing dynamics starting from $\ket{DS_{min},0}$. The parameters used are $\omega_{t}=0.1\gamma_{0}$, $\eta=0.01$, $d=0.3\lambda_{0}$, $\hat{q}=\hat{z}$, and $N=30$. In the time shown, approximately $0.7\%$ of the population has decayed. A small fraction, about $0.04\%$, transfers from the initial internal state to other internal states in the first vibrational mode, $n=1$, (summed over the all the internal states). This process generates entanglement entropy $S(t)$ and produces infidelity $I(t)$. Compared to the 3 atoms case in Fig. 5 which has the same $\eta$ and $\omega_{t}$, the infidelity here is much smaller, of order $10^{-4}$, compared to $10^{-2}$ for the waveguide atoms. This can be explained by the magnitude of the induced forces between the atoms. For a waveguide at $\lambda_{0}$ separation, the force is $\hbar\gamma_{0}k_{0}$. For free space at $d=0.3\lambda_{0}$ and when the polarization is perpendicular to the plane of the atoms, the force is much smaller, approximately $0.03\hbar\gamma_{0}k_{0}$. This is a particularly useful arrangement for storing quantum information, as the decay rate is small and the induced forces are weak. The force stays small for $d/\lambda_{0}$ in the range from approximately $0.25$ to $0.5$ but increases rapidly for smaller $d$, as later indicated in Fig. 8. For several other polarization directions, the induced force was found to be larger.

Unlike the periodic mixing seen for three atoms in Fig. 5, the mixing here is not periodic at early times. The lack of periodicity can be explained by the change of energy splitting between the interacting states. For atoms in a waveguide, the splitting is $\omega_{t}$, which causes oscillations with frequency close to $\omega_{t}$. While being absent in a waveguide at $d=\lambda_{0}$, another contribution to the splitting comes from the Lamb shift due the collective coupling in free space at $d=0.3\lambda_{0}$. The Lamb shift in this case is approx $-0.25\gamma_{0}$ for the initial state and $-0.32\gamma_{0}$ for a typical strongly interacting state in a higher vibrational mode. This uneven Lamb shift causes the energy splitting between the interacting states to diminish from $\omega_{t}=0.1\gamma_{0}$ (the trap splitting) to only $\approx 0.03\gamma_{0}$. This results in a slower oscillation and a revival time longer than the time window shown in Fig.  7.

Due to the increasing size of the Hilbert space with $N$, this simulation uses a restricted vibrational Hilbert space allowing at most one atom to vibrate in one direction at a time. Due to the large number of states in these calculations, we did not check the validity of this vibrational restriction. However, this restriction was checked for validity in the case of 3 atoms in Fig. 5, which was relatively easier to perform due to the fewer number of states. We found that the errors in the populations to be less than $0.1\%$ between this restriction and a convergent calculation that used $N_{vib}=5$. The errors were even smaller for cases with larger $\omega_{t}$ or smaller $\eta$. This suggests that such restriction results in a good approximation in the intermediate trapping regime where vibrational excitations in the $n=1$ state are less than $\sim 10\%$. Further details of the equations and approximations used are given in the App. B.