\section{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~\cite{Chang2018}. 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~\cite{Dicke1954,Lehmberg1970a,Lehmberg1970b}.
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~\cite{AntezzaPRL2009,*AntezzaPRA2009,BettlesPRL2016,BettlesPRA2016,Shahmoon2017,Rui2020} and to prepare topological edge modes~\cite{Bettles2017,Perczel2017}.

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~\cite{AsenjoGarcia2017PRX} and read it out using superradiant excitations~\cite{Porras2008}. However, for the creation and manipulation of quantum information one also requires a non-linear optical response, \ie the dependence of light-matter 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 non-linearities 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)~\cite{Lukin2001}. Rydberg excitations have also recently been considered for atomc arrays~\cite{Bekenstein2020,MorenoCardoner2021}

In this article, we propose and analyze an alternative 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 non-linear response of subwavelength arrays of atoms in free space~\cite{Masson2020,Williamson2020,Cidrim2020,Bettles2020} 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~\cite{Facchi2002}. 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~\cite{AntezzaPRL2009,Guimond2019}. 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 (\secref{sec:PinnedAtoms}). Within this simplified assumption, we present the main ideas and discuss how to perform single and two-qubit gates between the ground and the single-excitation subradiant state of atomic arrays.
In \secref{sec:Motion}, 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 \secref{sec:Discussion}, 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 \secref{sec:Conclusions}. 
Additional non-essential details are left to the appendices.

\section{Arrays of Pinned Atoms}\label{sec:PinnedAtoms}

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~\figref{fig:Fig1}.a). 
\begin{figure*}
	\includegraphics[width=2\columnwidth]{./Fig1}
	\caption{{\bf System description.} a) Illustration of two parallel arrays with lattice spacing $a$ and separation $l$. Internal level structure: Raman transition leads to effective two-level systems dynamics with damping $\Gamma_0=\varepsilon^2\gamma_1$ and effective dephasing $\kappa_0\equiv \varepsilon^2 \gamma_2$, where $\gamma_{1,2}$ is the spontaneous emission rate for the transition $\ket{e}\leftrightarrow\ket{g_{1,2}}$. b) Scheme for driving a subradiant (dark) mode of a subwavelength array. The detuning $\delta_f$ can be chosen to resonantly drive a collective mode of the array. c) Colored region: values of $\alpha$ and $a/\lambda_e$ for which a solution for $\beta$ in \eqnref{eq:q_effective} exists. Other parameters $\w_a\equiv \w_f-\Delta_f=2\w_e$, $\w_g=0.1\w_e$, and $K=\pi/a$.}\label{fig:Fig1}
\end{figure*}
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 $\ket{e}$ coupled by a dipole-allowed transition to the ground state levels $\ket{g_1}$ and $\ket{g_2}$ (\figref{fig:Fig1}.a). 
An external laser drives the atoms on the $\ket{e}\leftrightarrow\ket{g_2}$ transition with a detuning $\Delta$ and at a rate $\W$. For $\varepsilon\equiv\W/2\Delta\ll1$, the excited state $\ket{e}$ is never populated and the atom behaves as an effective two-level system with excited state $\ket{g_2}$ and ground state $\ket{g_1}$, characterised by an effective decay rate $\Gamma_0\equiv \varepsilon^2 \gamma_1$ and dephasing rate $\kappa_0 \equiv \varepsilon^2 \gamma_2$.
We further assume $\ket{e}\leftrightarrow\ket{g_1}$ to be a much stronger transition than $\ket{e}\leftrightarrow\ket{g_2}$, so that $\kappa_0\ll\Gamma_0$. In the following, we thus neglect the dephsing arising when driving the $\ket{e}\rightarrow\ket{g_2}$ transition and discuss later under which conditions this approximation can be justified (see~\secref{sec:Discussion}).
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 \secref{sec:Motion}~\footnote{This is however not necessary for the case of pinned atoms.}. 


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~\cite{Lehmberg1970a,*Lehmberg1970b} (see also~\appref{app:ME_derivation})
\be\label{eq:Ham_2Arrays_Pinned}
	\Hop_0 = \Hop_\text{A}+\Hop_\text{B} + \Hop_\text{AB}. 
\ee
The non-hermitian Hamiltonian $\Hop_\nu$ describes the dynamics of the atoms of array $\nu=\text{A,B}$, and reads
\be\label{eq:Ham}
\begin{split}
\frac{\Hop_\nu}{\hbar} \equiv  \w_g \sum_j \spl_{\nu j}\smi_{\nu j}+\varepsilon^2\sum_{i,j}\pare{J_{ij}-\im\frac{\Gamma_{ij}}{2}}\spl_{\nu j}\smi_{\nu i},
\end{split}
\ee 
where $\w_g$ is the splitting between $\ket{g_2}$ and $\ket{g_1}$, $\spl_{\nu j} = (\smi_{\nu j})^\dag \equiv \ket{g_{2,j}}_\nu\bra{g_{1,j}}$ and $j=1,\ldots N$ labels the atom at position $\RR_{\text{A}j}=(0,0, j a)^T$ ($\RR_{\text{B}j}=(0,l, j a)^T$) within array A (B).
Here, $J_{ij}$ ($\Gamma_{ij}/2$) is the coherent (dissipative) interaction between two atoms at sites $i$ and $j$ within the same array, and for $i\neq j$ it is given by~\citep{Lehmberg1970a,*Lehmberg1970b}
\be\label{eq:G_coupling}
	J_{ij}-\im\frac{\Gamma_{ij}}{2} = -\mu_0\w_e^2 \dd_{\nu i} \cdot \GG^0(\RR_{\nu i}-\RR_{\nu j},\w_e)\cdot \dd_{\nu j}
\ee
where $\GG^0(\RR,\w_e)$ is the free space electromagnetic Green's tensor of a point dipole, $\w_e$ the frequency of the $\ket{e}\leftrightarrow\ket{g_1}$ transition, and $\dd_i$ is the electric dipole moment of the atom  at site $i$. For $i=j$, we defined $\Gamma_{ii}=\gamma_1$ and we implicitly included the vacuum Lamb-shift $J_{ii}$ into the definition of $\w_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~\appref{app:Perp_Polarization}.
The last term in \eqnref{eq:Ham_2Arrays_Pinned} describes the interaction between the atoms in the two arrays and reads
\be\label{eq:Ham_int}
	\frac{\Hop_\text{AB}}{\hbar}= \varepsilon^2\sum_{i,j=1}^N\pare{g_{ij}-\im \frac{\gamma_{ij}}{2}}\pare{\spl_{\text{A}i}\smi_{\text{B}j}+\smi_{\text{A}i}\spl_{\text{B}j}}
\ee
where $g_{ij}$ ($\gamma_{ij}/2$) is the coherent (dissipative) part of the dipole coupling given by
\be\label{eq:g_coupling}
 g_{ij}-\im\frac{\gamma_{ij}}{2}\equiv -\mu_0\w_e^2 \dd_{\text{A}i} \cdot \GG(\RR_{\text{A}i}-\RR_{\text{B}j},\w_e)\cdot \dd_{\text{B}j}.
 \ee
