Exploiting the Photonic Non-linearity of Free Space Subwavelength Arrays of Atoms

Ordered ensembles of atoms, such as atomic arrays, exhibit distinctive features from their disordered counterpart. In particular, while collective modes in disordered ensembles show a linear optical response, collective subradiant excitations of subwavelength arrays are endowed with an intrinsic non-linearity. Such non-linearity has both a coherent and a dissipative component: two excitations propagating in the array scatter off each other leading to formation of correlations and to emission into free space modes. We show how to take advantage of such non-linearity to coherently prepare a single excitation in a subradiant (dark) collective state of a one dimensional array as well as to perform an entangling operation on dark states of parallel arrays. We discuss the main source of errors represented by disorder introduced by atomic center-of-mass fluctuations, and we propose a practical way to mitigate its effects.


I. INTRODUCTION
Ordered atomic ensembles (arrays) have recently attracted significant attention as a new paradigm for controlling light-matter interaction which shows novel features not shared by their disordered counterpart [1]. When the interatomic separation is subwavelength with respect to the characteristic atomic dipole transition, the optical response of atomic arrays shows a strong collective behaviour characterised by bright (superradiant) and dark (subradiant) excitations [2,3]. Superradiant states allow for efficient coupling of internal atomic states to light, while subradiant states permit long coherent storage of atomic excitations thanks to their reduced linewidth. Additionally, the array's collective response can be used to realise perfect reflection of light off the array [4][5][6][7][8] and to prepare topological edge modes [9,10].
This combination of features contains some of the basic elements for applications in quantum information technologies, as one could store quantum information in subradiant states [11] and read it out using superradiant excitations [12]. However, for the creation and manipulation of quantum information one also requires a nonlinear optical response, i.e. the dependence of lightmatter interaction on the system's internal state. Specifically, a non linear response allows to define qubits and to perform universal gate operations on them. Optical nonlinearities are usually obtained by adding a new capability to the system. In the context of an atomic ensemble, for example, a non-linear response is obtained promoting atoms to a Rydberg state: strong dipolar interactions between two Rydberg excitations can be exploited to inhibit further absorption of photons from the ensemble (Rydberg Blockade) [13]. Rydberg excitations have also recently been considered for atomc arrays [14,15] In this article, we propose and analyze an alterna-tive way to produce and coherently manipulate quantum information with ordered atomic ensembles, which does not rely on Rydberg excitations or other technologies. We show how one can harness the intrinsic nonlinear response of subwavelength arrays of atoms in free space [16][17][18][19] to perform different tasks. Specifically, we discuss (i) a procedure to transfer an excitation from the ground to the single-excitation most-subradiant collective state of the system and (ii) a procedure to prepare an entangled state shared by two parallel arrays. Together these two tasks allow to perform a universal set of quantum gates on atomic arrays. We discuss the nature of this intrinsic non-linear response and show that it has both a coherent and a dissipative component. The coherent component is rooted in the large dipole shift between closely spaced atoms; It dominates at smaller lattice spacing but diminishes for longer arrays. The dissipative component arises from the difference in the decay rate of the subradiant mode and the (enhanced) decay rate of the double excited mode, and allows for coherent manipulation via the Zeno effect [20]. Notably, for the case when atoms are pinned to their position, the dissipative non-linearity has more prominent effects for longer arrays. Finally, we analyze the impact of atomic center-of-mass fluctuations on the proposed scheme. We find they represent the main source of imperfections for they are responsible for a suppression of the intrinsic dissipative non-linearity [4,21]. We show numerically that such detrimental effects are reduced when the time scale of atomic motion is shorter than the internal decay rate, and propose a practical way in which this regime can be obtained.
The article is organised as follows. We first consider the case where atoms are pinned to their lattice positions (Sec. II). Within this simplified assumption, we present the main ideas and discuss how to perform single arXiv:2107.00566v3 [quant-ph] 23 Apr 2023 and two-qubit gates between the ground and the singleexcitation subradiant state of atomic arrays. In Sec. III, we consider the effect of the atoms' motion. We propose a way to enter the regime of fast atomic motion, and analyze the gate fidelities in this case. In Sec. IV, we discuss a possible experimental implementation of our proposal in the context of neutral atoms in optical lattice and analyze limitations and additional assumptions behind our model. We draw our conclusions in Sec. V. Additional non-essential details are left to the appendices.

II. ARRAYS OF PINNED ATOMS
We consider two parallel atomic arrays labelled "A" and "B" and placed at a distance l from each other. Each array contains N atoms separated by a lattice spacing a (see Fig. 1.a). In the following, we consider the case in which the arrays have unit filling and the atoms are pinned to their lattice positions. The atoms' internal structure is described by a lambda scheme with one excited state |e coupled by a dipole-allowed transition to the ground state levels |g 1 and |g 2 ( Fig. 1.a). An external laser drives the atoms on the |e ↔ |g 2 transition with a detuning ∆ and at a rate Ω. For ε ≡ Ω/2∆ 1, the excited state |e is never populated and the atom behaves as an effective two-level system with excited state |g 2 and ground state |g 1 , characterised by an effective decay rate Γ 0 ≡ ε 2 γ 1 and dephasing rate κ 0 ≡ ε 2 γ 2 . We further assume |e ↔ |g 1 to be a much stronger transition than |e ↔ |g 2 , so that κ 0 Γ 0 . In the following, we thus neglect the dephsing arising when driving the |e → |g 2 transition and discuss later under which conditions this approximation can be justified (see Sec. IV). The additional complication introduced by the lambda scheme, as opposed to using simple dipole-coupled levels, is instrumental for mitigating the effects of the motion as explained in Sec. III [22].
The effective dynamics of the atoms treated as an open quantum system can be derived from the total Hamiltonian describing the dynamics of the atom-light interaction in the dipole approximation. Tracing out the electromagnetic field within the Born-Markov approximation, the dynamics of the system is described by the following non-hermitian Hamiltonian [3] (see also Appendix A) The non-hermitian HamiltonianĤ ν describes the dynamics of the atoms of array ν = A,B, and readŝ where ω g is the splitting between |g 2 and |g 1 ,σ + νj = (σ − νj ) † ≡ |g 2,j ν g 1,j | and j = 1, . . . N labels the atom at position R Aj = (0, 0, ja) T (R Bj = (0, l, ja) T ) within array A (B). Here, J ij (Γ ij /2) is the coherent (dissipative) interaction between two atoms at sites i and j within the same array, and for i = j it is given by [3] J ij − i Γ ij 2 = − µ 0 ω 2 e d νi · G 0 (R νi − R νj , ω e ) · d νj (3) where G 0 (R, ω e ) is the free space electromagnetic Green's tensor of a point dipole, ω e the frequency of the |e ↔ |g 1 transition, and d i is the electric dipole moment of the atom at site i. For i = j, we defined Γ ii = γ 1 and we implicitly included the vacuum Lamb-shift J ii into the definition of ω g . Hereafter, we consider all atoms to be polarized parallel to the array's direction (z-axis) unless otherwise specified. A different atomic polarization leads to results qualitatively similar to the ones presented here, as shown in Appendix G. The last term in Eq. (1) describes the interaction between the atoms in the two arrays and readŝ where g ij (γ ij /2) is the coherent (dissipative) part of the dipole coupling given by Note that as a consequence of the Raman transition used to define the effective two level atoms, the dipole-dipole interactions in Eq. (2) and Eq. (4) are proportional to the Green's tensor evaluated at ω e . Hence, despite the effective two level system having a characteristic wavelength λ g ≡ ω g /2πc, the array exhibit strong collective dipolar effects only for a < λ e /2 (subwavelength condition [6,11,12,23]). To better understand the origin of the non-linear response of atomic arrays and how one can use it for preparing the arrays in a specific quantum state it is useful to consider first the case of a single array.

