Light–matter interaction in subwavelength scale is of broad interest in quantum information processing and quantum metrology [1,2]. Such interactions in subwavelength atomic arrays [3–28] can induce long-range nonlinear dipole–dipole interactions by mediating photons in radiation modes [3,5,29,30] or via long-range van der Waals interactions between Rydberg states [31]. Further controls of atomic arrays by magnetic fields bring remarkable physical phenomena, including super- and subradiant states [32], photon storage and retrieval [33–35], subradiance-protected quantum state transport [36], and many others [37–45]. These quantum information processing methods utilize subradiant quantum states with inhibited spontaneous decay and exhibit useful applications in quantum storage [46]. Previous research shows that single-excitation subradiant states in one-dimensional (1D) arrays of N atoms typically have decay rates scaling as N−3 with a lifetime of 7 orders of magnitude greater than that of a single atom for N∼200 [6,11,15]. Nevertheless, it generally desires a quantum state lasting long enough for the purpose of quantum storage [46]. A challenging problem is that subradiant states can diffuse into the superradiant subspace through the interaction or boundary effect, which limits their applications in quantum information processing.

Here we show that a judicious coherent control of finite-size atomic arrays can prevent single-excitation subradiant states from entering the superradiant subspace and thus realize subradiance with a lifetime of orders of magnitude longer than that of a single atom with the help of an external magnetic field. Specifically, we consider a well-established model with atomic arrays under the magnetic field, which lifts the degeneracy of excited states in each atom through the Zeeman shift (labeled as |+⟩ and |−⟩) [37,38,41]. Such systems with an infinite size can support subradiant collective excited states outside of the dissipative zone in the momentum space [6,11,15]. However, for the same system with a finite size, such excitation inevitably spreads to the boundaries and the wavepacket gets scattered into the dissipative zone. To overcome such obstacles, we instead apply a spatially linear magnetic field, which brings linear potentials with opposite slopes on |+⟩ and |−⟩ states [see Fig. 1(b)]. Similar 1D lattices under linear potentials support Bloch oscillations [47–58]. Fundamentally different from the previous work, Bloch oscillations in atomic arrays under a linear magnetic field support two-band Bloch oscillations with inter-band interactions such that the wavepacket of excited states oscillates on two bands along opposite directions alternatively. This unique picture can force the wavepacket spatially localized within a few atoms while preventing the excitation from entering into the dissipative zone. Hence it supports subradient oscillations surviving extremely long time with finite atoms. We also show the way of switching between superradiant and subradiant oscillations by tuning the external magnetic field, which provides the opportunity for reading out information from subradiant quantum states [59–62].

The developments of state-of-the-art technologies in various experimental platforms such as trapped neutral atoms at the subwavelength scale [19–28,63–68] pave the way for exploring exotic phenomena in atomic arrays. We study such arrays of atoms that are equally spaced at a distance a along the y-axis [yn=na for the nth atom in Fig. 1(a)]. Each atom has a ground state |g⟩ and degenerate excited states. We apply an off-plane magnetic field B(y) in the z-direction to lift the degeneracy of two excited states |±⟩. Under the single-excitation limit, the dynamics of the wave packet of the excited states is described by the non-Hermitian Hamiltonian after integrating out photonic modes [37,38,41], $$\begin{gathered} Heff=ℏ∑n=1N∑α=±(ωA−iγ02)|αn⟩⟨αn|+∑n=1NμBn(|+n⟩⟨+n|−|−n⟩⟨−n|) \\ +3πℏγ0cωA∑n≠m∑α,β=±Gαβ(yn−ym)|αn⟩⟨βm|, \end{gathered}$$(1)where ωA is the atomic transition frequency, γ0 is the atomic decay rate in the free space, μBn/ℏ is the Zeeman shift for the nth atom with the magnetic moment μ, and Gαβ(yn−ym) is the free-space dyadic Green’s function describing the electric field at yn emitted by the atom located at ym. With the relation |±⟩=∓(|x⟩±i|y⟩)/2, one obtains G++=G−−=−eik0r·(k02r2−ik0r+1)/(8πk02r3) and G+−=G−+=eik0r(k02r2+3ik0r−3)/(8πk02r3), where r=|yn−ym| and k0=ωA/c=2π/λ with λ being the wavelength and c being the vacuum speed of light [38,41]. Equation (1) shows a tight-binding lattice model where excitations on sites at two arms (|±⟩) have on-site potentials V±,n≡±μBn and are connected by complex and long-range photon-mediated hoppings. In particular, we consider the magnetic field Bn=nB0,(2)where B0 is a constant. Equation (2) gives V±,n=±nμB0, which leads to the effective constant force F±=∓μB0/a in opposite directions on the two excited states. Two effective electric fields in opposite directions on two separated arms in the lattice bring Bloch oscillations on two arms exhibiting symmetric patterns. However, once the two arms are connected [by terms including G±∓ in Eq. (1)], the dynamics of the Bloch-oscillation wavepackets on the two arms is influenced by each other, which can be analyzed in the field-free Bloch band picture.