Note that as a consequence of the Raman transition used to define the effective two level atoms, the dipole-dipole interactions in \eqnref{eq:Ham} and \eqnref{eq:Ham_int} are proportional to the Green's tensor evaluated at $\w_e$. Hence, despite the effective two level system having a characteristic wavelength $\lambda_g\equiv \w_g/2\pi c$, the array exhibit strong collective dipolar effects only for $a<\lambda_e/2$ (subwavelength condition~\cite{Porras2008,Kornovan2016,BettlesPRA2016,AsenjoGarcia2017PRX}).

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.

\subsection{Single Array: Driving Subradiant Excitations}

The dynamics of an isolated array of atoms in free space is described by \eqnref{eq:Ham}. Because we consider only one array, we drop the label $\nu=\text{A,B}$ and refer to the single array Hamiltonian as $\Hop_\text{1array}$ in this section. 
The eigenstates of \eqnref{eq:Ham} are simultaneous eigenstates of $\hat{N}_\text{e} \equiv \sum_j^N \spl_j\smi_j$, as both the hermitian and anti-hermitian part of \eqnref{eq:Ham} commute with $\hat{N}_\text{e}$.
Within the single excitation manifold and when the atomic array is sufficiently long ($N\gg1$), the eigenmodes of \eqnref{eq:Ham} 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~\citep{Mewton2007,Porras2008,Kornovan2016,AsenjoGarcia2017PRX}.
We define the eigenmode with associated quasi-momentum $k$ as $\ket{k}\equiv\hat{\Sigma}^+_k \ket{0}$, where $\ket{0}\equiv  \ket{g_1}^{\otimes N}$ is the ground state of \eqnref{eq:Ham}, and $\hat{\Sigma}_k^\pm\equiv\sum_j c_{k,j} \hat{\sigma}_j^\pm$ with $c_{k,j}\equiv\braket{e_j}{k}$.
For $a<\lambda_e/2$, there exist single excitation eigenstates of \eqnref{eq:Ham} with 
a quasi-momentum $q>\w_e/c\equiv 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 $\sim \Gamma_0/N^3$~\citep{AsenjoGarcia2017PRX,Zhang2019}, and are referred to as collective subradiant modes or dark modes. 
Dark modes also exist in higher $n$-excitation manifolds, provided $n\ll N$, and are similarly characterized by a scaling of the decay rate $\sim\Gamma_0/N^3$. 
Notably, the $n$-excitation states $\ket{n q} \equiv (\hat{\Sigma}_q^+)^n\ket{0}$ are not eigenstates of  \eqnref{eq:Ham}. 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 $\Delta_n \equiv \w_n-n\w_1$, where $\w_n$ ($\w_1$) is the energy of the most-subradiant $n$-excitation (single excitation) state. The dissipative component is represented by the enhanced decay rate of $\ket{n q}$ which scales as $\sim \Gamma_0/N$ to first order in  \eqnref{eq:Ham}~\citep{AsenjoGarcia2017PRX}. 
We will show that $\Delta_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~\footnote{Note that this cannot be simply achieved with a single external laser as subwavelength excitation lie outside of the light-cone of free space electromagnetic modes.}. This condition can be achieved using a detuned Raman transition via an additional excited state $\ket{f}$. As illustrated in~\figref{fig:Fig1}.b, we consider two driving lasers with wave vectors $\kk_a,\kk_b$ forming respectively angles $\alpha,\beta$ with the array's direction ($z$-axis). For a sufficiently large detuning $\Delta_f$ such that $\ket{f}$ is never populated, the effect of the Raman lasers on the array can be modelled by an effective driving Hamiltonian
\be\label{eq:Vop}
	\Vop_\text{1array} \equiv \im\hbar \W_0 \sum_{j=1}^N \sin(K_z z_j)\pare{e^{-\im\w_d t} \spl_j -\hc},
\ee
Here, $z_j\equiv aj$, $\w_d \equiv \w_a-\w_b$, and $\W_0$ is the single atom effective Rabi-frequency which depends on the detail of the two-photon transition (see~\appref{app:Driving}). 
\eqnref{eq:Vop} excites a collective spin wave of the array with an effective quasi momentum
\be\label{eq:q_effective}
	K_z \equiv k_a \cos\alpha  -k_b\cos\beta.
\ee
For $\w_f-\Delta_f = 2\w_e$, we show in \figref{fig:Fig1}.c under which conditions on the angle $\alpha$ and lattice spacing $a$ it is possible to match the driving with the most subradiant single-excitation state, namely $K_z=q_a\equiv\pi/a$.
The driving Hamiltonian \eqnref{eq:Vop} produces Rabi oscillations at frequency $\W_0\sqrt{N}$ between the atomic ground state and the collective state $\ket{\psi_{q_a}} \equiv \sqrt{2/(N+1)} \sum_{j=1}^N \sin(z_j q_a)\spl_j \ket{0}$. The state $\ket{\psi_{q_a}}$ approximates the actual dark mode $\ket{q_a}$ with an overlap-error scaling as $\sim 1/N^2$~\cite{AsenjoGarcia2017PRX}.

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 $\ket{\psi(t)}=\exp[-\im (\Hop_\text{1array}+V_\text{1array}) t/\hbar]\ket{0}$, where we chose $\w_d$ in resonance with the energy of the target state $\ket{q_a}$ (we refer to \appref{app:Numerics} for details on the method used for the simulation).
The error in the target state preparation is calculated as $\epsilon \equiv 1-\text{max}_t[\mathcal{F}(t)]$ where $\mathcal{F}(t)\equiv |\braket{q_a}{\psi(t)}|^2$~\footnote{This gives an upper bound to the Fidelity as given by the full master equation}.
In \figref{Fig:Fig2}, we illustrate the dependence of the dark state  preparation's error on both the lattice spacing $a/\lambda_e$ and the number of atoms $N$ in the array.
\begin{figure}
	\includegraphics[width=\columnwidth]{./Fig2.png}
	\caption{{\bf Single Array's Dark Mode Preparation.} a) Error for dark state preparation as function of the lattice spacing $a/\lambda_e$ (for different array's length $N$) and b) as function of the number of atoms $N$ (for different values of $a/\lambda_e$). c) Scaling with $N$ for $a=\lambda_e/4$ of different relevant quantities at the optimal driving frequency $\W_0^\text{(opt)}$: error $\epsilon$ (infidelity) (blue circles), total population in the ground state ($P_0$) and two excitation manifold ($P_2$) at the time which minimizes the error (red squares), error calculated with the effective three-level model Hamiltonian \eqnref{eq:Heff_app}  (pink diamond), and fit function $\varepsilon\approx 0.45\times N^{-0.6}$ (black dashed line). d) $N$-dependence of the optimal frequency $\W_0^\text{opt}$ at which the minimal error occurs for $a=\lambda_e/4$. Scaling $N^{-2.5}$ as guide for the eye (black-dashed line).}\label{Fig:Fig2}