A. Single Array: Driving Subradiant Excitations
The dynamics of an isolated array of atoms in free space is described by Eq. (2). Because we consider only one array, we drop the label ν = A,B and refer to the single array Hamiltonian asĤ 1array in this section. The eigenstates of Eq. (2) are simultaneous eigenstates ofN e ≡ N jσ + jσ − j , as both the hermitian and anti-hermitian part of Eq. (2) commute withN e . Within the single excitation manifold and when the atomic array is sufficiently long (N 1), the eigenmodes of Eq. (2) can be understood as spin waves of a definite quasi-momentum k, where the value k corresponds to the point in reciprocal space where the eigenmode wavefunction is peaked [11,12,23,24]. We define the eigenmode with associated quasi-momentum k as |k ≡Σ + k |0 , where |0 ≡ |g 1 ⊗N is the ground state of Eq. (2), and For a < λ e /2, there exist single excitation eigenstates of Eq. (2) with a quasimomentum q > ω e /c ≡ k e lying outside the light-cone of free space electromagnetic modes (hereafter q labels values of the quasi-momentum which lie outside the light cone unless otherwise specified). These collective states, which are intuitively understood as an excitation propagating along the array, have been shown to decay at a rate ∼ Γ 0 /N 3 [11,25], and are referred to as collective subradiant modes or dark modes. Dark modes also exist in higher n-excitation manifolds, provided n N , and are similarly characterized by a scaling of the decay rate ∼ Γ 0 /N 3 . Notably, the n-excitation states |nq ≡ (Σ + q ) n |0 are not eigenstates of Eq. (2). This results in a non-linear structure of the dark-mode spectrum which has both a coherent and dissipative component. The coherent component is represented by a non-linear spacing of the energy levels with the number of excitation and it is quantified by the difference ∆ n ≡ ω n − nω 1 , where ω n (ω 1 ) is the energy of the most-subradiant nexcitation (single excitation) state. The dissipative component is represented by the enhanced decay rate of |nq which scales as ∼ Γ 0 /N to first order in Eq. (2) [11]. We will show that ∆ n approaches zero for increasing N , while the dissipative non-linearity grows with the array's size.
To excite a collective dark mode of a subwavelength array, it is necessary for the driving to match energy and modulation of the target state [26]. This condition can be achieved using a detuned Raman transition via an additional excited state |f . As illustrated in Fig. 1.b, we consider two driving lasers with wave vectors k a , k b forming respectively angles α, β with the array's direction (z-axis). For a sufficiently large detuning ∆ f such that |f is never populated, the effect of the Raman lasers on the array can be modelled by an effective driving Hamiltonian Here, z j ≡ aj, ω d ≡ ω a − ω b , and Ω 0 is the single atom effective Rabi-frequency which depends on the details of the two-photon transition (see Appendix B). Eq. (6) excites a collective spin wave of the array with an effective quasi momentum For ω f − ∆ f = 2ω e , we show in Fig. 1.c under which conditions on the angle α and lattice spacing a it is possible to match the driving with the most subradiant single-excitation state, namely K z = q a ≡ π/a. The driving Hamiltonian Eq. (6) produces Rabi oscillations at frequency Ω 0 √ N between the atomic ground state and the collective state |ψ qa ≡ 2/(N + 1) N j=1 sin(z j q a )σ + j |0 . The state |ψ qa approximates the actual dark mode |q a with an overlaperror scaling as ∼ 1/N 2 [11].
To study the fidelity of the subradiant state preparation, we numerically simulate the evolution of the open system under the condition of no jump occurring. Specifically, we calculate the state at time t as |ψ(t) = exp[−i(Ĥ 1array +V 1array )t/ ]|0 , where we chose ω d in resonance with the energy of the target state |q a (we refer to Appendix H for details on the method used for the simulation). The error in the target state preparation is calculated as ≡ 1−max t [F(t)] where F(t) ≡ | q a |ψ(t) | 2 [27]. In Fig. 2, we illustrate the dependence of the dark state preparation's error on both the lattice spacing a/λ e and : error (infidelity) (blue circles), total population in the ground state (P0) and two excitation manifold (P2) at the time which minimizes the error (red squares), error calculated with the effective three-level model Hamiltonian Eq. (E3) (pink diamond), and fit function ε ≈ 0.45 × N −0.6 (black dashed line). d) N -dependence of the optimal frequency Ω opt 0 at which the minimal error occurs for a = λe/4. Scaling N −2.5 as guide for the eye (black-dashed line). the number of atoms N in the array. The data shown in Fig. 2.a,b are obtained optimizing the fidelity with respect to the Rabi-frquency in Eq. (6). Fig. 2.a,b shows the two different contribution of the array's non-linear response: the error decreases for smaller lattice spacing and for larger N as a consequence respectively of the dipole shift ∆ 2 and of enhanced decay of the double excited state as compared to the decay of the target state (Zeno effect). The fidelity F improves for smaller a/λ e at a rate which reduces with N , as expected from the reduction of ∆ 2 for longer arrays (Fig. 2.a). This effect is particularly evident when comparing the case of N = 2 with the other lines, but it is present also for larger array sizes as shown by the crossing of the results for N = 60 and for N = 100 with the results for N = 20 at a/λ e ≈ 0.01. This implies that for a fixed value of a/λ e , increasing N might at first lead to an increased error. Ultimately however the er-ror will improve for larger N at fixed a/λ e (see the case of a/λ e = 0.05 in Fig. 2.b). The improvement in the error with the array's size N is a signature of the Zeno effect. In particular, when N −7/2 Ω 0 /Γ 0 N −3/2 , it is possible to efficiently excite |q a while the population transfer to the doubly excited state |2q a , and hence to higher excited manifolds containing n > 2 excitations, is suppressed by the enhanced decay. The contribution of the Zeno effect to the non-linear response is prominent for larger interatomic separations. In Fig. 2.c, we extend the results in Fig. 2.b for the case of a = λ e /4 to much larger array sizes, showing an improvement os the error for larger N . This improvement is only due to the Zeno effect as proven in Fig. 2.d, where the optimal frequency Ω opt 0 as function of N shows a scaling ∼ 1/N 2.5 falling within the limit of the Zeno regime. At Ω opt 0 , the error has two major contributions coming from (i) the undesired transfer of population to levels other than the target and (ii) the finite decay rate Γ qa of the target state |q a . The latter alone would predict an error scaling ∼ N −1 as Γ qa /(Ω opt 0 √ N ) ∼ 1/N . The main limitation to the state preparation fidelity is thus mainly due to population transferred to the most-subradiant two-excitation state due to imperfect Zeno-blockade ( Fig. 2.c). The dynamics of the driven array within the Zeno regime can thus be captured by a simple three level model involving only the ground state |0 , the target state |q a , and the most subradiant two-excitation state |2 (see Appendix E), as shown by the pink diamond markers in Fig. 2.c.
We define logical qubit states of an array as |0 L ≡ ⊗ N |g 1 and |1 L ≡ |q a = j c qa,j |g 2,j , where c qa,j ≡ g 2,j |q a . Note also that, once the collective dark mode has been prepared, it is possible to turn off the dipole coupling between the levels |g 1 and |g 2 by turning off the laser which couples |g 2 to |e (ε = 0 in Eq. (2) and Eq. (4)). In this way, the dark state becomes a metastable state which can be stored for long in the array. Let us now show how, given two parallel arrays of atoms, one can perform entangling operations between the logical states of the arrays.