We first study the band structure of the infinite atomic arrays with constant magnetic field Bc [37,38]. The single-excitation Bloch modes are given by |Ψky⟩=∑n(C+,ky|+n⟩+C−,ky|−n⟩)eikyna=(C+,kyC−,ky)T,(3)where ky is the Bloch wave vector. The Hamiltonian in the momentum space on two arms has the form Hkℏ=(ωA−iγ02+χky++χky+−χky−+ωA−iγ02+χky−−)+(iμB0ℏaddky00−iμB0ℏaddky),(4)where χkyαβ=3πγ0cωA∑n≠0eikynaGαβ(na) and diagonalization of the first term plus the shifts caused by the constant magnetic field gives the two-band structure. Our results show that the oscillatory mode for an extended Gaussian distribution of the single-excitation wavepacket approximately evolves in a single band. The second term in Eq. (4) would transform the distribution on one band into two bands and gives the acceleration force in one band and the tunneling strength to the other band, which surpasses the former one if the two arms have comparable probabilities. We thus bring in a phenomenological force averaging on two arms, ⟨F⟩≡μB0a∑n⟨Ψ|(|+n⟩⟨+n|−|−n⟩⟨−n|)|Ψ⟩=μB0a∑n(P+−P−), to estimate the acceleration force.

We choose a=0.1λ, which brings the first Brillouin zone ky∈[−5k0,5k0]. We plot band structures in ky-space for different Bc in Fig. 1(c). For each Bc, there are two bands, and it exhibits large collective decay, i.e., the dissipative zone in the free-space light cone (|ky|<k0). Outside of the dissipative zone, subradiant states are supported. For Bc=0, there are two degenerate points near ky=±2.3k0, while two bands are separated in the entire ky-space for other Bc. We denote the lower (upper) band as band I (II) and define P±I(II) as the probability intensity of an eigenstate on band I (II) projected onto the arm of |±⟩ states. P+II and P+I are shown in Figs. 1(d) and 1(e), respectively. One notes that P+I(II)+P−I(II)=1. When there is no magnetic field, P±I=P±II=0.5 in the entire ky-space, indicating that probabilities for exciting |±⟩ states are equivalent due to the degeneracy of |±⟩. However, for m≠0, it shows asymmetric distributions of P±I(II). In particular, for Bc>0 on band I (II), we can see that P+I<0.5 (P+II>0.5), meaning that less (more) density of the state is located on the arm of |+⟩ states. Distributions of P±I(II) are highly symmetric because the influence of |±⟩ is symmetric under the positive/negative magnetic field. Moreover, we find that, for Bc>0 at ky=±2.3k0, it shows P+I=0 (P+II=1) on band I (II). We note that band structures in Fig. 1 are calculated by simply summing the Green’s function in real space over all lattice sites, which artificially induces the inaccurate eigenvalues at ky near ±k0 in the free-space light cone due to slow convergence. However, this inaccuracy in band structures does not affect the later analysis throughout this paper. Also, if desired, it can be overcome by summations of Green’s function in momentum space with an appropriate regularizing method [37,38].