\end{figure}
The data shown in \figref{Fig:Fig2}.a,b are obtained optimizing the fidelity with respect to the Rabi-frquency in \eqnref{eq:Vop}.
\figref{Fig:Fig2}.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 $\Delta_2$ and of enhanced decay of the double excited state as compared to the decay of the target state (Zeno effect).
The fidelity $\mathcal{F}$ improves for smaller $a/\lambda_e$ at a rate which reduces with $N$, as expected from the reduction of $\Delta_2$ for longer arrays (\figref{Fig:Fig2}.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/\lambda_e\approx0.01$. This implies that for a fixed value of $a/\lambda_e$, increasing $N$ might at first lead to an increased error. Ultimately however the error will improve for larger $N$ at fixed $a/\lambda_e$ (see the case of $a/\lambda_e=0.05$ in \figref{Fig:Fig2}.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}\ll \W_0/\Gamma_0 \ll N^{-3/2}$, it is possible to efficiently excite $\ket{q_a}$ while the population transfer to the doubly excited state $\ket{2 q_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 \figref{Fig:Fig2}.c, we extend the results in \figref{Fig:Fig2}.b for the case of $a=\lambda_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 \figref{Fig:Fig2}.d, where the optimal frequency $\W_0^\text{opt}$ as function of $N$ shows a scaling $\sim 1/N^{2.5}$ falling within the limit of the Zeno regime.
At $\W_0^\text{opt}$, the error $\epsilon$ has two major contributions coming from (i) the undesired transfer of population to levels other than the target and (ii) the finite decay rate $\Gamma_{q_a}$ of the target state $\ket{q_a}$. The latter alone would predict an error scaling $\epsilon\sim N^{-1}$ as $\Gamma_{q_a}/(\W_0^\text{opt}\sqrt{N}) \sim 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 (\figref{Fig:Fig2}.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 $\ket{0}$, the target state $\ket{q_a}$, and the most subradiant two-excitation state $\ket{2}$ (see~\appref{app:Heff_3levels}), as shown by the pink diamond markers in \figref{Fig:Fig2}.c.

We define logical qubit states of an array as $\ket{0}_\text{L}\equiv \otimes^N \ket{g_1}$ and $\ket{1}_L \equiv \ket{q_a}= \sum_j c_{q_a,j}\ket{g_{2,j}}$, where $c_{q_a,j}\equiv \braket{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 $\ket{g_1}$ and $\ket{g_2}$ by turning off the laser which couples $\ket{g_2}$ to $\ket{e}$ ($\varepsilon=0$ in \eqnref{eq:Ham} and \eqnref{eq:Ham_int}). 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.

\subsection{Two Parallel Arrays: Entangling $\sqrt{\text{iSWAP}}$-Gate via dipole-dipole interactions}\label{ssec:PinnedAtoms_Parallel}

The dynamics of two parallel arrays is described by the Hamiltonian in \eqnref{eq:Ham_2Arrays_Pinned}.
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, \eqnref{eq:Ham_int}, couples only state with the same quasi-momentum and thus we can write the single excitation Hamiltonian of two parallel arrays as $\Hop_0^{(1)}=\sum_k \Hop_k$, where $k=\pi/[a(N+1)]\ldots \pi N/[a(N+1)] $ and
\be\label{eq:H_k_block}
\Hop_k\equiv
\begin{pmatrix}
 \w_k & g_k-\im\gamma_k/2\\
 g_k-\im\gamma_k/2 & \w_k+\delta
\end{pmatrix},
\ee
on the basis $\{\ket{k 0},\ket{0k}\}$. Here, $\ket{k0}$ ($\ket{0k}$) is the state with one excitation of momentum $k$ in the array A (B), and $g_k$ ($\gamma_k/2$) is the coherent (dissipative) coupling between $\ket{k0}$ and $\ket{0k}$. 
For an infinite array, an exact formula can be derived for $g_k-\im\gamma_k/2$ (see~\appref{app:Parallel_Arrays}). 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
\be\label{eq:g_k}
	g_k =-\frac{3\Gamma_0}{k_e a}\pare{1-\frac{k^2}{k_e^2}}K_0\pare{l\sqrt{k^2-k_e^2}},
\ee
and $\gamma_k=0$, where $K_0(x)$ is the modified Bessel function of the second kind.
\eqnref{eq:g_k} holds to a high degree of accuracy also for finite arrays provided $N\gg 1$. In particular, the coupling is coherent and falls off approximately exponentially with the array separation $l$ as described by \eqnref{eq:g_k}.
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 \figref{fig:Fig3}.a we show that already for arrays as small as $N=6$ atoms for $l=a=\lambda_e/4$ the coupling between different $k$-blocks is negligibly small.
\begin{figure*}
	\includegraphics[width=2\columnwidth]{./Fig3}
	\caption{{\bf Parallel Arrays of Pinned Atoms.} a) Left (right) Panel: absolute value of the imaginary (real) part of the matrix elements of $\Hop_2^{(1)}$ on the basis $\{\ket{0k},\ket{k0}\}_k$, namely $H_{kk'}\equiv \bra{\phi_k}\Hop_2^{(1)}\ket{\phi_{k'}}$ where $\ket{\phi_k}
	\in\{\ket{0,k},\ket{k,0}\}$. Diagonal $2\times2$ blocks corresponds to $\Hop_k$, while of diagonal elements are cross coupling between states of different momentum. The light line is marked by red dashed lines dividing subradiant and superradiant sectors. b) Scheme for selective AC Stark Shift to detune one array with respect to the other. c) Error for preparing the state $\ket{q_a 0}$ as function of $a/\lambda_e$ for two parallel array of $N=40$ atoms. Other parameters: $l/a=1$ and $\delta/\Gamma_0=100$. d) Error $\sqrt{\text{iSWAP}}$-gate as function of $a/\lambda_e$ for $N=40$ atoms. Different markers corresponds to different separation $l/a$ between the arrays as specified by the legend. e) Error $\sqrt{\text{iSWAP}}$-gate as function of $N$ for $a/\lambda_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 $\ket{10}_\text{L},\ket{01}_\text{L}$, and $\ket{11}_\text{L}$ is assumed to be an exponential decay at rates $\Gamma_{q_a}$ and $2\Gamma_{q_a}$.}\label{fig:Fig3}
\end{figure*}
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 quasi-momentum eigenstates $\ket{k}$. Indeed, the cross coupling is larger around the light cone separating bright from dark states where this approximation is less accurate~\cite{AsenjoGarcia2017PRX}.

In \eqnref{eq:H_k_block}, we included a detuning $\delta$ between the two arrays. For arrays separated by $l<\lambda_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 $\delta$ (\figref{fig:Fig3}.b). Alternative ways using electrostatic or magnetostatic field gradient could also be envisioned. 
When the arrays are separated by a distance comparable to $\lambda_e$, they are collectively driven according to the Hamiltonian (in a frame rotating at $\w_d$)
\be\label{eq:Vop_2array}
	\Vop_0 \equiv \im\hbar\W_0 \sum_{j=1}^N \sin(K_z z_j)\pare{ \sy_{\text{A}j} +\sy_{\text{B}j}}.
\ee
Hence, the possibility of detuning one array with respect to the other is instrumental for selectively addressing one array.
In \figref{fig:Fig3}.c, we show the scaling of the fidelity for preparing the state $\ket{10}_\text{L}\equiv\ket{q_a 0}$ of array A as function of the lattice spacing $a/\lambda_e$ for $l=a$ when array B is detuned by $\delta=100\Gamma_0$. The state preparation error (blue circles in \figref{fig:Fig3}.c) decreases for smaller lattice separation $a/\lambda_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/\lambda_e$.  For $a/\lambda_e\lesssim 0.08$ the error increases if one further reduces the lattice separation because the detuning $\delta$ is now comparable to or smaller than the dipole coupling between $\ket{10}_\text{L}$ and $\ket{01}_\text{L}$ leading to substantial population transfer to $\ket{01}_\text{L}$.
This explanation can be confirmed with an effective four-level model obtained by projecting the full Hamiltonian \eqnref{eq:Ham} on the subspace spanned by the levels $\{\ket{00}_\text{L},\ket{10}_\text{L},\ket{01}_\text{L},\ket{\psi_2}\}$, where $\ket{\psi_2}$ is the most subradiant two-excitation state of the parallel-arrays system. The error and the population transferred to $\ket{01}_\text{L}$ as calculated by this simple model are shown in \figref{fig:Fig3}.c by the solid blue and dashed red line respectively.
For larger separation $l$ the coupling between $\ket{10}_\text{L},\ket{01}_\text{L}$ decreases leading to better state preparation fidelities at smaller values of $a/\lambda_e$ as compared to the case of $l=a$ shown in~\figref{fig:Fig3}.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 $\{\ket{00}_\text{L},\ket{01}_\text{L},\ket{10}_\text{L},\ket{11}_\text{L}\}$.