B. Two Parallel Arrays: Entangling √ iSWAP-Gate via dipole-dipole interactions
The dynamics of two parallel arrays is described by the Hamiltonian in Eq. (1). As for the case of a single array, the Hamiltonian can be diagonalized separately for each number of excitations. For the case of a single excitation and infinite arrays, it is possible to prove that the dipole interaction between the two arrays, Eq. (4), couples only state with the same quasimomentum and thus we can write the single excitation Hamiltonian of two parallel arrays asĤ on the basis {|k0 , |0k }. Here, |k0 (|0k ) is the state with one excitation of momentum k in the array A (B), and g k (γ k /2) is the coherent (dissipative) coupling between |k0 and |0k . For an infinite array, an exact formula can be derived for g k − iγ k /2 (see Appendix F). In this limit, for the case of atoms polarised along the array's axis, the coupling between subradiant modes (k > k e ) of infinite arrays reads and γ k = 0, where K 0 (x) is the modified Bessel function of the second kind. Eq. (9) holds to a high degree of accuracy also for finite arrays provided N 1. In particular, the coupling is coherent and falls off approximately exponentially with the array separation l as described by Eq. (9). While the decoupling between states of different quasi-momentum is strictly true only for infinitely long arrays (and for the particular case of N = 2), in Fig. 3.a we show that already for arrays as small as N = 6 atoms for l = a = λ e /4 the coupling between different kblocks is negligibly small. Such a strong suppression of the cross coupling between modes associated to different quasi-momenta k can be traced back to the accuracy of approximating the exact system's eigenmodes with quasimomentum eigenstates |k . Indeed, the cross coupling is larger around the light cone separating bright from dark states where this approximation is less accurate [11].
In Eq. (8), we included a detuning δ between the two arrays. For arrays separated by l < λ e , such a selective detuning could be obtained via an AC Stark-Shift. Using a standing wave and placing the array A in a node of the field, only the atoms in the array B pick up a shift δ ( Fig. 3.b). Alternative ways using electrostatic or magnetostatic field gradient could also be envisioned. When the arrays are separated by a distance comparable to λ e , they are collectively driven according to the Hamiltonian (in a frame rotating at ω d ) Hence, the possibility of detuning one array with respect to the other is instrumental for selectively addressing one array. In Fig. 3.c, we show the scaling of the fidelity for preparing the state |10 L ≡ |q a 0 of array A as function of the lattice spacing a/λ e for l = a when array B is detuned by δ = 100Γ 0 . The state preparation error (blue circles in Fig. 3.c) decreases for smaller lattice separation a/λ e owing to the larger dipole shift as for the case of a single array. Accordingly, the population transfer to the two excitation manifold also decreases with a/λ e . For a/λ e 0.08 the error increases if one further reduces the lattice separation because the detuning δ is now comparable to or smaller than the dipole coupling between |10 L and |01 L leading to substantial population transfer to |01 L . This explanation can be confirmed with an effective four-level model obtained by projecting the full Hamiltonian Eq. (2) on the subspace spanned by the levels {|00 L , |10 L , |01 L , |ψ 2 }, where |ψ 2 is the most subradiant two-excitation state of the parallel-arrays system. The error and the population transferred to |01 L as calculated by this simple model are shown in Fig. 3.c by the solid blue and dashed red line respectively. For larger separation l the coupling between |10 L , |01 L decreases leading to better state preparation fidelities at smaller values of a/λ e as compared to the case of l = a shown in Fig. 3.c. We note that by applying this procedure sequentially to both arrays, it is possible to initialise the system into any states of the computational basis {|00 L , |01 L , |10 L , |11 L }.
The structure of the two excitation manifold of two parallel arrays is more complicated. Two-excitation eigenstates of Eq. (1) can still be classified between bright and dark modes depending on their decay rate, however the exact form of such states shows a strong dependence on the separation l between the arrays. For large separation l a, the two arrays are non-interacting and the most-subradiant two excitation state is simply given by |ψ 2 = |11 L ≡ |q a , q a and it is characterised by a decay rate Γ 11 which tends to 2Γ qa in the limit of l/a → ∞. As the arrays are brought closer, the dipole-dipole inter-actionĤ AB becomes stronger thus substantially modifying the form of the most-subradiant two-excitations state (see Appendix F). In particular for l ∼ a, |11 L strongly couples to a large number of states in the two-excitations manifolds, exhibiting a generally complicated dynamics.
Let us now discuss how one can realize an entangling gate between the computational states |00 L , |01 L , |10 L , and |11 L of two parallel arrays. The resonant dipole-dipole interaction Eq. (4) between two parallel arrays naturally implements a √ iSWAP gate on the time scale T g ≡ π/4g qa . The system's dynamics on the computational subspace {|00 L , |10 L , |01 L , |11 L } realises the following mapping Compared to the truth table of the ideal √ iSWAP-gate (Γ = 0 and ξ = 1), the realisation in Eq. (11) shows two main sources of imperfection which come from (i) the error 1 due to finite decay Γ qa of the computational states |10 L and |01 L and (ii) the error due to loss of population from |11 L . At small interatomic separation a/λ e 1, Γ qa T g 1 and the error 1 can be estimated

d) Error
√ iSWAP-gate as function of a/λe for N = 40 atoms. Different markers corresponds to different separation l/a between the arrays as specified by the legend. e) Error √ iSWAP-gate as function of N for a/λe = 0.10. In both panel d) and e), black dashed lines with star markers represent the case case in which the dynamics of the states |10 L , |01 L , and |11 L is assumed to be an exponential decay at rates Γq a and 2Γq a .
For fixed l/a, Eq. (12) decreases both for longer arrays and for smaller lattice spacing a/λ e due respectively to the scaling Γ qa ∼ Γ 0 /N 3 and to the decrease of T g . At large separations l/a, 1 increases according to [K 0 (l/a)] −1 ∼ l/a exp(l/a) as expected for an increased gate time. The second source of errors is due to the loss of population from the state |11 L which we model by a parameter |ξ| < 1 in Eq. (11). The dynamics of |11 L strongly depends on the separation l between the arrays. At large separation (l a), |11 L is an eigenstate of the parallel arrays, and it decays exponentially with a characteristic rate Γ 11 2Γ qa ∼ 2Γ 0 /N 3 . In this regime, we can approximate |ξ| exp(−Γ 11 T g ). For small separation (l ∼ a), |11 L couples strongly to other states in the two excitations manifold leading to large population transfer outside the computational subspace of the parallel arrays.
To esimate the performance of the proposed implementation of the √ iSWAP-gate, we calculate the average gate fidelity according to [28,29] F G ≡ 1 20 where M ≡PÛ † 0ÛP . Here,P is the projector on the computational subspace of the parallel arrays,Û ≡ exp(−iĤ 0 T g ) , andÛ 0 is the ideal gate which acts on the computational subspace according to the truth table in Eq. (11) with ξ = 1 and Γ = 0, and as the identity on all the other states. We numerically evaluate Eq. (13) as a function of the system parameters N, a/λ e , and l/a and show the results for the gate error tot = 1 − F G (infidelity) in Fig. 3.d,e. In Fig. 3.d, we plot the average gate error as a function of the lattice separation a/λ e for N = 40 and different array separations l/a. For large separation l/a, the gate error is well described by the spontaneous decay from the computational states |10 L , |01 L , and |11 L (black dashed lines in Fig. 3.d) as expected from the decoupling of |11 L from other states in the two-excitations manifold. Accordingly, the error improves when reducing a/λ e as a consequence of the reduced gate operation time (see Eq. (9)). In particular, when ΓT g 1 the total error reads tot 3Γ qa T g /5 and scales as in Eq. (12). For the case of arrays of N = 40 atoms as shown in Fig. 3.d, the decoupling of |11 L from the two excitation manifold is found to be accurate already for l/a ≥ 4. For smaller arrays, such a decoupling occurs already for smaller values of l/a as shown in Appendix F. In the opposite regime, i.e. when arrays are placed close to each other, |11 L couples strongly to different states in the two-excitation manifold. This leads to a large population transfer outside the computational subspace and hence to an increased error. This effect is dominant in the case of l/a = 1, 2 in Fig. 3.d. In this case, the strong coupling leads to an oscillation of population between |11 L and |S 2 = (|2q a , 0 + |0, 2q a )/ √ 2 [30]. The case l = 3a in Fig. 3.d, represents an intermediate case where the decoupling of |11 L is not perfect and a small fraction of population is exchanged with |S 2 . For this reason, the curve shows a slight improvement for smaller a/λ e due to a reduction in the gate operation time, which however saturates to a value determined by the population transferred from |11 L to states outside the computational subspace. In Fig. 3.e, we plot the average gate error as a function of N for a/λ e = 0.10 and different array separations l/a. As for the case of Fig. 3.d, for large separations between the arrays, the error is solely due to the spontaneous decay from the computational states |10 L , |01 L , and |11 L as shown by the agreement between the colored markers and the black dashed lines in Fig. 3.e at l/a = 4, 5. At short separation the error is instead again dominated by the population transfer from |11 L to states outside the computational subspace. We observe a trend, particularly evident for the case l/a = 3, which shows a first improvement in the error for short arrays followed by a later increase in the gate error. As we show in Appendix F, this can be interpreted as coming from the reduction of the dipole shift between |11 L and |S 2 with growing array's length, which leads to a larger exchange of population between the two states.
The picture presented here for arrays of pinned atoms has not considered the detrimental effects which come from the atoms' center-of-mass motion around their trapping position in the lattice. As part of the system's nonlinear response comes from the destructive interference of the emitted field, we expect the atomic fluctuation to drastically reduce the intrinsic dissipative non-linearity. In the following section, we discuss the motional effects and their impact on the fidelity of single and two qubit operations.

III. ARRAYS OF FLUCTUATING ATOMS: MOTIONAL EFFECTS
We now consider the fluctuating motion of the array's atoms around their trapping position (see Fig. 4.a). In this case, the open dynamics of an array of fluctuating atoms is described by the non-hermitian Hamiltonian where we modelled the mechanical degrees of freedom as harmonic oscillators of frequency ω T , and m is the atomic mass. The operatorr j (p j ) represents the centerof-mass position (momentum) displacement relative to the trap centre R j ≡ (0, lj y , aj z ) T , where j ≡ (j y , j z ), j y = 0, 1 labels the array A,B, and j z = 1, . . . , N the position within each array. The second line in Eq. (14) represents coherent and dissipative atomic dipole-dipole interactions where now the interaction strength G(r j ,r i ) is an operator which acts on the atomic center-of-mass degrees of freedom and reads The exponential operator appearing as the last factor in Eq. (15) is a result of the Raman transition used to define the effective two level system {|g 1 , |g 2 }. Similarly, in the presence of atomic fluctuations, the driving Hamiltonian reads (see Appendix B) In general, the dynamics described by Eq. (14) and Eq. (16) leads to correlations between the internal and external degrees of freedom, which quickly complicate the simulation of the time evolution of the full system. In two limiting cases, the study of the dynamics can be greatly simplified [12]: (i) the slow atomic-motion regime ω T Γ qa , and (ii) the fast atomic-motion regime ω T Γ 0 . The slow-motion limit describes the situation where the time scale of the center-of-mass dynamics is much longer than the one of the slowest internal dynamics, typically proportional to the decay rate of the subradiant modes. Under this condition, the atoms can be considered frozen at their current positions during the internal evolution, and the system's dynamics can be approximated solely by the internal dynamics where the coupling Eq. (15) is evaluated with the substitutionr j → r j where r j is the particular value of the displacement which determines the instantaneous position of the j-th atom [12]. The average dynamics of the system is thus determined by solving the evolution for different realisations of the atomic position and then averaging the results. Applying this method to the one-and two-qubit gates described in Sec. II, one finds poor fidelities.
The fast-motion limit describes the opposite situation in which the center-of-mass dynamics is much faster than the internal one. While this regime is usually challenging for neutral atoms in optical trap, it might be possible to meet such condition for an appropriate choice of the parameters of the Raman transition defining the effective two level atom, such that ω T /γ 1 ε 2 . In this case, the evolution of the system can still be approximated solely by the internal dynamics described by the non-hermitian Hamiltonian in Eq. (14) (and associated Jump operator) averaged over many different realisations of the atoms' positions [12]. This second limiting case leads to better fidelities for the one-and two-qubit operations between parallel arrays as we shall now prove.