We next introduce the simulation method for studying the dynamics of the single excitation wavepacket in the atomic arrays. We use the Schrödinger equation d|Ψ(t)⟩/dt=−iHeff|Ψ(t)⟩/ℏ and the excitation wavepacket of the atomic array |Ψ(t)⟩=∑n[C+,n(t)|+n⟩+C−,n(t)|−n⟩]e−iωAt, where |C±,n|2 gives the excitation probability of the |±⟩ state in the nth atom. Two hundred one atoms (n=−100,…,0,…,100) are considered in the arrays. B0=0.2ℏγ0/μ is used in simulations. We assume that the wavepacket of the system is initially prepared at a superposition state on two arms following C±,n(0)=ψ±exp[ikcan−(n−nc)2/200], where nc is the spatial center of the initial excitation, kc is the initial momentum, and ψ± gives the initial ratio of excitation amplitudes on two arms of the |±⟩ states. Such an initial state in the weak excitation limit can be pumped with phase-controlled schemes [32,33,59,69] or by applying a spatial modulation of the atomic detuning [18].

For finite atomic arrays, the key ingredient for maintaining a long-standing subradiant state is to prevent the wavepacket from entering into the dissipative zone either by oscillations or scattering due to the boundary effects. One might notice that the well-established Bloch oscillation may lead the wavepacket oscillation within a finite spatial region and propagating unidirectionally in the ky-space [50,53,70], so by frequently altering the direction of the constant force, it is possible to accomplish the task for avoiding entering the dissipative zone. Nevertheless, here we show that one can realize such a task in atomic arrays under a spatially linear but temporally constant magnetic field [Eq. (2)], which leads to effective constant forces in opposite directions on the two excited states. Once the two arms are connected via photon-mediated dipole–dipole interactions, the effective electric force that acts on the single-excitation wavepacket depends on the wavepacket distributions between the two arms. Our results show that an extended Gaussian distribution of the single-excitation wavepacket approximately evolves in a single band driven by the phenomenological force ⟨F⟩, which can be estimated by the eigenstate distributions of the moving wavepacket in a band structure under a constant magnetic field felt by the wavepacket center. The ability of altering the effective force direction during Bloch oscillations provides the mechanism to build subradiant states that avoid entering the dissipative zone inside the free-space light cone and are spatially localized far away from the boundary of the finite array.

To see the long-standing subradiant state, we excite the wavepacket centered at n=0, i.e., nc=0, with ψ+=ψ−=0.168 and kc=1.5k0, which satisfies ∑n[|C+,n(0)|2+|C−,n(0)|2]=1. The evolutions in the simulation are plotted for the time up to 10TB, where TB=2πℏ/μB0=10π/γ0. Figures 2(a1) and 2(b1) depict the dynamics of |C+,n(t)|2 and |C−,n(t)|2, respectively, and one sees a sustained Bloch oscillation pattern with the shape of the wavepacket deforming gradually. The excitation of the |+⟩ (|−⟩) states in the atomic arrays is mainly located in the finite region with y<0 (y>0) so there is no boundary effects due to finite atoms.