The structure of the two excitation manifold of two parallel arrays is more complicated. Two-excitation eigenstates of \eqnref{eq:Ham_2Arrays_Pinned} 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\gg a$, the two arrays are non-interacting and the most-subradiant two excitation state is simply given by $\ket{\psi_2}=\ket{11}_\text{L}\equiv\ket{q_a,q_a}$ and it is characterised by a decay rate $\Gamma_{11}$ which tends to $2\Gamma_{q_a}$ in the limit of $l/a\rightarrow \infty$. As the arrays are brought closer, the dipole-dipole interaction $\Hop_\text{AB}$ becomes stronger thus substantially modifying the form of the most-subradiant two-excitations state (see~\appref{app:Parallel_Arrays}). In particular for $l\sim a$, $\ket{11}_\text{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 $\ket{00}_\text{L}$, $\ket{01}_\text{L}$, $\ket{10}_\text{L}$, and $\ket{11}_\text{L}$ of two parallel arrays.
The resonant dipole-dipole interaction \eqnref{eq:Ham_int} between two parallel arrays naturally implements a $\sqrt{\text{iSWAP}}$ gate on the time scale $T_\text{g} \equiv \pi /4g_{q_a}$.
The system's dynamics on the computational subspace $\{\ket{00}_\text{L},\ket{10}_\text{L},\ket{01}_\text{L},\ket{11}_\text{L}\}$ realises the following mapping
\be\label{eq:TT_real_sqSWAP}
\begin{split}
\ket{00}_\text{L} \rightarrow&\, \ket{00}_\text{L},\\
\ket{10}_\text{L} \rightarrow&\, e^{-\Gamma_{q_a} T_\text{g}}\pare{\ket{10}_\text{L}-\im \ket{01}_\text{L}}/\sqrt{2},\\
\ket{01}_\text{L} \rightarrow&\, e^{-\Gamma_{q_a} T_\text{g}}\pare{\ket{01}_\text{L}-\im \ket{10}_\text{L}}/\sqrt{2},\\
\ket{11}_\text{L} \rightarrow&\, \xi \ket{11}_\text{L}.
\end{split}
\ee
Compared to the truth table of the ideal $\sqrt{\text{iSWAP}}$-gate ($\Gamma=0$ and $\xi=1$), the realisation in \eqnref{eq:TT_real_sqSWAP} shows two main sources of imperfection which come from (i) the error $\epsilon_1$ due to finite decay $\Gamma_{q_a}$ of the computational states $\ket{10}_\text{L}$ and $\ket{01}_\text{L}$ and (ii) the error due to loss of population from $\ket{11}_\text{L}$.
At small interatomic separation $a/\lambda_e\ll1$, $\Gamma_{q_a} T_\text{g}\ll1$ and the error $\epsilon_1$ can be estimated as 
\be\label{eq:error_1}
	\epsilon_1 \simeq \Gamma T_\text{g} \sim \pare{\frac{a/\lambda_e}{N}}^3 \spare{K_0(\pi l/a)}^{-1}.
\ee
For fixed $l/a$, \eqnref{eq:error_1} decreases both for longer arrays and for smaller lattice spacing $a/\lambda_e$ due respectively to the scaling $\Gamma_{q_a}\sim \Gamma_0/N^3$ and to the decrease of $T_\text{g}$. 
At large separations $l/a$, $\epsilon_1$ increases according to $[K_0(l/a)]^{-1} \sim \sqrt{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 $\ket{11}_\text{L}$ which we model by a parameter $|\xi|<1$ in \eqnref{eq:TT_real_sqSWAP}. The dynamics of $\ket{11}_\text{L}$ strongly depends on the separation $l$ between the arrays.
At large separation ($l\gg a$), $\ket{11}_\text{L}$ is an eigenstate of the parallel arrays, and it decays exponentially with a characteristic rate $\Gamma_{11}\simeq 2\Gamma_{q_a} \sim 2\Gamma_0/N^3$. In this regime, we can approximate $|\xi|\simeq \exp(-\Gamma_{11} T_\text{g})$.
For small separation ($l\sim a$), $\ket{11}_\text{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 $\sqrt{\text{iSWAP}}$-gate, we calculate the average gate fidelity according to~\cite{Nielsen2002,Pedersen2007}
\be\label{eq:F_swap}
	\mathcal{F}_\text{G}\equiv \frac{1}{20}\pare{\text{Tr}[\mathcal{M}\mathcal{M}^\dag]+|\text{Tr}[\mathcal{M}]|^2},
\ee
where $\mathcal{M}\equiv\Pop \Udop_0 \Uop\Pop$. Here, $\Pop$ is the projector on the computational subspace of the parallel arrays, $\Uop\equiv \exp(-\im\Hop_0 T_\text{g})$ , and $\Uop_0$ is the ideal gate which acts on the computational subspace according to the truth table in \eqnref{eq:TT_real_sqSWAP} with $\xi=1$ and $\Gamma=0$, and as the identity on all the other states.
We numerically evaluate \eqnref{eq:F_swap} as a function of the system parameters $N,a/\lambda_e,$ and $l/a$ and show the results for the gate error $\epsilon_\text{tot}=1-\mathcal{F}_\text{G}$ (infidelity) in \figref{fig:Fig3}.d,e.
In \figref{fig:Fig3}.d, we plot the average gate error as a function of the lattice separation $a/\lambda_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 $\ket{10}_\text{L},\ket{01}_\text{L}$, and $\ket{11}_\text{L}$ (black dashed lines in~\figref{fig:Fig3}.d) as expected from the decoupling of $\ket{11}_\text{L}$ from other states in the two-excitations manifold. 
Accordingly, the error improves when reducing $a/\lambda_e$ as a consequence of the reduced gate operation time (see~\eqnref{eq:g_k}). In particular, when $\Gamma T_\text{g}\ll1$ the total error reads $\epsilon_\text{tot}\simeq 3 \Gamma_{q_a} T_\text{g}/5 $ and scales as in \eqnref{eq:error_1}. For the case of arrays of $N=40$ atoms as shown in~\figref{fig:Fig3}.d, the decoupling of $\ket{11}_\text{L}$ from the two excitation manifold is found to be accurate already for $l/a\geq 4$. For smaller arrays, such a decoupling occurs already for smaller values of $l/a$ as shown in \appref{app:Parallel_Arrays}.
In the opposite regime, \ie when arrays are placed close to each other, $\ket{11}_\text{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~\figref{fig:Fig3}.d. In this case, the strong coupling leads to an oscillation of population between $\ket{11}_\text{L}$ and $\ket{S_2}= (\ket{2q_a,0}+\ket{0,2q_a})/\sqrt{2}$~\footnote{We found the coupling strength between $\ket{11}_\text{L}$ and $\ket{S_2}$ to be $g_2\simeq 2g_{q_a}$ for all the values $a/\lambda_e$ in the range shown. This is particularly dramatic as at the  $T_\text{g}$ the population in $\ket{S_2}$ is maximised.}. 
The case $l=3a$ in~\figref{fig:Fig3}.d, represents an intermediate case where the decoupling of $\ket{11}_\text{L}$ is not perfect and a small fraction of population is exchanged with $\ket{S_2}$. For this reason, the curve shows a slight improvement for smaller $a/\lambda_e$ due to a reduction in the gate operation time, which however saturates to a value determined by the population transferred from $\ket{11}_\text{L}$ to states outside the computational subspace.
In \figref{fig:Fig3}.e, we plot the average gate error as a function of $N$ for $a/\lambda_e=0.10$ and different array separations $l/a$. 
As for the case of \figref{fig:Fig3}.d, for large separations between the arrays, the error is solely due to the spontaneous decay from the computational states $\ket{10}_\text{L},\ket{01}_\text{L}$, and $\ket{11}_\text{L}$ as shown by the agreement between the colored markers and the black dashed lines in~\figref{fig:Fig3}.e at $l/a=4,5$.
At short separation the error is instead again dominated by the population transfer from $\ket{11}_\text{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 \appref{app:Parallel_Arrays}, this can be interpreted as coming from the reduction of the dipole shift between $\ket{11}_\text{L}$ and $\ket{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 non-linear 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.

\section{Arrays of Fluctuating Atoms: Motional Effects}\label{sec:Motion}

We now consider the fluctuating motion of the array's atoms around their trapping position (see~\figref{Fig:Fig4}.a).
\begin{figure}
	\includegraphics[width=\columnwidth]{./Fig4}
	\caption{{\bf Arrays of moving atoms.} a) Atoms localized around their trapping positions $\RR_\bold{i}$ with a gaussian probability distribution of width $\sigma<a$. b) Decay rate $\Gamma_{q_a}$ of the most subradiant mode for a single array as function of the number $N$ of atoms in the array as calculated by diagonalising the averaged Hamiltonian over $100$ realisation of the atomic positions (colored markers). Horizontal black dashed line corresponds to $(\sigma k_e)^2\Gamma_0$. c) Dependence on $N$ of the optimal error for preparing the most subradiant single excitation state of a single array for different values of $a$ (colored markers) and $\sqrt{2}\sigma=0.05a$.}\label{Fig:Fig4}
\end{figure}
In this case, the open dynamics of an array of fluctuating atoms is described by the non-hermitian Hamiltonian (see~\appref{app:ME_derivation})
\be\label{eq:Ham_motion}
\begin{split}
\Hop \equiv &  \sum_\jj \bigg[\frac{\hat{\pp}_\jj^2}{2m} +\frac{1}{2}m\w_\text{T}^2\pare{ \hat{\rr}_\jj-\RR_\jj}^2+\hbar\w_g\hat{\sigma}_{21}^\jj\hat{\sigma}_{12}^\jj\bigg]\\
& + \hbar\varepsilon^2\sum_{\bold{i},\jj}G(\hat{\rr}_\jj,\hat{\rr}_\bold{i})\hat{\sigma}_{21}^\jj\hat{\sigma}_{12}^\bold{i},
\end{split}
\ee
where we modelled the mechanical degrees of freedom as harmonic oscillators of frequency $\w_\text{T}$, and $m$ is the atomic mass.
The operator $\hat{\rr}_\jj$ ($\hat{\pp}_\jj$) represents the center-of-mass position (momentum) displacement relative to the trap centre $\RR_\jj\equiv(0,lj_y,aj_z)^T$, where $\jj\equiv(j_y,j_z)$, $j_y=0,1$ labels the array A,B, and $j_z=1,\ldots,N$ the position within each array. 
The second line in \eqnref{eq:Ham_motion} represents coherent and dissipative atomic dipole-dipole interactions where now the interaction strength $G(\hat{\rr}_\jj,\hat{\rr}_\bold{i})$ is an operator which acts on the atomic center-of-mass degrees of freedom and reads
\be\label{eq:G_ij_op}
\begin{split}
	G(\hat{\rr}_\jj,\hat{\rr}_\bold{i})\equiv& -\im \frac{\mu_0 \w_e^2}{\hbar}|d|^2G^0_{zz}(\RR_\jj+\hat{\rr}_\jj-\RR_\bold{i}-\hat{\rr}_\bold{i},\w_e)\\
	&\times e^{\kk_\text{L}\cdot(\hat{\rr}_\jj -\hat{\rr}_\bold{i})}.
\end{split}
\ee
The exponential operator appearing as the last factor in \eqnref{eq:G_ij_op} is a result of the Raman transition used to define the effective two level system $\{\ket{g_1},\ket{g_2}\}$. 
Similarly, in the presence of atomic fluctuations, the driving Hamiltonian reads (see~\appref{app:Driving})
\be\label{eq:V_fluct}
	\Vop = \hbar\W_0 \sum_{\jj} \sin[\pare{\kk_a-\kk_b}\cdot (\RR_\jj +\hat{\rr}_\jj)] \sy_\jj.
\ee
In general, the dynamics described by \eqnref{eq:Ham_motion} and \eqnref{eq:V_fluct} 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~\cite{Porras2008}: (i) the slow atomic-motion regime $\w_\text{T}\ll \Gamma_{q_a}$, and (ii) the fast atomic-motion regime $\w_\text{T}\gg \Gamma_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 \eqnref{eq:G_ij_op} is evaluated with the substitution $\hat{\rr}_\jj \rightarrow \rr_\jj$ where $\rr_\jj$  is the particular value of the displacement which determines the instantaneous position of the $\jj$-th atom~\cite{Porras2008}.
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 \secref{sec:PinnedAtoms}, 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 $\w_\text{T}/\gamma_1 \gg \varepsilon^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 \eqnref{eq:Ham_motion} (and associated Jump operator) averaged over many different realisations of the atoms' positions~\cite{Porras2008}. This second limiting case leads to better fidelities for the one- and two-qubit operations between parallel arrays as we shall now prove.

\subsection{Regime of Fast Atomic Motion}\label{sec:FastMotion}

In the limit of fast atomic motion, $\w_\text{T}\gg \Gamma_0$, we approximate \eqnref{eq:Ham_motion} by the following non-hermitian Hamiltonian
\be\label{eq:Ham_fam}
	\Hop_\text{fam} \equiv \hbar \sum_{\bold{i},\jj} \pare{\w_g\delta_{\bold{i}\jj}+\varepsilon^2\tilde{G}_{\bold{i}\jj} } \hat{\sigma}_{21}^\bold{i}\hat{\sigma}_{12}^\jj
\ee
where we defined the averaged coupling over the atomic positions $\tilde{G}_{\bold{i}\jj}$  according to
\be\label{eq:avg_G_ij}
	\tilde{G}_{\bold{i}\jj}\equiv \int_{\mathbb{R}^3}\!\!\text{d}\rr_\bold{i} \,\text{d}\rr_\jj G(\rr_\bold{i},\rr_\jj) P(\rr_\bold{i})P(\rr_\jj),
\ee
where we assumed the atoms' positions to be independent variables, normally distributed according to the probability distribution 
$P(\rr) \equiv \exp(-\rr^2/2\sigma^2)/(\sqrt{2\pi}\sigma)^3$.
In the Lamb-Dicke regime $\sigma k_e\ll 1$, one can approximate \eqnref{eq:avg_G_ij} as (\appref{app:Fast_Atomic_Fluctuations})
\be\label{eq:approx_Gij_avg}
	\tilde{G}_{\bold{i}\jj} = \left \{ 
	\begin{array}{lcr}
	-\im \gamma_1/2 & \text{for} & \bold{i}=\jj\\
	(1-2\sigma^2k_e^2) G_{\bold{i}\jj} & \text{for} & \bold{i}\neq\jj
	\end{array}
	\right . ,
\ee
where $G_{\bold{i}\jj}$ is the coupling between pinned atoms as given in \eqnref{eq:G_coupling} and \eqnref{eq:g_coupling}. 
Substituting \eqnref{eq:approx_Gij_avg} into \eqnref{eq:Ham_fam}, we obtain (in a frame rotating at $\w_g$)
\be\label{eq:Ham_fam_approx}
\begin{split}
	\Hop_\text{fam} =&(1-2\sigma^2k_e^2)\hbar\sum_{\bold{i},\jj} \varepsilon^2 G_{\bold{i}\jj} \hat{\sigma}_{21}^\bold{i}\hat{\sigma}_{12}^\jj\\
	&   -\im \sigma^2k_e^2 \hbar\Gamma_0 \sum_\bold{i}\hat{\sigma}_{21}^\bold{i}\hat{\sigma}_{12}^\bold{i}.
\end{split}
\ee
According to \eqnref{eq:Ham_fam_approx}, the fast atomic motion renormalizes the coherent and dissipative coupling between the atoms by a factor $(1-2\sigma^2k_e^2)$, and adds an independent decay rate $\sim \sigma^2k_e^2\Gamma_0$ for each atom. The main effect of the atomic motion is thus to suppress the scaling $\sim \Gamma_0/N^3$ of the dark mode decay rate, which saturates at the constant value $\sim \sigma^2k_e^2\Gamma_0$ as shown in~\figref{Fig:Fig4}.b for the case of a single array.
This effect has been pointed out before in~\cite{AntezzaPRL2009,Guimond2019}. The expression we obtain in \eqnref{eq:Ham_fam_approx} differs from the corresponding ones in~\cite{AntezzaPRL2009,Guimond2019} by a factor of two. This increased noise originates from the scattering of the Raman-laser which adds an additional contribution $\sim k_e^2 \sigma^2$ to the atomic center-of-mass diffusion (see~\appref{app:Fast_Atomic_Fluctuations}). 

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 \eqnref{eq:V_fluct}, averaged over the atomic fluctuations, reduces to the pinned-atom expression \eqnref{eq:Vop}, with a renormalized driving frequency $\tilde{\W}_0=\W_0 \exp(-k_e^2\sigma^2)\simeq (1-\sigma^2k_e^2)\W_0$.
Proceeding as in \secref{sec:PinnedAtoms}, we simulate the Sch\"odinger 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 calculate $\tilde{G}_{\bold{i}\jj}$ by averaging \eqnref{eq:G_ij_op} over one hundred realizations of the atomic positions assumed to be distributed around the lattice sites according to $P(\rr)$.
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 $\W_0$.
We present the results of the numerical simulations in~\figref{Fig:Fig4}.c.
As a consequence of the saturation of the dark-mode decay rate (\figref{Fig:Fig4}.b), the state-preparation fidelity decreases at large $N$ for a fixed lattice spacing $a/\lambda_e$ as evidenced in~\figref{Fig:Fig4}.c, where we consider the case $\sqrt{2}\sigma/a=0.05$. For larger values of $\sigma$ 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 \figref{Fig:Fig4}.c). 
The results in \figref{Fig:Fig4}.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\gtrsim40$ for the case $\sqrt{2}\sigma=0.01a$ (\figref{Fig:Fig4}.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 \figref{Fig:Fig5}.
In \figref{Fig:Fig5}.a, we compare,  for different values of $\sigma$, the minimal error for preparing the state $\ket{10}_\text{L}$ of two parallel arrays of $N=20$ atoms separated by a distance $l=a$ as function of the lattice spacing $a/\lambda_e$.
\begin{figure}
	\includegraphics[width=\columnwidth]{./Fig5}
	\caption{{\bf Gate Error for parallel arrays of fluctuating atoms.} a) Error for the preparation of the state $\ket{10}_\text{L}$ as function of $a/\lambda_e$, for different values of $\sigma$. Other parameters: $N=20$, $l=a$, $\delta=100\Gamma_0$. b) $\sqrt{\text{iSWAP}}$-gate error as function of $N$ for $a/\lambda_e=0.10$ and $2\sigma=0.1a$. Different colored markers corresponds to different values of $l$ (see legend). c) Left panel: $\sqrt{\text{iSWAP}}$-gate error as function of $N$ for $l=2a$ and $a/\lambda_e=0.10$. Right panel: $\sqrt{\text{iSWAP}}$-gate error as function of $a/\lambda$ for $l=2a$ and $N=12$. In both panels different values of $\sigma$ corresponds to different coloured markers (see legend). In all panel we calculated the coupling in \eqnref{eq:Ham_fam} by averaging over 100 realisation according to $P(\rr)$.}\label{Fig:Fig5}
\end{figure}
As expected, the error grows with $\sigma$ as a consequence of the contribution $(\sigma k_e)^2 \Gamma_0$ to the dark-mode decay rate. This effect is particularly prominent at large values of $a/\lambda_e$ where the main contribution to the non-linearity comes from the Zeno effect. For small $a/\lambda_e\lesssim 0.05$, values corresponding to different $\sigma$ yield similar results because the main contribution to the error comes from populating the state $\ket{01}_\text{L}$, as for the case of pinned atoms (see discussion in \secref{ssec:PinnedAtoms_Parallel}).
In \figref{Fig:Fig5}.b, we compare the $\sqrt{\text{iSWAP}}$-gate error, calculated according to \eqnref{eq:F_swap}, as a function of $N$ for different value of the separation $l$ between the arrays, assuming a position noise characterized by $\sqrt{2}\sigma=0.10a$. 
The gate fidelity deteriorates for increasing $l$, much faster than for the case of pinned atoms [Cf. \figref{fig:Fig3}.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 \figref{Fig:Fig5}.c, we study in more detail the particular case of $l=2a$, which yields better fidelities than other values of $l$ in~\figref{Fig:Fig5}.b~\footnote{We consider only integer values of $l/a$. Better results can be obtained if this assumption is relaxed.}.
In particular, we compare the scaling of the error with both $N$ (left panel) and $a/\lambda_e$ (right panel) for different values of position noise $\sigma$. 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 \figref{Fig:Fig5}.c). 
The left panel of \figref{Fig:Fig5}.c shows that at larger lattice separation positional noise always deteriorate the gate fidelity, as expected for a reduced Zeno blockade. However, for $a\lesssim 0.10\lambda_e$, a value of $\sqrt{2}\sigma=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 $\ket{11}_\text{L}$ and $\ket{S_2}$, which leads to a larger population in $\ket{11}_\text{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 \eqnref{eq:Ham_fam} do not yield correct estimates for $a/\lambda_e\lesssim 0.1$ (hatched region in \figref{Fig:Fig5}.c) as we argue in the following. 

\subsection{Limitation of the current description of the atomic fluctuations and alternative approaches}\label{sec:Limitation_AtomicMotion}

The treatment of the atomic fluctuations presented in the previous section does not take into account corrections due to a finite velocity of the atoms~\cite{Cooper1974,Power1974,Berman1997,Zhu2016,Bettles2020}. 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~\footnote{The perturbative approached followed here was presented in~\cite{Guimond2019}. In the present case the form of the dipole coupling rate is difference due to the Raman transition.}.  
As in~\secref{sec:FastMotion}, we consider a subwavelength atomic array ($a<\lambda_e/2$) where atoms are well localized around their trapping position, namely $r_0\ll a$, with $r_0\equiv \sqrt{\hbar/2m \w_\text{T}}$ the center-of-mass zero point motion. Under this assumption and for low thermal phonon number $n_\text{th} < 1$, the system is in the Lamb-Dicke regime, $\eta \equiv \sigma k_e\ll1$, for $\sigma\equiv r_0\sqrt{2n_\text{th}+1}$ and we approximate the coupling between internal and external degrees of freedom by a power expansion of $G(\rr)$ to second order in $\eta$~\cite{Guimond2019},
\be\label{eq:Ham_Motion_Expansion}
	\Hop\simeq \Hop_0 + \Hop_\text{I1}+\Hop_\text{I2}.
\ee
Here, $\Hop_0$ describes the dynamics in the absence of mechanical effects of light and is simply given by \eqnref{eq:Ham_2Arrays_Pinned} with the addition of the mechanical energy of the atoms. 
The term $\Hop_\text{I1}$ ($\Hop_\text{I2}$) represents the first (second) order correction in the atomic center-of-mass displacement. They have the general form
\bea
	\Hop_\text{I1} &\equiv& \varepsilon^2 \sum_{\bold{i},\jj} \pare{\hat{\rr}_\bold{i}-\hat{\rr}_\jj}\cdot \grad G(\RR_\jj,\RR_\bold{i}) \hat{\sigma}_{21}^\jj\hat{\sigma}_{12}^\bold{i} ,\label{eq:Ham_I1}\\
	\Hop_\text{I2} &\equiv& \frac{\varepsilon^2}{2} \sum_{\jj\bold{i}} \spare{\pare{\hat{\rr}_\bold{i}-\hat{\rr}_\jj}\cdot \grad}^2G(\RR_\bold{i},\RR_\jj)  \hat{\sigma}_{21}^\jj\hat{\sigma}_{12}^\bold{i}.\label{eq:Ham_I2}
\eea
Let us here remark that \eqnref{eq:Ham_Motion_Expansion} reduces to \eqnref{eq:Ham_fam_approx}, as expected, if one traces over the center-of-mass degrees of freedom assuming the atoms in a thermal state of their mechanical motion.

The driving Hamiltonian \eqnref{eq:V_fluct} up to second order in the center-of-mass fluctuations reads
\be\label{eq:V_Motion_Expansion}
	\Vop \simeq \Vop_0 + \Vop_1 + \Vop_2.
\ee
Here, $\Vop_0$ is the driving Hamiltonian for arrays of pinned atoms as given in \eqnref{eq:Vop_2array}. The term $\Vop_1$ ($\Vop_2$) represents the first (second) order correction in the Lamb-Dicke parameter. They read (see~\appref{app:LambDicke})
\be\label{eq:V1_main}
\begin{split}
\Vop_1 \equiv& -\im \hbar\W_0 \sum_{j=1}^N \cos(K_z aj_z)\Big[\eta_z\pare{\bdop_{zj}+\bop_{zj}}\\
& +\eta_x \pare{\bdop_{xj}+\bop_{xj}}\Big]\sy_j,
\end{split}
\ee
and
\be\label{eq:V2_main}
\begin{split}
\Vop_2 \equiv& \im \hbar\W_0 \sum_{j=1}^N \sin(K_z aj_z)\Big[\eta_z^2\pare{\bdop_{zj}+\bop_{zj}} ^2\\
&+2\eta_z\eta_x\pare{\bdop_{zj}+\bop_{zj}}\pare{\bdop_{xj}+\bop_{xj}}\\
&+\eta_x^2 \pare{\bdop_{xj}+\bop_{xj}}^2 \Big]\sy_j.
\end{split}
\ee
Here, we assumed the driving laser to be directed orthogonal to the $y$-axis. Moreover, we defined $\eta_{x,z}\equiv K_{x,z}r_0$, where $\xop_j \equiv r_0(\bdop_{xj}+\bop_{xj})$, $\zop_j \equiv r_0(\bdop_{zj}+\bop_{zj})$, and $K_x \equiv k_a \sin\alpha -k_b \sin\beta$, while $K_z$ is given in \eqnref{eq:q_effective}.

We now proceed to study the dynamics generated by \eqnref{eq:Ham_Motion_Expansion} and \eqnref{eq:V_Motion_Expansion} in the fast-motion regime $\w_T\gg\Gamma_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_\text{th}=0$).
In particular, we compare the results obtained from exact diagonalisation of the total Hamiltonian \eqnref{eq:Ham_Motion_Expansion} and \eqnref{eq:V_Motion_Expansion}, and the results obtained from an effective model in which the atomic motion has been adiabatically eliminated keeping contribution up to second order in $\eta$ (\appref{app:LambDicke}). 
This effective model contains both contributions from processes which do not change the motional state of the atoms (first order in $\Hop_\text{I2}$ and $\Vop_2$) and from second order processes via intermediate excited states of the atomic motion (second order in $\Hop_\text{I1}$ and $\Vop_1$).
In~\figref{fig:Fig6}.a, we show the error for preparing the state $\ket{10}_\text{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 $\w_\text{T}/\Gamma_0=100$. We assume one array to be detuned from the target one by $\delta$.
\begin{figure}
	\includegraphics[width=\columnwidth]{./Fig6}
	\caption{{\bf Gate error for arrays of moving atoms: perturbative treatment.} a) Error for preparing the dark state $\ket{10}_\text{L}$ for two parallel arrays separated by $l=a$. Solid markers refer to exact digonalisation (ED) results while empty markers refer to the effective model (Eff.) where the motion is adiabatically eliminated in its ground state to second order in $\eta$. b) Total average center-of-mass phonon population for $\delta/\Gamma_0=500$ (circles) and $\delta/\Gamma_0=100$ (squares). For both panels a) and b) we used $\eta=2\pi\times 0.05 a/\lambda_e$.
	c)  $\sqrt{\text{iSWAP}}$-gate error for parallel arrays of fluctuating atoms separated by $l=2a$ as function of $a/\lambda_e$. Different values of atomic center of mass fluctuations $r_0/a$ are compared (see legend).