A. Regime of Fast Atomic Motion
In the limit of fast atomic motion, ω T Γ 0 , we approximate Eq. (14) by the following non-hermitian Hamiltonian where we defined the averaged coupling over the atomic positionsG ij according tõ where we assumed the atoms' positions to be independent variables, normally distributed according to the probability distribution P (r) ≡ exp(−r 2 /2σ 2 )/( √ 2πσ) 3 . In the Lamb-Dicke regime σk e 1, one can approximate Eq. (18) as (Appendix C) where G ij is the coupling between pinned atoms as given in Eq. (3) and Eq. (5). Substituting Eq. (19) into Eq. (17), we obtain (in a frame rotating at ω g ) According to Eq. (20), the fast atomic motion renormalizes the coherent and dissipative coupling between the atoms by a factor (1 − 2σ 2 k 2 e ), and adds an independent decay rate ∼ σ 2 k 2 e Γ 0 for each atom. The main effect of the atomic motion is thus to suppress the scaling ∼ Γ 0 /N 3 of the dark mode decay rate, which saturates at the constant value ∼ σ 2 k 2 e Γ 0 as shown in Fig. 4.b for the case of a single array. This effect has been pointed out before in [4,21]. The expression we obtain in Eq. (20) differs from the corresponding ones in [4,21] by a factor of two. This increased noise originates from the scattering of the Raman-laser which adds an additional contribution ∼ k 2 e σ 2 to the atomic center-of-mass diffusion (see Appendix C).
Let us now consider the fidelity for preparing the dark mode of a single atomic array in the regime of fast atomic motion. The driving Hamiltonian Eq. (16), averaged over the atomic fluctuations, reduces to the pinned-atom expression Eq. (6), with a renormalized driving frequencỹ Proceeding as in Sec. II, we simulate the Schödinger evolution generated by the total averaged Hamiltonian (including the driving) and calculate the error for preparing the most subradiant single-excitation state of the system. We calculateG ij by averaging Eq. (15) over one hundred realizations of the atomic positions assumed to be distributed around the lattice sites according to P (r). For the driving Hamiltonian there is no need to average the Rabi frequency over different realizations, because the effect of the motion leads to an overall renormalization factor in the driving, which is irrelevant after optimizing Ω 0 . We present the results of the numerical simulations in Fig. 4.c. As a consequence of the saturation of the dark-mode decay rate ( Fig. 4.b), the state-preparation fidelity decreases at large N for a fixed lattice spacing a/λ e as evidenced in Fig. 4.c, where we consider the case √ 2σ/a = 0.05. For larger values of σ the dependence on N is qualitatively the same, albeit with larger noise, and the small improvement at small N is lost (see cases N = 4, 8, 12 in Fig. 4.c). The results in Fig. 4.b,c suggest that, due to the suppression of the Zeno effect in the presence of atomic fluctuations, it is not advantageous to use larger arrays.
Let us now analyze the effects of fast atomic motion on the fidelity of the single-and two-qubit gates in a system of two parallel arrays. Because the decay rate of the single excitation dark mode is clamped by the fluctuations after N 40 for the case √ 2σ = 0.01a ( Fig. 4.b), we investigate the gate fidelity only for small arrays up to a maximum of N = 40 atoms. The results of the numerical simulation for the gate fidelity of two parallel arrays are shown in Fig. 5. In Fig. 5.a, we compare, for different values of σ, the minimal error for preparing the state |10 L of two parallel arrays of N = 20 atoms separated by a distance l = a as function of the lattice spacing a/λ e . As expected, the error grows with σ as a consequence of the contribution (σk e ) 2 Γ 0 to the darkmode decay rate. This effect is particularly prominent at large values of a/λ e where the main contribution to the non-linearity comes from the Zeno effect. For small a/λ e 0.05, values corresponding to different σ yield similar results because the main contribution to the error comes from populating the state |01 L , as for the case of pinned atoms (see discussion in Sec. II B). In Fig. 5.b, we compare the √ iSWAP-gate error, calculated according to Eq. (13), as a function of N for different value of the separation l between the arrays, assuming a position noise characterized by √ 2σ = 0.10a. The gate fidelity deteriorates for increasing l, much faster than for the case of pinned atoms [Cf. Fig. 3.e]. This behaviour is consistent with the increased decay rate of the array's dark mode, because we expect the free space decay to have a larger contribution to the total error at larger l due to an increased gate time T g . In Fig. 5.c, we study in more detail the particular case of l = 2a, which yields better fidelities than other values of l in Fig. 5.b [31]. In particular, we compare the scaling of the error with both N (left panel) and a/λ e (right panel) for different values of position noise σ. While larger noise typically leads to an increase in the gate error, we find that for N = 12 moderate noise might improve the gate fidelity as compared to the pinned-atom case (see gray arrow in the right panel of Fig. 5.c). The left panel of Fig. 5.c shows that at larger lattice separation positional noise always deteriorate the gate fidelity, as expected for a reduced Zeno blockade. However, for a 0.10λ e , a value of √ 2σ = 0.05a leads an improvement in the gate error over the pinned-atom case. This improvement is attributed to a change in the effective coupling rate between the states |11 L and |S 2 , which leads to a larger population in |11 L at t = T g . We do not further investigate this effect, as we believe the values of the error calculated from the effective model Eq. (17) do not yield correct estimates for a/λ e 0.1 (hatched region in Fig. 5.c) as we argue in the following.