We then perform fast Fourier transform (FFT) on the simulation results C±,n(t) and calculate the projection of the excited wavepacket in ky-space onto the two bands [excitation probabilities PI(ky,t) and PII(ky,t), respectively]. Note that the Bloch band picture is valid only with constant magnetic field. Here, the linear magnetic field and resulting constant force produce perturbation, so we can calculate the projection by taking into account the uniform potentials at the value of V±,n=±n¯μB0, where n¯≡⟨Ψ|∑n(|+n⟩n⟨+n|+|−n⟩n⟨−n|)|Ψ⟩ is the mean position of the wavepacket. This analysis in the Bloch band picture with uniform on-site potentials reveals dynamical features different from a conventional Bloch oscillation. Figures 2(c1) and 2(d1) show the evolution of excitation probabilities of band I and band II and PI and PII at each ky. The choice of ψ± makes the initial excitation on band I, and throughout the evolution, PII is about 1 order of magnitude smaller than PI. Remarkably, one sees that PI does not evolve unidirectionally on ky with time, which gives the fundamental difference from the conventional Bloch oscillations. The evolution direction of PI reverses every time (referred to as the reverse time) before it enters the free-space light cone in the ky-axis [dash lines in Figs. 2(c1) and 2(d1)], which leads to the important consequence of avoiding a large collective decay of the excited wavepacket. We further show the total excitation probabilities of the |+⟩ state and the |−⟩ state (P+, P−) in the real space and the total excitation probability Pt=P++P− in Fig. 2(e1). We find that the times when P+=P− are exactly the same as the reverse times because the phenomenological force ⟨F⟩ becomes zero and subsequently changes its direction. In other words, when the center of the wavepacket arrives at n=0, P+≈P− and hence ⟨F⟩ changes its direction, resulting in the reverse of PI(ky) in the momentum space. Moreover, Pt exhibits no decay in this dissipative atomic system, which is actually lasting ∼108γ0−1 given by the simulation with a much longer evolution time. Therefore, as long as the excited wavepacket in the atomic arrays does not enter the dissipative zone inside the free-space light cone, the life time of the excitation gives the subradiant feature, which is much longer than the free-space decay time. Such extreme subradiant oscillations can be optimized by choosing initial excitations. In particular, we find that Pt can last more than ∼1013γ0−1 for an initial state centered at kc=3k0 on band I. Based on the eigenvalue analysis of the effective Hamiltonian for an array under a linear magnetic field, one can determine the decay rate of the most subradiant eigenstate, which is consistent with the property of the optimized subradiant Bloch oscillations (see Appendix A for details).

As one may notice, the initial choice of the excited wavepacket determines how the wavepacket evolves in ky-space and hence how the wavepacket decays. For example, if we choose ψ+=0.168, ψ−=−0.168, and kc=4k0 to excite band I or ψ+=0.168, ψ−=−0.168, and kc=2k0 to excite band II, we find the total probability of the wavepacket undergoes no decay even though the evolution patterns are different (see Fig. 5 in Appendix A). However, in the case of ψ+=ψ−=0.168 and kc=4k0, with band II being excited, the dramatically different dynamics is as shown in Figs. 2(a2)–2(e2). One notes that the excitation undergoes a rapid decay after the partial oscillation during the initial time duration ∼11γ0−1. The fundamental difference here is that the excitation probability PII moves into the free-space light cone along the ky-axis before ⟨F⟩ changes its direction. This result leads to a very fast decay within ∼8γ0−1 for the total probability Pt dropping to ∼0.02. In real space, the excitation on the |−⟩ arm decreases fast, while the excitation on the |+⟩ arm increases initially with an oscillation towards the positive y-axis but decays rapidly afterwards.

The two-band Bloch oscillation in the atomic arrays under the linear magnetic field brings unique opportunity for manipulating the decay of the excited wavepacket. For instance, we can control the collective decay exhibiting quantization-like decay with stable plateaus, as shown in Fig. 3(c1). To achieve it, we choose ψ+=ψ−=0.168 and kc=k0, so the wavepacket on band I centered at k0 is initially excited. The simulation results are summarized in Figs. 3(a1)–3(c1). The wavepacket experiences decay at t∼0 because the initial wavepacket is partially inside the free-space light cone. However, the wavepacket initially moves toward the positive direction in the ky-axis [see PI(ky)+PII(ky) in Fig. 3(b1)], so Pt ends in a stable plateau at t∼2γ0−1 with a value ∼0.67. At t∼27γ0−1, the excited wavepacket is reached at ky=−k0, and a reversal happens due to the change of ⟨F⟩, which results in a rapid drop of Pt to ∼0.48. Periodic oscillations are exhibited, together with the oscillations of Pn=|C+,n|2+|C−,n|2, in the real space. The total excitation probability Pt therefore shows the overall decay tendency with the periodic plateaus [see Fig. 3(c1)]. Initial excitations inside the free-space light cone can exhibit similar oscillations but quicker decay. Particularly, for an initial wavepacket at kc=−0.5k0 on band II, Pt decays to ∼0.02 within t∼0.85γ0−1, or a collective superradiant decay rate of ∼4.7γ0. Moreover, excitations with nc≠0 can bring different Block oscillation phenomena because the central momenta at reverse times are changed during the evolution of the wavepacket (see Appendix A for details).