For all the plot the following values of additional parameters have been used : $n_\text{th}=0$, $N=2$, phonon Hilbert space dimension $d=2$, and $\w_\text{T}/\Gamma_0=100$.}\label{fig:Fig6}
\end{figure}
The results of exact diagonalization and of the effective model show excellent agreement at large lattice spacings $a/\lambda_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 (\figref{fig:Fig6}.b). 
In~\figref{fig:Fig6}.c, we show the dependence on $a/\lambda_e$ of the error in performing a $\sqrt{\text{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/\lambda_e$ at which the center-of-mass motion is excited depends on how tight the trap is: for smaller $r_0/a$ (tighter trap) fluctuations are excited for smaller $a/\lambda_e$.

\figref{fig:Fig6}.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 \eqnref{eq:Ham_I1}, \eqnref{eq:Ham_I2}, \eqnref{eq:V1_main} and \eqnref{eq:V2_main}.
Precisely, the excitation of the atomic fluctuations is predominantly due to the atomic recoil involved in the dipole-dipole interaction (\eqnref{eq:Ham_I1} and \eqnref{eq:Ham_I2}) 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_{\bold{i}\jj}$ for small separation. 
Because of the enhancement in the the coupling rate in \eqnref{eq:Ham_I1} and \eqnref{eq:Ham_I2} for small atomic separation, the perturbative approach to the atomic fluctuations developed here is valid only for sufficiently small $\eta$ such that $\eta k_e^{-1}|\pa{\alpha}G(\RR)/G(\rr)|\ll1$ and $\eta^2 k_e^{-2}\pa{\alpha}\pa{\beta} G(\RR)/G(\rr)\ll1$ for $\alpha,\beta=x,y,z$.
Furthermore the critical values of $\eta$ for which the perturbative approach is justified depends on $N$: for fixed $\eta$ one observes the appearence of unphysical eigenvalues of \eqnref{eq:Ham_Motion_Expansion} 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\gtrsim 0.1\lambda_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~\secref{sec:FastMotion}.

\section{Discussion}\label{sec: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 $\ket{e}\leftrightarrow\ket{g_1}$, (ii) a particular atomic structure which comprises two distinct Raman transitions connecting $\ket{g_2}$ to $\ket{g_1}$ via the intermediate levels $\ket{e}$ and $\ket{f}$ where $\w_f \geq \w_e$, (iii) a negligible dephasing rate of $\ket{g_2}$ ($\kappa_0\ll\Gamma_0$), (iv) fast atomic center-of-mass fluctuations $\w_\text{T} \gg \Gamma_0$, and (v) the Lamb-Dicke regime for the atomic motion.
Let us now show how these requirements can be met using ultracold alkaline-earth atoms in an optical lattice. The relevant level structure for such atoms is shown in \figref{fig:Fig7} where we marked transitions and levels used to implement the scheme illustrated in~\figref{fig:Fig1}.
\begin{figure}
	\includegraphics[width=0.9\columnwidth]{./Fig7}
	\caption{{\bf Relevant level structure (not to scale) of alkaline-earth-metal atoms.} We highlight the transitions and levels used to implement our scheme with the same notation of \figref{fig:Fig1}.}\label{fig:Fig7}
\end{figure}
The ground-state manifold levels $\ket{g_{1,2}}$ are encoded in the levels $^3P_0$ and $^3P_2$ respectively, while the excited state $\ket{e}$ is encoded in $^3D_1$, and the auxiliary state $\ket{f}$ used for the excitation of the array's collective dark modes is encoded in $^3S_1$.
Decay from $^3P_0$ to $^1S_0$ is a forbidden transition which happens at a rate $\approx 10 \text{mHz}$ for Sr, while decay from $^3P_2$ to $^1S_0$ has an even longer predicted lifetime~\cite{Ludlow2015}. The two levels $\ket{g_1}$ and $\ket{g_2}$ can thus be considered stable.
Alkaline-earth atoms exhibit long range dipole interactions on the transition $^3P_0- ^3\!D_1$, which combined with the possibility of creating optical lattices using transitions from $^3P_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 $\ket{e}$ and $\ket{g_1}$ is $2.6~\mu\text{m}$ which allows to attain a subwavelength array with lattice spacing $a/\lambda_e \approx 0.08$ for an optical lattice with wavelength $\lambda_\text{opt}\simeq 400~\text{nm}$~\cite{Olmos2013}.
Additionally, the transition between $^3P_0- ^3\!\!D_1$ has a linewidth $\gamma_1/2\pi\approx 290 \text{kHz}$ which is broader than the linewidth $\gamma_2/2\pi\approx 10 \text{kHz}$ for $^3P_2- ^3\!D_1$~\cite{Zhou2010} resulting in a weak effective dephasing $\kappa_0/\Gamma_0 \approx 3\times 10^{-3}$ of the level $\ket{g_2}$ after the elimination of the level $\ket{e}$.
The Raman-driving of deep-subwavelenght arrays can be realised with $\w_f/\w_e \simeq 3.8$ using the rapid transitions between $^3S_1-  ^3\!P_2$ ($\gamma/2\pi\approx 45~\text{MHz}$) and between $^3S_1- ^3\!P_0$ ($\gamma/2\pi\approx 10~\text{MHz}$)~\cite{Zhou2010}.
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~\secref{sec:Limitation_AtomicMotion} 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 $\sim 10-10^3 ~\text{kHz}$ for standard cold-atoms experiments~\cite{Ludlow2015}. Note that the recoil energy of the transition $\ket{e}\leftrightarrow\ket{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 $\eta=r_0 k_e \ll1$ with respect to the relevant dipole transition.
For these values of the trap depth, the trap frequency is of order $\w_\text{T}\sim 100\text{kHz}$ which is comparable to the decay rate $\gamma_1$ of the transition $\ket{e}\leftrightarrow \ket{g_1}$. Hence, it should be feasible to enter the fast atomic motion regime ($\w_\text{T}/\Gamma_0 < 1$) as soon as $\W/\Delta\lesssim 0.1$.
Lattices of cold alkaline-earth atoms thus represent a promising platform for the realization of the scheme proposed in \secref{sec:PinnedAtoms} and \secref{sec:Motion}.

It is important to point out that the dephasing arising from the Raman driving (\figref{fig:Fig1}.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 equal the dephasing rate $\kappa_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 $\ket{e}\leftrightarrow\ket{g_1}$, and a two-photon driving for the transition $\ket{e}\leftrightarrow\ket{g_2}$ as proposed in~\cite{Porras2008}. 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~\secref{sec:Motion}).

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~\cite{Betzig2006,Hell2007}. 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. 

\section{Conclusions}\label{sec: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~\cite{Plankensteiner2015,He2020a,He2020b}.  
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 center-of-mass fluctuations place a bound on the size $N$ of the array as well as on the lattice spacing $a/\lambda_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.
Is 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~\cite{Sorensen1999,Sorensen2000}. 
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 dipole-dipole interactions~\cite{Shahmoon2019,Shahmoon2020}.
Another interesting direction, which we are currently pursuing, is the development of better theoretical models which would allow to extend the results of \secref{sec:Motion} to larger arrays while accurately taking into account possible excitation of the center-of-mass 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. 

\begin{acknowledgments}
  C.C.R.\ would like to thank Ana Asenjo-Garcia, Giuseppe Calaj\'o, Henrik Dreyer, Giacomo Giudice, Johannes Kn\"orzer, Daniel Malz, Nicola Pancotti, and Daniel Robaina for enlightening discussions. The moral support of Claudio Benzoni during these particularly difficult times is also gratefully acknowledged.
  C.C.R is particularly thankful to Daniel H\"ummer for providing insightful comments on the manuscript.
  C.C.R.\ and J.I.C.\ acknowledge funding from ERC Advanced Grant QUENOCOBA under the EU Horizon 2020 program (Grant Agreement No. 742102). T.S. was supported by the NSFC (Grants No. 11974363).
\end{acknowledgments}