B. Limitation of the current description of the atomic fluctuations and alternative approaches
The treatment of the atomic fluctuations presented in the previous section does not take into account correc- tions due to a finite velocity of the atoms [19,[32][33][34][35]. Finite atomic velocity might lead to excitation of the atomic motional state, which would result in a reduction in the fidelity of the one and two qubit gates between arrays. In the following, we take into account effects of finite atomic velocity using a perturbative treatment of the atomic fluctuations [36]. As in Sec. III A, we consider a subwavelength atomic array (a < λ e /2) where atoms are well localized around their trapping position, namely r 0 a, with r 0 ≡ /2mω T the center-of-mass zero point motion. Under this assumption and for low thermal phonon number n th < 1, the system is in the Lamb-Dicke regime, η ≡ σk e 1, for σ ≡ r 0 √ 2n th + 1 and we approximate the coupling between internal and external degrees of freedom by a power expansion of G(r) to second order in η [21], Here,Ĥ 0 describes the dynamics in the absence of mechanical effects of light and is simply given by Eq. (1) with the addition of the mechanical energy of the atoms. The termĤ I1 (Ĥ I2 ) represents the first (second) order correction in the atomic center-of-mass displacement. They have the general form Here,V 0 is the driving Hamiltonian for arrays of pinned atoms as given in Eq. (10). The termV 1 (V 2 ) represents the first (second) order correction in the Lamb-Dicke parameter. They read (see Appendix D) Here, we assumed the driving laser to be directed orthogonal to the y-axis. Moreover, we defined η x,z ≡ K x,z r 0 , wherex j ≡ r 0 (b † xj +b xj ),ẑ j ≡ r 0 (b † zj +b zj ), and K x ≡ k a sin α − k b sin β, while K z is given in Eq. (7).
We now proceed to study the dynamics generated by Eq. (21) and Eq. (24) in the fast-motion regime ω T Γ 0 for the case of single-and two-qubit gates between parallel arrays. We consider here the case where the atoms are initially prepared in their motional ground state (n th = 0). In particular, we compare the results obtained from exact diagonalisation of the total Hamiltonian Eq. (21) and Eq. (24), and the results obtained from an effective model in which the atomic motion has been adiabatically eliminated keeping contribution up to second order in η (Appendix D). This effective model contains both contributions from processes which do not change the motional state of the atoms (first order inĤ I2 andV 2 ) and from second order processes via intermediate excited states of the atomic motion (second order inĤ I1 andV 1 ). In Fig. 6.a, we show the error for preparing the state |10 L for a system of two parallel arrays of N = 2 atoms separated by a distance l = a and trapped with a center-of-mass frequency ω T /Γ 0 = 100. We assume one array to be detuned from the target one by δ. The results of exact diagonalization and of the effective model show excellent agreement at large lattice spacings a/λ e . For smaller interatomic separations, the difference between the exact and effective model is due the excitation of the center-of-mass motion from atomic recoil (Fig. 6.b).
In Fig. 6.c, we show the dependence on a/λ e of the error in performing a √ iSWAP operation between two arrays of N = 2 atoms separated by a distance l = 2a. We ascribe the decrease of fidelity at small interatomic separations to excitation of the center-of-mass motion. We observe that the values of a/λ e at which the center-ofmass motion is excited depends on how tight the trap is: for smaller r 0 /a (tighter trap) fluctuations are excited for smaller a/λ e . Fig. 6.a-c show that, at small interatomic separation, the main source of error is due to the excitation of center-of-mass fluctuations caused by the terms Eq. (22), Eq. (23), Eq. (25) and Eq. (26). Precisely, the excitation of the atomic fluctuations is predominantly due to the atomic recoil involved in the dipole-dipole interaction [Eq. (22) and Eq. (23)] because these terms represents the largest coupling between internal and external degrees of freedom. The mathematical origin of this large coupling is in the faster divergence of the derivatives of G ij for small separation. Because of the enhancement in the the coupling rate in Eq. (22) and Eq. (23) for small atomic separation, the perturbative approach to the atomic fluctuations developed here is valid only for sufficiently small η such that ηk −1 e |∂ α G(r)/G(r)| 1 and η 2 k −2 e ∂ α ∂ β G(r)/G(r) 1 for α, β = x, y, z. Furthermore the critical values of η for which the perturbative approach is justified depends on N : for fixed η one observes the appearence of unphysical eigenvalues of Eq. (21) with positive imaginary part as N is increased.
The perturbative model, while not sufficiently accurate to correctly describe the effects of motion for long arrays, validates the results of the fast atomic motion regime for a 0.1λ e . Furthermore, it indicates that for smaller interatomic separations the center-of-mass motion is expected to play a major role in the sytem dynamics, a feature not captured by the model in Sec. III A.

IV. DISCUSSION
The results presented in the previous sections are obtained under several assumptions: (i) the possibility to trap atoms in an optical lattice which is subwavelength with respect to the transition |e ↔ |g 1 , (ii) a particular atomic structure which comprises two distinct Raman transitions connecting |g 2 to |g 1 via the intermediate levels |e and |f where ω f ≥ ω e , (iii) a negligible dephasing rate of |g 2 (κ 0 Γ 0 ), (iv) fast atomic centerof-mass fluctuations ω T Γ 0 , and (v) the Lamb-Dicke regime for the atomic motion. Let us now show how these requirements can be met using ultracold alkalineearth atoms in an optical lattice. The relevant level structure for such atoms is shown in Fig. 7 where we marked transitions and levels used to implement the scheme illustrated in Fig. 1. The ground-state manifold levels |g 1,2 are encoded in the levels 3 P 0 and 3 P 2 respectively, while the excited state |e is encoded in 3 D 1 , and the auxiliary state |f used for the excitation of the array's col- lective dark modes is encoded in 3 S 1 . Decay from 3 P 0 to 1 S 0 is a forbidden transition which happens at a rate ≈ 10mHz for Sr, while decay from 3 P 2 to 1 S 0 has an even longer predicted lifetime [37]. The two levels |g 1 and |g 2 can thus be considered stable. Alkaline-earth atoms exhibit long range dipole interactions on the transition 3 P 0 − 3 D 1 , which combined with the possibility of creating optical lattices using transitions from 3 P J to higher excited states allows for the creation of deep-subwavelength arrays. In the case of bosonic Strontium, for instance, the transition wavelength between |e and |g 1 is 2.6 µm which allows to attain a subwavelength array with lattice spacing a/λ e ≈ 0.08 for an optical lattice with wavelength λ opt 400 nm [38]. Additionally, the transition between 3 P 0 − 3 D 1 has a linewidth γ 1 /2π ≈ 290kHz which is broader than the linewidth γ 2 /2π ≈ 10kHz for 3 P 2 − 3 D 1 [39] resulting in a weak effective dephasing κ 0 /Γ 0 ≈ 3 × 10 −3 of the level |g 2 after the elimination of the level |e . The Raman-driving of deep-subwavelenght arrays can be realised with ω f /ω e 3.8 using the rapid transitions between 3 S 1 − 3 P 2 (γ/2π ≈ 45 MHz) and between 3 S 1 − 3 P 0 (γ/2π ≈ 10 MHz) [39]. Finally, let us consider the requirement on the atomic motion. We assumed the atoms to be sufficiently cold initially to be well localized at their optical lattice sites (in Sec. III B we even assumed them in their motional ground state). This condition can be achieved via known cooling schemes for alkali-earth atoms which requires trapping depths of the order ∼ 10 − 10 3 kHz for standard cold-atoms experiments [37]. Note that the recoil energy of the transition |e ↔ |g 1,2 is weaker than for the scattering of photons from the optical lattice, thus ensuring the system is in the Lamb-Dicke regime η = r 0 k e 1 with respect to the relevant dipole transition. For these values of the trap depth, the trap frequency is of order ω T ∼ 100kHz which is comparable to the decay rate γ 1 of the tran-sition |e ↔ |g 1 . Hence, it should be feasible to enter the fast atomic motion regime (ω T /Γ 0 < 1) as soon as Ω/∆ 0.1. Lattices of cold alkaline-earth atoms thus represent a promising platform for the realization of the scheme proposed in Sec. II and Sec. III.
It is important to point out that the dephasing arising from the Raman driving ( Fig. 1.a) imposes, in principle, an upper limit to the length N of an array, which corresponds to the case in which the decay rate of the most subradiant single excitation state of the array equals the dephasing rate κ 0 . For arrays longer than this upper limit, dephasing plays an important role and cannot be neglected. However, we note that incoherent processes caused by dephasing in the Raman scheme could be avoided by using a cycling transition for the |e ↔ |g 1 , and a two-photon driving for the transition |e ↔ |g 2 as proposed in [12]. Additionally, the limit on the array length imposed by the dephasing is generally negligible while considering the effect of the atomic motion which has a much stronger effects on the optical response of long arrays (see Sec. III).
Finally, let us remark that the scheme proposed here can be extended to more than two parallel arrays by generalizing the idea of selective detuning to several parallel arrays. This can be done by borrowing techniques used in super-resolved fluorescence microscopy [40,41]. For instance, the doughnut-shaped Laguerre-Gaussian mode (p, l) = (0, 1) has a dark central region which is not diffraction limited. Illuminating the system with such a beam and placing the target array in the dark spot allows to selectively tune the other arrays out of resonance.

V. CONCLUSIONS
Our findings can be summarised in four main points. First, we have discussed how to coherently excite dark modes of subwavelength arrays using a Raman laser. This techniques represents a novel alternative to optical phase imprinting techniques [42][43][44]. Second, we described how to realise a universal set of gates based on dipole-blockade between qubit states. We found that the intrinsic non-linear response of collective dark modes of arrays of pinned atoms leads to fidelities of 99% for sufficiently large arrays or small interatomic separation. Third, we considered the effects of atomic motion, and showed that its detrimental effects on the non-linear optical response of atomic arrays can be partially mitigated in the fast atomic motion regime. Additionally, the centerof-mass fluctuations place a bound on the size N of the array as well as on the lattice spacing a/λ e . Surpassing this bound by either considering longer arrays or shorter atomic separations does not lead to an improvement in the fidelity of gate operations or worse might even increase their errors. It is worth mentioning that while working in the fast atomic motion limit allows to reduce detrimental motional effects on the array's collective response, it comes at the cost of an increased gate time and of additional diffusion of the atomic center-of-mass motion due to scattering of the driving photons. At our present understanding, this cost is however a necessary one to pay to partially recover the collective non-linear response of subwavelength atomic arrays. It is an interesting question for future investigations whether a clever pulse scheme can be devised such as to produce the desired internal dynamics while disentangling it from the external one after a pulse cycle, as for the case of trapped ions [45,46]. To this aim one could think to use the collective mechanical modes of an array of atoms which arise due to the mechanical forces mediated by the dipoledipole interactions [47,48]. Another interesting direction, which we are currently pursuing, is the development of better theoretical models which would allow to extend the results of Sec. III to larger arrays while accurately taking into account possible excitation of the center-ofmass motion. Finally, we showed that realisation of the proposed scheme for manipulating quantum information with arrays of cold alkaline-earth atoms seems possible in the near future.

Appendix A: Derivation of the Effective Atomic Dynamics
In this Appendix, we derive the effective internal dynamics of the atoms starting from the full Hamiltonian describing the interaction between the atoms and the free-space electromagnetic field within the dipole approximation. We include the atomic center-of-mass motion for the case of harmonically trapped atoms. These general expressions reduce to the one of Sec. II, if the position of the atom is treated as a simple c-number.
The total Hamiltonian of the system when the atoms are driven by an external laser on the transition |g 2 ↔ |e readsĤ tot ≡Ĥ rad +Ĥ at +Ĥ int +Ĥ L (t) (A1) In Eq. (A1), we definedĤ rad (Ĥ at ) the energy of the field (atoms),Ĥ int the atom-field interaction in the dipole approximation, andĤ L (t) the external laser driving. Specifically, these different terms are defined aŝ kε k |e j g 1,j | + g (2) kε k |e j g 2,j | + H.c. , Here, k L (ω L = c|k L |) is the laser wave vector (frequency), ω k ≡ c|k|, and upon introducing the dipole d ν for the transition |e ↔|g ν (ν = 1, 2), we defined the coupling constants where ε L (E L ) is the polarisation (modulus) of the driving electric field and V the quantisation volume. The derivation of Eq. (14) starting from Eq. (A1) is based on the following two steps. First, we adiabatically eliminate the excited state |e and obtain an effective Hamiltonian describing the dynamics of two level systems of levels |g 1 , |g 2 interacting with the electromagnetic field [49][50][51]. The adiabatic elimination can be carried out independently for each atom. The Schrödinger equation for the single-atom state |Ψ(r) ≡ ψ 1 (r)|g 1 + ψ 2 (r)|g 2 + ψ e (r)|e reads (in a frame rotating at the driving frequency) i ψ 2 (r) =Ĥ cm ψ 2 (r) + Ω 2 e ik·r ψ e (r) + k,ε k g (2) kε k e iω L tâ k,ε k e ik·r ψ e (r), (A9) i ψ 1 (r) =Ĥ cm ψ 1 (r) + k,ε k g (1) kε k e i(ωe−∆)tâ kε k e ik·r ψ e (r), where ∆ ≡ ω e − (ω g + ω L ) andĤ cm ≡p 2 /2m + mω 2 Tr 2 /2. We setψ e (r) = 0 in Eq. (A8) and solve for ψ e (r) neglecting the contribution of ofĤ cm under the assumption ∆ ω T . We obtain an effective equation for the dynamics of ψ 1 (r) and ψ 2 (r). Generalising this procedure to the case of N atoms is straightforward, and one obtains that a Hamiltonian which yields the effective equations of motion for ψ 1j (r) and ψ 2j (r) readŝ kε k e i(ωe−∆)tσj 21 + g (2) kε k e iω L tσj where we redefined ω g to include the AC Stark shift Ω 2 /2∆. We stress that such procedure yields a light matter coupling which can be tuned via Ω/∆. Second, we eliminate the photonic degrees of freedom assumed to be in the vaccum state and obtain a Born-Markov master equation for the atomic variables. Master equations for the internal and center-of-mass dynamics of atoms interacting with electromagnetic field have been derived for single or independent atoms in the context of laser cooling [52,53]. In the case of strong dipole interactions, such master equation can be extended to include coherent coupling and interference [54]. Following a similar derivation, the master equation for the effective two-level atom reads Here,Ĥ is the non-hermitian Hamiltonian of the system and readŝ where we defined the operators Taking the continuum limit for k and making use of the identity 2 describes the dipole emission pattern, n is a unit vector which points in the direction of k, and ε k ⊥ k. Before proceeding we need to make an important remark. As it is, Eq. (A16) does not lead to the correct (collective) shift induced by the electromagnetic field on the atoms as originally derived in [55] using other techniques. Comparing Eq. (A16) with the equivalent expressions in [3], which yields results consistent with [55], one notes that the integral on ω k should be extended to the whole real line. The difference between Eq. (A16) and Eq. (15) in [3] originates from the rotating wave approximation we assumed in Eq. (A4) [56,57]. The largest contribution to the frequency integral in Eq. (A16), comes from values of ω k around ω e − ∆ ω e . In the optical domain such frequency is very large as compared to any other frequencies in the system. It is thus justified to extend the limit of integration In a similar way, starting from Eq. (A15) one can show that where we approximated ω e + ∆ ω e and ω e − ω g − ∆ ω e − ω g ≡ ω e2 , and where the components of the free-space electromagnetic Green's function, G 0 αβ ≡ e α · G 0 (r, ω) · e β (for α, β = x, y, z), read The last two terms in Eq. (A12) are the superoperators representing quantum jumps associated respectively to a decay from |g 2 to |g 1 or a dephasing of |g 2 . They read where the integral is carried out on the unit sphere, γ 1 ≡ |d 1 | 2 ω 3 e /(3πε 0 c 3 ), and γ 2 ≡ |d 2 | 2 ω 3 e2 /(3πε 0 c 3 ). Let us now interpret the processes in Eq. (A20) and Eq. (A21). The process in Eq. (A20) describes the cor-related emission of a photon between each pair of atoms i, j associated to the corresponding recoil of the atomic center-of-mass wavefunction. We observe that the atoms undergo two different recoils: one associated with the absorption of a photon from the laser on the transition |e ↔ |g 2 and a second recoil associated with the spontaneous emission of a photon on the transition |e ↔ |g 1 . The first recoil happens always along the same direction fixed byk L ≡ k L /|k L | where |k L | k e2 , while the second occurs around a random direction as prescribed by the dipole emission pattern Ω(n). Similar interpretation holds for Eq. (A21). In the following, we assume γ 2 γ 1 . Within this approximation the term proportional to F (r i ,r j ) in Eq. (A13) and the contribution J 22 (ρ) are negligible and the atomic evolution can be approximated by the non hermitian Hamiltonian and the action of stochastic quantum jump according to Eq. (A20).

Limit of Pinned Atoms
The case of atoms pinned to their lattice site is readily obtained from Eq. (A22) by settingr j = 0. The center-ofmass motion is thus decoupled from the internal motion and the mechanical energy contribution in Eq. (A22) can be neglected. Following this procedure one obtains the non-hermitian Hamiltonian where G ji is simply given by Eq. (A17) where the centerof-mass fluctuations r j have been neglected. The complex coupling rate G ji is not the same as the free-space dipolecoupling rate between two two-level systems at positions R i and R j , because of the phase factor e k L ·(R j −R i ) due to the Raman driving [Cfr. Eq. (A17) and Eq. (22) and Eq.(26) in [3] (beware that we used a different notation as compared to [3])]. However, G ji reduces to the usual form of the free-space dipole coupling for the case of pinned atom when we consider the Raman laser in Eq. (A5) to be directed orthogonal to the plane containing the arrays, i.e. k L e x , for arrays lying in the yz-plane. We assumed this to be the case in all results presented in Sec. II.

Appendix B: Derivation of the Driving Hamiltonian
We now discuss in detail the conditions to drive single excitations dark modes of a single atomic array. The results generalise easily to the case of parallel arrays discussed in the text. We consider an additional highly excited level |f for each atom, and drive a two-photon transition between |g 1 and |g 2 through |f (see Fig. 1.b).
When the two driving lasers are detuned from the intermediate level |f , this latter is negligibly populated during the process and can be adiabatically eliminated, yielding the effective driving Hamiltonian (in a frame rotating at the driving frequency ω d ) where k a (k b ) is the wave-vector associated to the laser driving the transition |f ↔ |g 1 (|f ↔ |g 2 ),σ y j ≡ i(σ j 12 −σ j 21 ) and Ω 0 ≡ Ω a Ω b /2∆ f is the effective Rabi frequency of the two-photon transition. Note that to derive Eq. (B1), we assumed the system to be drive both from the left and from the right, a situation which can be achieved using a mirror to reflect the driving laser as illustrated in Fig. 1.b.

Limit of Pinned Atoms
For the case of atoms pinned to their lattice position, Eq. (B1) reduces tô where we further assumed the array to be aligned along the z-axis, R j = Z j e z , Z j ≡ aj z , and we defined K z ≡ k a cos α − k b cos β as in Eq. (7). Here, Fig. 1.b). For a given interatomic separation a, the values of α and β which allows to tune the driving in Eq. (B2) to match the most subradiant single excitation state, can be obtained by setting K z = π/a. We write and assuming ω f = pω a + ∆ f (p ∈ Q), we solve for β varying a/λ 0 and α. The results are shown in Fig. 2.b for p = 2 demonstrating that a Raman transition through a higher level allows for exciting dark modes for arrays with small a/λ e . Let us finally remark that larger values of p allow to drive subradiant modes of arrays with smaller interatomic separation.