Taking advantage of the relation between reverse times and the atomic position in the linear magnetic field, one can control the momentum of the wavepacket at the reverse time. The magnetic field acts as an effective field on the single particle and breaks the spatial translation symmetry, which leads to the particle movement in the momentum space. Change of the magnetic field would alter the effective force and the momentum space evolution of the wavepacket. We give an example in Figs. 3(a2)–3(c2), where the spatial distribution of the linear magnetic field relative to the atomic position is tuned to control the evolution of the wavepacket that is initially excited at k0 on band I [the same as that in Figs. 3(a1)–3(c1)]. In the numerical simulation, we tune the zero point of the magnetic field from y=0 to y=30a over a time period from 0 to 20γ0−1. As a result, after Pt drops to the first plateau, the wavepacket oscillates around the new center y=30a [Fig. 3(a2)] and reverses at ∼±1.4k0 in the momentum space to maintain subradiance [Fig. 3(b2)]. The excitation can also be driven to the dissipative zone by moving the zero point of the magnetic field back to y=0 over a time period from 180γ0−1 to 200γ0−1 and decay radiatively, shown by descending plateaus in the total excitation probability.

The proposed atomic arrays can be realized in various experimental systems. Dependent on different experimental platforms, atomic arrays can be composed with atom ensembles in the gas phase, such as neutral particles, including Rydberg atoms trapped in optical lattices [19,21–23,63,64] and optical tweezers [20,24–28,65–68], or atom-like particles in the solid phase, such as color defects in diamond nanophotonic devices [71,72], and excitons in atomically thin semiconductors [73]. Considering the possible imperfection of atomic positions and the inhomogeneous broadening in the atomic array, numerical results show that positional disorder suppresses the positive influence of the magnetic field on subradiance (see Appendix B for details). This limitation is expected to be overcome by the development of experimental control techniques. On the other hand, the effect of decay suppression can be found in arrays with different atomic numbers and interatomic distances to possibly decrease experimental challenges. The performance of subradiance relies on the magnitude of the linear magnetic field. The parameters used above present better initial state selections and characteristic performance of Bloch oscillations, which may further be optimized for extreme subradiance. Moreover, it has been shown recently that dispersionless or trapped subradiant photon states can be created in one-dimensional emitter chains [74]. Different from Ref. [74], our study shows that, in the model with the external magnetic field it is possible to decrease the subradiant decay rate to smaller than that in the atomic arrays without the magnetic field and also points out the possibility of quantum-state switching from the external-field control.

In summary, we have explored single-excitation dynamics in 1D atomic arrays under a linear magnetic field. Resulting two-band Bloch oscillations are studied. We found that the evolution of the wavepacket in the momentum space depends on the phenomenological force ⟨F⟩, which is determined by local properties of the band structure and thus influenced by the position of the wavepacket in the real space. Compared to the conventional Bloch oscillations, the evolution direction of the momentum reverses as long as the wavepacket passes the position where ⟨F⟩ changes its sign. Taking advantage of this unique property, we show the capability of generating exotic subradiant oscillations that lasts orders of magnitude greater than the atomic spontaneous decay time. Our study takes advantage of the external control from the magnetic field and therefore points toward fundamental opportunities in switching quantum states between superradiant and subradiant states and realizing an extreme subradiant quantum state that could be useful in important applications for quantum storage.

Acknowledgment. L.Y. thanks the sponsorship from the Yangyang Development Fund.