Appendix C: Motional Averaging in the limit of Fast Atomic Fluctuations
In the limit ω T γ 1 , the dynamics of an array of fluctuating atoms on the time scale of the internal dynamics can be approximated by an effective master equation for the sole internal degrees of freedom, where the coupling coefficients have been averaged over the center-of-mass motional state. Specifically, the non hermitian Hamiltonian Eq. (A22) in the fast-atomic-motion limit is approximated asĤ whereG ij ≡ G(r i ,r j ) cm , and the average is taken with respect to the probability distribution of the center of mass position of atoms at sites i and j. As in the main text, we assume the position fluctuation of the atoms to be independently and equally distributed according to In the following, we evaluate the expression for the averaged coupling rateG ij for the case of a single array of atoms polarised along the array direction in the limit of tightly trapped atoms, σ a. We start from Eq. (A16) and take the average over the atoms fluctuation with respect to the distribution Eq. (C2). Using the result we obtain, separating real (J ij ) and imaginary (Γ ij ) parts, For the case of i = j, one hasΓ ii = γ 1 ∀i, and the divrgent Lamb-shiftJ ii . For the case i = j is it convenient to analyze Eq. (C4) and Eq. (C5) separately. Starting from the dissipative coupling rateΓ ij , one proceed by changing to spherical coordinates, from which the radial integral can be immediately performed to yield where we defined k e ≡ k e (cos φ sin θ, sin φ sin θ, cos θ) T , k L ≡ k L e x k e e x , assumed the array to be aligned along the z-axis. In the limit of σk e 1, keeping only terms up to order (k e σ) 2 we obtaiñ where Γ ij is the dissipative coupling rate between two pinned atoms at site i and j as given by the imaginary part of Eq. (3). To evaluate the averaged collective dipole shift Eq. (C5), it is convenient to rewrite it in Cartesian coordinates The integral in Eq. (C8) has been approximated in [4,10] for the case of an infinitely long array. By noticing that e 2k 2 e σ 2J ij is independent on σ one can show that a good approximation forJ ij reads, where in the last passage we assumed k 2 e σ 2 1, and J ij is the coherent coupling rate between two pinned atoms at site i and j given in Eq. (3).
Let us point out that the method of averaging over the atomic fluctuations albeit yielding for the couplingG ij similar results as the ones obtained in [4,10] following a renormalisation procedure, it follows a fundamentally different approach. In particular, in [4,10] the free space coupling is properly renormalized yielding finite results forJ ii . The procedure followed here is not a renormalisa-tion: the coupling is still between point-like atoms, which however are randomly distributed according to Eq. (C2).

Averaged Driving Hamiltonian
Under the assumption of fast atomic motion, we can approximate the driving Hamiltonian Eq. (B1), aŝ (C11) The effect of atomic motion, in the limit ω T γ 1 , is thus a renormalisation of the Rabi frequency Ω 0 .

Perturbative Expansion of the Driving Hamiltonian
Within the Lamb-dicke regime, we can expand the driving modulation in Eq. (B1) up to second order in power of r j asV V 0 +V 1 +V 2 . The three contributions to the driving Hamiltonian read Here, we assumed the two lasers to be directed orthogonal to the y-axis, As illustrated in Fig. 1

Adiabatic elimination of the center-of-mass motion: Effective Hamiltonian for the internal dynamics
In the regime ω T Γ 0 , we can approximate the dynamics of Eq. (D1) and Eq. (D8-D10) by adiabatically eliminating the center-of-mass motion. Combining the expression for the Hamiltonian and the driving up to second order in the Lamb-Dicke parameter we obtainĤ t =Ĥ ≡Ĥ i +V i contains the array and driving Hamiltonian up to order η i (i = 0, 1, 2). Assuming the array to be in its motional ground state, we calculate the effective Hamiltonian for the atoms' internal dynamics aŝ Here, we introduced the projectorsP = 1 int ⊗ |0 0| andQ ≡ 1 int ⊗ n =0 |n n|, where 1 int is the identity operator acting on the internal atomic degrees of freedom and |n is the Fock state containing n phonons in the array. Let us now derive an expression for the different terms on the right-hand side of Eq. (D13). We define,Ĥ eff,0 ≡P (Ĥ 0 +Ĥ I2 )P , which is easily obtained from Eq. (D2) and Eq. (D6), and readŝ Eq. (D14) corresponds to Eq. (D7) for σ = r 0 , and is thus equivalent to the limit of fast moving atom Eq. (17). The second term on the right-hand side of Eq. (D13) is defined asV eff ≡P (V 0 +V 2 )P , and it can be written aŝ whereσ j y ≡ i(σ j 12 −σ j 21 ). The contribution of Eq. (D15) corresponds to an averaging of the driving modulation over the atomic center-of-mass state. The last term in Eq. (D13),Ŵ eff ≡ −PĤ (1) totP , represents the coupling between different collective states mediated by the atomic motion. It readŝ where we defined G α (R ij ) ≡ k −1 e ∂ α G ij . The interpretation of the different contributions inŴ eff is the following. The term Eq. (D16) is an off-diagonal contribution to the dipole coupling in Eq. (D2) for it couples different eigenstates of PĤ 0P within the same n-excitation manifold. The term in Eq. (D17) can be interpreted as a distorsion of the driving modulation, which couples the initial state to states other than the targeted dark mode. The last term Eq. (D18) is correction to second order in the drivingV 1 and thus excites two photons transitions.

Appendix E: Effective Three-Level Dynamics of a Driven Array of Pinned Atoms
In the following, we show that the complex many body dynamics of driven single arrays of pinned atoms is effectively confined to a handful of levels. We do so by deriving an effecting model for the driven dynamics.
Within the blockade regime the dynamics of a single driven array is effectively confined to the three-state manifold containing the ground state |0 and the most subradiant one-and two-excitation state. The remaining one-and two-excitation states are not populated due to the mode matching of the driving, while the population transferred to the three excitation manifold is negligible due to the non-linearity of the array (Dipole blockade and Zeno effect). We note that the array's non-hermitian Hamiltonian Eq. (2) can be diagonalized in terms of right |n µ and left |n µ eigenvectors (n = 1, . . . , N ). The set of left and right eigenvectors, which do not form separately orthogonal basis, satisfy the bi-orthogonality relation where we univocally fix the normalization of the left and right eigenvector by the additional requirement n µ |n µ = 1 ∀n, µ. In the following, it is convenient to define |n 1 ≡ |n (|n 1 ≡ |n ) to represent the mostsubradiant right (left) eigenstate with n excitations. The effective three-states non-hermitian Hamiltonian to second order in the couplingV is obtained as [59] H eff =P (Ĥ 1array +V 1array )P whereĤ 1array is the Hamiltonian of a single array as given in Eq. (2),V 1array is given in Eq. (6), we defined the pro-jectorsP ≡ |0 0| + 1 1 + 2 2 ,Q = 1 −P , and the average complex energy of the three-level manifold E 0 , assuming the driving laser to be resonant with |1 . On the bi-orthogonal basis {|0 , |1 , |2 , |1 , |2 }, the effective three-state non-hermitian Hamiltonian in the frame rotating at the driving frequency thus readŝ Here, we defined the coupling rates Ω nm ≡ n|V 1array |m and we neglected the complex energy shift of the ground state as well as a direct coupling between the gound state and the two excitation state. Both these latter quantities are in fact much smaller than the remaining matrix elements already for arrays of few atoms, and decrease for increasing array size. The diagonal elements in Eq. (E3) already contain the correction due to coupling to higher excitation manifolds. They read where ∆ n,µ is the detuning of the driving from the energy of the µ-th state in the n-excitation manifold ofĤ 1array , and Γ n,µ its decay rate. In Eq. (E4) and Eq. (E5) the indexes µ label states within each n-excitation manifold ordered for increasing values of the decay rate Γ n,µ . We also defined δ n ≡ ∆ n,1 and Γ n ≡ Γ n,1 . The effective three-level model in Eq. (E3) captures the correct dynamics of the system in the Zeno regime as shown by the pink diamond in Fig. 2.c. Furthermore, via the effective description in Eq. (E3) we can understand the discontinuous behaviour as function of N of the total population in the ground state and the subradiant twoexcitations state. We point out that the state obtained by creating two subradiant collective excitations on top of each other in the array, 2/(N + 1) sin(z j k N )σ + j , has a large overlap with the two-excitation most-subradiant state of the array even in the limit N 1 (Fig. 8.a). Population transfer from the single-excitation subradiant state to the two-excitation subradiant state under the action of the driving is prevented by the detuning δ 2 of the latter with respect to the laser driving. The dependence of δ 1,2 on the array's length N is discontinuous as shown in Fig. 8.b. This discontinuity is a consequence of the discreteness of the array and it is reflected by the discontinuous jumps in the total population transferred to |0 and |2 observed in Fig. 2.c.
Before concluding this section, let us remark that the effective Hamiltonian derived according to Eq. (E2) is generally different from the one that would have been obtained, had we performed an adiabatic elimination on the full master equation describing the array's dynamics. The two approaches give consistent results only if the additional correction obtained following the latter (more rigorous) method can be neglected. The good agreement Effective three-level model. a) Overlap | 2 |(Ŝ + ) 2 |0 | 2 as function of N . b) Dependence of δ1 (blue circles) and δ2 (red squares) on the array length N . The arrows point at the sudden jumps in the detuning which corresponds to the jumps in the population shown in Fig. 2.c. In both panels we set a = λe/4. between the results obtained from Eq. (E3) and the simulation of the full Hamiltonian in Eq. (2), justifies a posteriori the simple approach followed here.

Appendix F: Description of Two Parallel Arrays of Pinned Atoms
The non-hermitian Hamiltonian describing the open dynamics (conditioned on no jump occurring) for two parallel arrays of atoms is given in Eq. (1). As for the case of a single array, the hermitian and anti-hermitian parts of Eq. (1) commute with the total number of excitationŝ N e . The Hamiltonian of two parallel arrays can thus be block diagonalized in blocks which act only on a manifold of given number of excitations n, Here,Ĥ (n) 0 is the Hamiltonian acting on the n-excitations manifold. In the following, we look separately at the eigenvalues and eigenvectors of one-and two-excitation block ofĤ 0 .

Single Excitation Structure
We consider first the case of periodic boundary conditions, where the single-excitation blockĤ where the two frequencies ω B,k and ω A,k may differ if the atoms of one array are detuned from the ones of the other array. The expression for (g k − iγ k /2) appearing in Eq. (F3) is given in Table I for different direction of the atomic polarization (assuming all atoms to be polarized parallel to each other). The k-space expression of the couplings can be derived as follows. Let us first define the mixed representation of the free-space electromagnetic Green's function as Substituting Eq. (F4) into and using the identity where Q ∈ {2nπ/a} n∈Z is the reciprocal lattice vector of the atomic array, one obtains Before proceeding to evaluateG 0 αβ (x, y, k + Q) we remark that it is sufficient to keep only the term Q = 0 in Eq. (F7). Higher values of Q represent contribution coming from modes beyond the first Brillouin zone (umklapp processes) which can be neglected since a < λ e /2 and ω e Γ 0 . The coupling of an excitation of quasimomentum k between two parallel arrays within the single excitation manifold thus reads and depends on the polarization direction α, β of the atoms in the arrays. The 2D-fourier transform I. Momentum space expression of the conservative g k and dissipative γ k coupling between two arrays within the single excitation subspace. Here Jn(ρ) (Yn(ρ)) and Kn(ρ) are respectively the n-th order Bessel function of the first (second) kind and the modified Bessel function of the second kind. We defined ρ ≡ l √ k 2 − k 2 e when k > ke and ρ ≡ l √ k 2 e − k 2 when k < ke.
Polarization |k| > ke |k| < ke where p ⊥ ≡ p x e x + p y e y and r ⊥ = xe x + ye y , can be evaluated by expressing the integral in polar coordinate and evaluating the polar integral first. To evaluate the radial integral it is necessary to distinguish between the two cases (i) |k| > k e and (ii) |k| < k e . The first case represent the coupling between collective mode lying outside the light cone and thus its dissipative contribution vanish. The second case corresponds to mode within the light cone and it thus have a strong dissipative contribution (see Table I). In the subradiant part of the spectrum γ k = 0, thus the eigenstates of Eq. (F3) are the usual dressed state with frequency ω ± = δ/2 ± δ 2 + 4g 2 k /2 and decay rate Γ k . Within the light cone γ k and g k are comparable and one needs to diagonalized a 2×2 complex symmetric Hamiltonian. The results of such diagonalization is similar to the previous case upon substituting g k with g k − iγ k /2.
For open boundary conditions, the decoupling between collective excitations with different wave number is not exact and in the subradiant sector the coupling acquires a small dissipative part ( Fig. 9.a). Remarkably for l = a ≤ λ e /4 the block diagonal picture still holds accurately for array as short as N = 6 (see Fig. 3.a in the main text). In Fig. 9.b, we compare the dependence of g k for k = q a on the separation l for different array sizes N with the case N = ∞ in Table I. For small separations the agreement is good, while at larger separation shorter arrays deviates from the ideal case. For finite array size, the dressed state of Eq. (F3) in the subradiant sector have different decay rate due to a non-zero γ k . In the limit γ k g k , δ, the dressed complex energies to first order in γ k read E ± ≡ δ ± g k /Ω k + i(Γ 0 ± γ k g k /Ω k ) where Ω k ≡ δ 2 + 4g 2 k .
FIG. 9. Single excitation coupling rate between dark modes of parallel arrays. Scaling of g k − iγ k /2 for atoms polarized along ez for a = λe/4. a) γq a in function of N for l = a. b) gq a in function of l for different array length N (markers) and for an infinite array (dashed black line).

Two Excitation Structure and Dynamics
Let us now study the two excitation structure of two parallel atomic arrays. Here, we focus in particular on the subradiant sector of the two-excitations manifold. Due to the dipole interaction Eq. (4) between the arrays, we cannot make general statements on the form of the first few dark modes as it was done for the case of an isolated array [11]. In particular, the structure of the subradiant states strongly depends on the separation l between the arrays. In the limit of separation l/a 1, the arrays do not interact. The two-excitations most-subradiant state of the system is thus given by the state |11 L = |q a , q a , which contains one subradiant excitation in each array and is a characterised by a decay rate Γ 11 = 2Γ qa . When the arrays are brought closer they start to interact thus modifying the structure of the ground state. In particular, at small separation l ∼ a, we expect the two excitations to exhibit anti-bunching behaviour between the two arrays due to strong dipole-dipole interaction. In Fig. 10.a, we plot the overlap between |11 L and the two-excitations most-subradiant state of the system |ψ 2 . In Fig. 10.b we plot the probability distribution of |ψ 2 on the sites of each array and show that while for l = a it shows anti-bunching behaviour (bottom panel) at large separation converges towards |11 L (top panel).
The dynamics of |11 L during the √ iSWAP-gate operation depends strongly on the array's separation l and length N . Generally |11 L couples to a large number of two-excitation eigenstates of the system, and the resulting dynamics is quite complicated. To elucidate the evolution of |11 L and gain a qualitative understanding of its dynamics, it is convenient to implement a Lanczos transformation [60]. According to such transformation, we can imagine the dynamics of |11 L as the dynamics of a chain of neighbouring coupled sites (see Fig. 11.a). Each site j is represented by a state |φ j reached by the dynamics such that states corresponding to different sites are orthogonal. In particular, |φ 0 ≡ |11 L . On the basis |φ j the dynamics of the system starting in |φ 0 is which describes the nearest neighbour hopping of excitations at rate t j between the sites j and j + 1, as well as the energy δ j of each site. In the particular case considered here, the state directly coupled to |11 L by Eq. (1) reads |φ 1 = |S 2 ≡ (|2q a , 0 + |0, 2q a )/ √ 2. and the coupling rate is simply given by t 1 = L 11|Ĥ AB |φ 1 . In general, the dynamics of |11 L cannot be reduced to an effective two-level dynamics on the subspace {|11 L , |S 2 }, for the coupling between |S 2 and |φ 3 is typically strong. Such reduction to an effective two level dynamics is justified only at large separation l/a as shown by the agreement between the black dashed markers and the colored markes in Fig. 3.d,e. The two level dynamics represented by the first block in Eq. (F10) gives however a qualitative explanation to the trend observed in Fig. 3.e. In fact, the decrease in the gate fidelity over the array's size N at small separation l/a can be understood as coming from a reduction in the detuning δ 1 (Fig. 11.b). FIG. 13. Numerical details for subradiant state excitation in a driven single array at the optimal Rabi frequency in Fig. 2.d. a) Comparison of the Infidelity obtained with full diagonalization ofH (blue circles), Krylov-Schur diagonalization truncating after the second (red squares) and third (magenta triangle) exictation manifold. b) Total population in the ground state (red squared), two-excitations manifold (blue circles), three-excitations manifold (magenta triangles), and states other than the target in the one-excitation manifold (black stars). c) Left (right) panel: simulation of the driven array infidelity (total population in other manifold) in function of N for different values of M . In all panels we set a = λe/4. |ψ(0) = |0 at t ∼ 1/Ω opt 0 overlaps with the first handful of states ordered for increasing decay rate (see Fig. 13.c). We stress that, as the Ritz spectrum approximate the exact spectrum up to an estimated tolerance at least 10 −10 , the main error in calculating the transition amplitudes is due to the truncation in Eq. (H4). The result displayed in Fig. (2-3) have been calculated with M = 30.
The simulation for the fast motion regime Fig. (4 -5), are obtained using the same truncation as described above but with the Hamiltonian Eq. (17) instead. Due to a reduced array length used for those simulations, there was no need to use the Krylov-Schur method.