Generation of Photonic Matrix Product States with Rydberg Atomic Arrays

Zhi-Yuan Wei^1,2, Daniel Malz^1,2, Alejandro González-Tudela³, J. Ignacio Cirac^1,2 ¹ Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, D-85748 Garching, Germany ²Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 München, Germany ³ Instituto de Física Fundamental IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain

(23rd June 2025)

Abstract

We show how one can deterministically generate photonic matrix product states with high bond and physical dimensions with an atomic array if one has access to a Rydberg-blockade mechanism. We develop both a quantum gate and an optimal control approach to universally control the system and analyze the photon retrieval efficiency of atomic arrays. Comprehensive modeling of the system shows that our scheme is capable of generating a large number of entangled photons. We further develop a multi-port photon emission approach that can efficiently distribute entangled photons into free space in several directions, which can become a useful tool in future quantum networks.

I Introduction

The generation of multiphoton quantum states is at the core of many quantum technologies, including computing [1], cryptography [2], networks [3], or sensing [4]. The standard method utilizes parametric down-conversion (PDC) [5], where the desired state is generated by heralding its presence by post-selection. While very successful in several scenarios, this method possesses certain limitations, notably the exponential dependence of the success probability on the number of photons. These limitations can be overcome, for instance, by using emitters coupled to photonic waveguides, where one first prepares an entangled state of the emitters, which is then mapped into the multiphoton state using collective decay [6, 7]. This requires an effective non-linear interaction (coherent or dissipative) among the emitters, and the number of photons is limited by the number of emitters. An alternative way of generating multiphoton states at a higher rate is to do it sequentially using a photon source with long coherence time [8, 9, 10, 11, 12, 13]. The class of states that can be generated in this way coincides with the set of matrix product states (MPS) [14], a class of tensor network states that have been extensively studied in condensed matter physics [15, 16, 17, 18]. Some of the sequential generation protocols have recently been experimentally realized in quantum dots [19] and circuit QED [20, 21]. Going beyond one-dimensional MPS, it has been proposed to produce a subclass of projected entangled pair states (PEPS) [22, 23, 18] by allowing the emitted photons to travel back and interact with the photon source again [24, 25, 26], coupling multiple MPS sources [27, 28, 29, 30], or applying linear optical operations on certain resource states [31].

Most of the above-mentioned schemes produce tensor network states of photonic qubits with a limited bond dimension (typically, $D=2$ ), which contains important states like the GHZ state [32] and the cluster state [33]. A notable exception is the scheme studied in Ref. [25], which is capable of probabilistically producing tensor network states with small $D>2$ utilizing the PDC process in an optical loop. Nevertheless, It is very desirable to deterministically produce tensor network states with moderately higher bond dimensions, as the ground states of a large variety of many-body systems are encoded there [17, 34, 35, 36, 37]. Tensor network states of high bond dimensions can also be instrumental for metrological purposes [38, 39]. Moreover, tensor network states with higher physical dimensions (photonic qudit) further allows one to explore high-dimensional entanglement [40] with many applications, for instance, to produce bosonic quantum error-correcting codes [41], or act as high-dimensional resource states for state preparation [42] and computation [43, 44, 45, 46, 47].

A very attractive way of generating photonic states is by using atomic ensembles and arrays in free space [48, 49, 50, 51, 52, 53, 54, 55]. The most significant advantage is that that system does not require the presence of a cavity or a waveguide to collect the photons; interference effects are harnessed instead to emit the photons in a determined direction with a high probability, a phenomenon that is enhanced by collective effects (superradiance) in three spatial dimensions. In Ref. [49] a scheme was proposed to efficiently generate photonic states in atomic arrays. The basic idea is to combine those interference/collective phenomena with a Rydberg blockade mechanism [56, 57], so that in a first step an entangled atomic state is prepared, which is mapped in a second step into photonic degrees of freedom. This results in one or two outgoing photonic wavepackets in a state that is predetermined by the procedure.

In this paper we show how, in a similar setup, it is possible to generate quantum photonic states in a sequential way, where the output states consist of many wavepackets outgoing in predetermined directions and whose quantum state is an arbitrary MPS. This generalizes the class of states that can be produced on a similar setup [12] (a subclass of MPS with $D=2$ ) with a simpler atomic level scheme. We illustrate the method with simple examples, like the cluster state [33] (an MPS with $D=2$ ) and the generalized GHZ states [42] (MPS with high bond dimension and physical dimension). We analyze the performance of the source by careful simulation of the impact of imperfections in this scheme and characterizing the photon retrieval efficiencies of atomic arrays with various geometries. We predict that with this scheme and realistic conditions it should be possible to generate several tens of entangled photons, and extend the setup to send the photons to different directions by controlling the phase-matching direction of the collective excitations. Thus, this setup provides us with an excellent platform for the sequential generation of photonic MPS, which can efficiently distribute the entanglement among multiple ports.

The rest of the manuscript is structured as follows. In Section II, we present our scheme to generate arbitrary MPS using a Rydberg-blockaded atomic array. In Section III, we develop both a quantum gate approach and a quantum optimal control approach [58, 59] to implement arbitrary unitaries in our system, illustrate their usage by constructing unitaries required for generating the cluster state and the generalized GHZ state, and analyze the impact of imperfections during the unitary evolution process. In Section IV, we characterize the photon retrieval efficiencies of atomic arrays with various geometries and the impact of imperfections during the photon emission process. Then we compute the number of entangled photons that can be created with our scheme in Section V and extend it to operate as a multi-port device in Section VI. We propose a possible experimental realization in Section VII and summarize our work in Section VIII.

II MPS generation with a Rydberg array

The general protocol for the sequential generation of MPS is presented in Ref. [14]. In this section, we first briefly review this protocol in Section II.1, and then introduce our Rydberg-blockaded atomic array setup in Section II.2.

II.1 MPS generation protocol

Let us consider a photon source with Hilbert space ${\cal H}_{{\rm{src}}}^{\left({D,d}\right)}$ consisting of a $D$ -dimensional ancilla with Hilbert space ${\cal H}_{D}$ coupled to a $d$ -dimensional emitter with Hilbert space ${\cal H}_{d}$ , that ${{\cal H}^{\left({D,d}\right)}_{{\rm{src}}}}={{\cal H}_{D}}\otimes{{\cal H}_{d}}$ . We assume that we can implement any arbitrary unitary operation ${U_{\rm src}}\in\rm{SU}({D\cdot d})$ in ${{\cal H}_{{\rm{src}}}^{\left({D,d}\right)}}$ , and we can trigger a photon emission process ${M_{P}}\mathrel{\mathop{\mathchar 58\relax}}{{\cal H}_{d}}\to{{\cal H}_{d}}\otimes{{\cal H}_{{\rm{\rm ph}}}}$ that generates a $d$ -dimensional photonic qudit with Hilbert space ${{\cal H}_{\mathrm{\rm ph}}}$ :

{M_{P}}\mathrel{\mathop{\mathchar 58\relax}}{\left|i\right\rangle_{d}}\to{\left|0\right\rangle_{d}}{\left|i\right\rangle_{{\rm{\rm ph}}}},\qquad i\in\left({0,...,d-1}\right).

(1)

The sequential photon generation protocol starts with an initial state ${\left|{{\varphi_{I}}}\right\rangle_{D}}{\left|0\right\rangle_{d}}$ without excitations on the emitter. In each round, we first apply a unitary on the photon source ${U_{\left[i\right]}}\in{{\cal H}_{{\rm{src}}}^{\left({D,d}\right)}}$ and then trigger a photon emission $M_{P}$ , which produces a photonic qudit. Thus, in this protocol, the unitaries always act on states of form ${\left|\varphi\right\rangle_{D}}{\left|0\right\rangle_{d}}$ as

{U_{\left[i\right]}}{\left|\varphi\right\rangle_{D}}{\left|0\right\rangle_{d}}=\sum\limits_{{j}=0}^{d-1}{V_{\left[i\right]}^{{j}}}{\left|\varphi\right\rangle_{D}}{\left|j\right\rangle_{d}}.

(2)

After operating the protocol for $n$ rounds (shown in Fig. 1(b)), we get the final state

\left|\Psi\right\rangle={M_{P}}{U_{\left[n\right]}}...{M_{P}}{U_{\left[1\right]}}\left|{{\varphi_{I}}}\right\rangle.

(3)

Assuming that after the $n$ -th photon emission, the photon source disentangles from the photonic state, such that $\left|\Psi\right\rangle={\left|{{\varphi_{F}}}\right\rangle_{D}}{\left|0\right\rangle_{d}}\otimes\left|{{\psi_{{\rm{MPS}}}}}\right\rangle$ , we get a $n$ -qudit photonic MPS of the form [14] (see Appendix A for details)

\left|{{\psi_{{\rm{MPS}}}}}\right\rangle_{\rm ph}=\sum\limits_{{i_{1}},...,{i_{n}}=0}^{d-1}{{}_{D}\left\langle{{\varphi_{F}}}\right|V_{\left[n\right]}^{{i_{n}}}...V_{\left[1\right]}^{{i_{1}}}{{\left|{{\varphi_{I}}}\right\rangle}_{D}}{{\left|{{i_{n}}...{i_{1}}}\right\rangle}_{{\rm{\rm ph}}}}}.

(4)

If one is able to generate arbitrary unitaries in ${{\cal H}^{\left({D,d}\right)}_{{\rm{src}}}}$ , then the $D$ -dimensional matrices $\{{V_{\left[k\right]}^{{i_{k}}}}\}$ can be arbitrary, as long as the isometry condition $\sum\nolimits_{{i_{k}}=0}^{1}{V_{\left[k\right]}^{{i_{k}}{\dagger}}V_{\left[k\right]}^{{i_{k}}}}={I_{D}}$ is fulfilled [14]. This shows that in principle, the whole family of MPS can be generated in this way. Furthermore, by including another $m$ -level ancilla, it is possible to create superposition of the $m$ -level ancilla and photonic MPSs of the form (see Appendix A for details)

\left|{{\psi_{{m,\rm{ph}}}}}\right\rangle=\sum\limits_{i=0}^{m-1}{{\alpha_{i}}{{\left|i\right\rangle}_{{m}}}{{\left|{\psi_{{\rm{MPS}}}^{i}}\right\rangle}_{{\rm{\rm ph}}}}}.

(5)

This additional ancilla can be either a physical system or a part of the Hilbert space of a high-dimensional ancilla. Such superposed states find many applications in quantum networks [60], and realize the functionality of quantum random access memory [61].

A good physical platform to generate MPS of high bond and physical dimensions with the above protocol requires one to

1.

have a Hilbert space or its subspace of structure ${\cal H}_{{\rm{src}}}^{\left({D,d}\right)}$ with high dimensions,
2.

be able to efficiently implement unitary operations [cf. Eq. 2],
3.

be able to realize the photon emission $M_{P}$ [cf. Eq. 1].

In the following subsections we show that a Rydberg blockaded [56, 57] atomic array is well-suited for this task.

II.2 System setup

We consider an array of $N$ atoms (sketched in Fig. 1(a)), with three hyperfine ground states $\left|{g}\right\rangle$ , $\left|{l}\right\rangle$ and $\left|{q}\right\rangle$ , a Rydberg state $\left|{r}\right\rangle$ , and an auxiliary excited state $\left|e\right\rangle$ used for the photon emission. The transitions between all these states are controlled with lasers (denoted with coloured arrows in Fig. 1(a)). We assume our system is inside the Rydberg blockade radius, such that at most one Rydberg excitation can be present in the system. In addition, the system is operated in the low-excitation regime, which means that most atoms are in the ground state $\left|g\right\rangle$ . We denote the many-body ground state of the system by $\left|0\right\rangle={\left|g\right\rangle^{\otimes N}}$ .

Refer to caption — Figure 1: Generation of photonic MPS with a Rydberg-blockaded array. (a) Schematic plot and single-atom level scheme of our protocol. We consider an atomic array within the Rydberg blockade radius. For each atom, there are three hyperfine ground states $\left|g\right\rangle,\left|l\right\rangle,\left|q\right\rangle$ coupled to an intermediate Rydberg state $\left|{{r}}\right\rangle$ via lasers denoted by arrows of different colours. A fast decaying excited state $\left|e\right\rangle$ is used to convert the collective excitation(s) on $\left|l\right\rangle$ to photon(s). (b) The sequential photon generation protocol. Starting from the atomic initial state $\left|{{\varphi_{I}}}\right\rangle$ , the photonic state is generated by successively applying the unitary operation ${U_{\left[i\right]}}$ followed by the photon emission process $M_{P}$ for $n$ times, which yields an $n$ -photon matrix product state. (c) The structure of collective Hilbert space ${\mathcal{H}_{\rm col}}$ in the low-excitation regime (without $\left|e\right\rangle$ ), with the basis defined in Eq. 12. The arrows indicate the same laser transitions with the same coloring as in (a). The subspace used in our scheme for generating the cluster state is shown in pink (see Section III).

Let us first consider the unitary evolution process in our protocol under ideal conditions, which means (i) the Rydberg blockade condition is perfect, and (ii) there are no decoherence processes. Since the state $\left|e\right\rangle$ is not involved in the unitary operations, we define $\sigma_{\alpha\beta}^{i}={\left|\alpha\right\rangle_{i}}\left\langle\beta\right|$ and write the Hamiltonian for our atomic array as follows

H\left(t\right)={H_{0}}+{H_{rg}}\left(t\right)+{H_{rl}}\left(t\right)+{H_{rq}}\left(t\right),

(6)

where ${H_{0}}$ contains the energy of the atomic levels (with $\hbar=1$ )

{H_{0}}={\omega_{r}}\sum\limits_{i=1}^{N}{\sigma_{rr}^{i}}+{\omega_{q}}\sum\limits_{i=1}^{N}{\sigma_{qq}^{i}}+{\omega_{l}}\sum\limits_{i=1}^{N}{\sigma_{ll}^{i}},

(7)

and we have three driving lasers, denoted by $\mathrm{L_{g}}$ , $\rm{L_{q}}$ and $\mathrm{L_{l}}$ , which are used to connect the levels $\left|g\right\rangle$ , $\left|q\right\rangle$ and $\left|l\right\rangle$ with the Rydberg level $\left|r\right\rangle$ during the unitary process. We define the collective excitation operator

S_{\alpha,{\bf{k}}}^{\dagger}=\sum\limits_{i=1}^{N}{{u_{i}}{e^{i{\bf{k}}\cdot{{\bf{r}}_{i}}}}\sigma_{\alpha g}^{i}},\qquad\alpha=r,q,l,

(8)

where the ${u_{i}}\sim O({1/\sqrt{N}})$ is the normalized laser profile for $\mathrm{L_{g}}$ with $\sum\nolimits_{i}{{{\left|{{u_{i}}}\right|}^{2}}}=1$ . The Hamiltonian ${H_{rg}}\left(t\right)$ corresponding to the laser $\mathrm{L_{g}}$ is

{H_{rg}}\left(t\right)=\frac{1}{2}\left({{\Omega_{rg}}\left(t\right){e^{-i{\omega_{rg}}t}}S_{r,{{\bf{k}}_{rg}}}^{\dagger}+{\rm{H}}.{\rm{c}}.}\right).

(9)

Starting from $\left|0\right\rangle$ , ${H_{rg}}(t)$ leads to a Rabi evolution between $\left|0\right\rangle$ and $S_{r,{{\bf{k}}_{rg}}}^{\dagger}\left|0\right\rangle$ with collective Rabi frequency ${\Omega_{rg}}\left(t\right)$ . The lasers $\mathrm{L_{l}}$ and $\mathrm{L_{q}}$ have plane-wave profiles, with Hamiltonians

\begin{array}[]{l}{H_{r\alpha}}\left(t\right)=\frac{1}{2}\sum\limits_{i=1}^{N}{\left({{\Omega_{r\alpha}}\left(t\right)\sigma_{r\alpha}^{i}{e^{i\left({{{\bf{k}}_{r\alpha}}\cdot{{\bf{r}}_{i}}-{\omega_{r\alpha}}t}\right)}}+{\rm{H}}.{\rm{c}}.}\right)},\\ \qquad\alpha=l,q.\end{array}

(10)

In the low-excitation regime, we can apply the Holstein-Primakoff approximation [62] to our system to treat the collective excitation operators $S_{l,{\bf{k}}}^{\dagger}$ and $S_{q,{\bf{k}}}^{\dagger}$ as bosonic operators, whereas the Rydberg collective excitation operator $S_{r,{\bf{k}}}^{\dagger}$ can be treated as a spin-1/2 creation operator due to the Rydberg blockade. We have

S_{r,{{\bf{k}}_{rg}}}^{\dagger}\approx{\sigma_{r}^{\dagger}},\quad S_{q,{{\bf{k}}_{rg}}-{{\bf{k}}_{rq}}}^{\dagger}\approx a_{q}^{\dagger},\quad S_{l,{{\bf{k}}_{rg}}-{{\bf{k}}_{rl}}}^{\dagger}\approx a_{l}^{\dagger}.

(11)

In our MPS generation protocol, we fix the laser momenta of $\mathrm{L_{g},L_{q}}$ and $\mathrm{L_{l}}$ , and hence omit the momentum dependence of the approximated operators (right hand side of Eq. 11) for notational simplicity. The lasers $\mathrm{L_{g},L_{q}}$ and $\mathrm{L_{l}}$ couple states in a collective Hilbert space ${{\cal H}_{{\rm{col}}}}$ consisting of a qubit formed by presence or absence of a collective Rydberg excitation with Hilbert space ${\cal H}_{r}$ , and two oscillators formed by collective excitations on $\left|q\right\rangle$ and $\left|l\right\rangle$ with Hilbert spaces ${\cal H}_{q}$ and ${\cal H}_{l}$ , such that ${{\cal H}_{{\rm{col}}}}={{\cal H}_{r}}\otimes{{\cal H}_{q}}\otimes{{\cal H}_{l}}$ . A basis of ${{\cal H}_{{\rm{col}}}}$ is spanned by

\left|{{s_{r}}{m_{q}}{n_{l}}}\right\rangle=\frac{1}{{\sqrt{{m_{q}}!{n_{l}}!}}}{({\sigma_{r}^{\dagger}})^{{s_{r}}}}{({a_{q}^{\dagger}})^{{m_{q}}}}{({a_{l}^{\dagger}})^{{n_{l}}}}\left|0\right\rangle.

(12)

Here, $s_{r}=0,1$ and $m_{q},n_{l}$ are positive integers. For notational simplicity, we further omit the label(s) of the mode(s) with no excitations, for instance we define $\left|{{s_{r}}{m_{q}}}\right\rangle\equiv\left|{{s_{r}}{m_{q}}{0_{l}}}\right\rangle$ . The structure of ${{\cal H}_{\rm col}}$ and the laser-induced transitions is depicted in Fig. 1(c).

The collective Hilbert space ${{\cal H}_{{\rm{col}}}}$ contains the the structure ${{\cal H}^{\left({D,d}\right)}_{{\rm{src}}}}$ needed for the MPS generation protocol in Section II.1. We identify the first $D$ Fock states in ${\mathcal{H}}_{q}$ as a $D$ -level ancilla, and the first $d$ Fock states in ${\mathcal{H}}_{l}$ as a $d$ -level quantum emitter. Thus this platform fulfills the requirement 1 for a good MPS source listed in Section II.1.

Using the definitions of the collective operators [Eq. 11], we can rewrite the system Hamiltonian [Eq. 6] as

H_{\rm col}\left(t\right)=H_{0}^{\rm col}+H_{rg}^{\rm col}\left(t\right)+H_{rl}^{\rm col}\left(t\right)+H_{rq}^{\rm col}\left(t\right),

(13)

where

{H_{0}^{\rm col}}={\omega_{r}}\sigma_{r}^{\dagger}{\sigma_{r}}+{\omega_{l}}a_{l}^{\dagger}{a_{l}}+{\omega_{q}}a_{q}^{\dagger}{a_{q}},

(14)

{H_{rg}^{\rm col}}\left(t\right)=\frac{{{\Omega_{rg}}\left(t\right)}}{2}\left({{e^{-i{\omega_{rg}}t}}\sigma_{r}^{\dagger}+{\rm{H}}.{\rm{c}}.}\right),

(15)

{H_{r\alpha}^{\rm col}}\left(t\right)=\frac{{{\Omega_{r\alpha}}\left(t\right)}}{2}\left({{e^{-i{\omega_{r\alpha}}t}}{a_{\beta}}\sigma_{r}^{\dagger}+{\rm{H}}.{\rm{c}}.}\right)\;\;\alpha=l,q.

(16)

By combining the universal controllability for the Jaynes-Cummings model [63, 64] and the Lemma 5.5 of Ref. [65], one can show that any arbitrary unitary in ${{\cal H}_{\rm col}}$ can be generated with time evolution under $H_{\rm col}\left(t\right)$ , and the optimal time cost to implement a unitary in $H_{\rm src}^{\left({D,d}\right)}$ scales as ${T_{D,d}}\sim{\left({Dd}\right)^{2}}$ (see Appendix H). Thus this platform fulfills the requirement 2 listed in Section II.1. We also estimate the scaling of implementing the unitary as a function of the dimension in Appendix H.

One can also engineer a on-demand photon emission process ${M_{P}}\mathrel{\mathop{\mathchar 58\relax}}{{\cal H}_{l}}\to{{\cal H}_{l}}\otimes{{\cal H}_{\rm ph}}$ by coupling $\left|l\right\rangle$ to an excited state $\left|e\right\rangle$ with a plane-wave profile laser ${\rm L}_{\rm e}$ , which triggers the emission of the photon(s) [49]

{M_{P}}\mathrel{\mathop{\mathchar 58\relax}}\left|{{s_{r}}{m_{q}}{n_{l}}}\right\rangle\to\left|{{s_{r}}{m_{q}}}\right\rangle{\left|{{n_{l}}}\right\rangle_{{\text{\rm ph}}}}

(17)

This fulfills the requirement 3 listed in Section II.1. Here, we assume that the atomic population only goes from $\left|e\right\rangle$ to $\left|g\right\rangle$ , which for instance can be realized by choosing $\left|e\right\rangle\leftrightarrow\left|g\right\rangle$ as a cycling transition [49]. Such high-efficiency photon retrieval of a spin-wave stored in an atomic array has been demonstrated experimentally [66, 67]. More details of the photon retrieval process are presented in Section IV.

With all three requirements for an efficient MPS source fulfilled, this setup can be used to generate photonic MPS with high bond dimension $D$ and physical dimension $d$ satisfying $D+d-2\ll N$ (required by the low-excitation regime) using the generic protocol introduced in Section II.1.

III Implementation of unitary operations

The Hamiltonian written in the collective basis [Eq. 13] is equivalent to a central qubit coupled to two oscillators. Universal control of qubit-oscillator systems has been studied in various setups, in particular the Jaynes-Cummings model of a single qubit coupled to an oscillator. While state preparation can be done using the Law-Eberly scheme [68], efficiently implementing arbitrary unitary operations is more challenging [69, 70, 64, 71, 72, 72]. In this section, we present efficient unitary control methods for this setup. We first provide a quantum gate approach that utilizes the AC Stark shift to implement a set of universal gates needed for generating photonic MPS of bond dimension $D=2$ and physical dimension $d=2$ in Section III.1. Then, we present a quantum optimal control scheme that achieves universal control in ${{\cal H}_{\rm col}}$ in Section III.2. We demonstrate these two approaches by constructing unitaries required for generating the cluster state [33] and the generalized GHZ state [42]. We further analyze the impact of imperfections during the photon generation protocol in Section III.3.

III.1 Quantum gate approach

The class of MPS of $D=2$ and $d=2$ contains many interesting states like the GHZ state [32] and the cluster state. To generate this class of states with our scheme, it is sufficient if both ancilla and emitter are two-level qubits. Thus, ${\cal H}^{\left({2,2}\right)}_{{\rm{src}}}={\rm{span}}\{{\left|0\right\rangle,\left|{{1_{l}}}\right\rangle,\left|{{1_{q}}}\right\rangle,\left|{{1_{q}}{1_{l}}}\right\rangle}\}$ (the basis of ${\cal H}^{\left({2,2}\right)}_{{\rm{src}}}$ is depicted in pink in Fig. 1(c), and cf. Eq. 12 for the definition of basis).

The procedure to implement a unitary in ${\cal H}^{\left({2,2}\right)}_{{\rm{src}}}$ consists of three steps.

1.

Apply a $\pi$ -pulse of the driving $\mathrm{L_{l}}$ to transfer the population from $\left|l\right\rangle$ to the Rydberg level $\left|r\right\rangle$ as $\left|{{1_{l}}}\right\rangle\to\left|{{1_{r}}}\right\rangle$ and $\left|{{1_{q}}{1_{l}}}\right\rangle\to\left|{{1_{r}}{1_{q}}}\right\rangle$ .

Implement the desired operation within a space

{\cal H}_{{\rm{mid}}}={\rm{span}}\left\{{\left|0\right\rangle,\left|{{1_{r}}}\right\rangle,\left|{{1_{q}}}\right\rangle,\left|{{1_{r}}{1_{q}}}\right\rangle}\right\}.

Here ${\cal H}_{{\rm{mid}}}$ is similar to ${\cal H}^{\left({2,2}\right)}_{{\rm{src}}}$ , but with the collective excitation(s) on $\left|l\right\rangle$ transformed to the Rydberg state $\left|r\right\rangle$ .

3.

Apply another $\pi$ -pulse with the laser $\mathrm{L_{l}}$ to transfer the population from $\left|r\right\rangle$ back to $\left|l\right\rangle$ as $\left|{{1_{r}}}\right\rangle\to\left|{{1_{l}}}\right\rangle$ and $\left|{{1_{r}}{1_{q}}}\right\rangle\to\left|{{1_{q}}{1_{l}}}\right\rangle$ .

The key step is to implement arbitrary unitaries in ${{\cal H}_{\rm mid}}$ , for which we need to construct a universal two-qubit gate set in ${{\cal H}_{\rm mid}}$ using the lasers ${\rm L}_{g}$ and ${\rm L}_{q}$ . $\mathrm{L_{g}}$ can realize arbitrary rotations in ${\cal H}_{r}$ that couples $\left|0\right\rangle\leftrightarrow\left|{{1_{r}}}\right\rangle$ and $\left|{{1_{q}}}\right\rangle\leftrightarrow\left|{{1_{r}}{1_{q}}}\right\rangle$ , whereas $\mathrm{L_{q}}$ couples $\left|{{1_{r}}}\right\rangle\leftrightarrow\left|{{1_{q}}}\right\rangle$ . However, ${\rm L}_{q}$ also leads to population leakage from $\left|{1_{r}1_{q}}\right\rangle$ to $\left|{{2_{q}}}\right\rangle$ , which is outside of ${{\cal H}_{\rm mid}}$ .

To avoid such leakage error, we propose to use the AC Stark shift [73, 74] induced by applying ${\rm L}_{q}$ with a large detuning $\Delta_{\rm ST}\gg{\Omega_{rq}}$ . This allows us to construct a phase gate of the form ${S_{T}}\left(\theta\right)=\textrm{diag}\left({1,{e^{-i\theta}},{e^{i\theta}},{e^{-2i\theta}},{e^{2i\theta}}}\right)$ with the basis $\{{\left|0\right\rangle,\left|{{1_{r}}}\right\rangle,\left|{{1_{q}}}\right\rangle,\left|{{1_{r}}{1_{q}}}\right\rangle,\left|{{2_{q}}}\right\rangle}\}$ . Using ${S_{T}}\left(\theta\right)$ together with single-qubit rotations, we can construct the SWAP gate and CNOT gate in ${{\cal H}_{\rm mid}}$ (see Appendix B for details). The SWAP gate, CNOT gate and the single-qubit rotations on ${{\cal H}_{r}}$ realize a universal gate set in ${{\cal H}_{\rm mid}}$ , and thus lead to the universal control of ${\cal H}^{\left({2,2}\right)}_{{\rm{src}}}$ using the quantum gate approach.

For example, to generate the cluster state of the form Eq. 4 with elements

\begin{gathered}{V^{0}}=\frac{1}{{\sqrt{2}}}\left({\begin{array}[]{*{20}{c}}1&0\\ 1&0\end{array}}\right),\quad{V^{1}}=\frac{1}{{\sqrt{2}}}\left({\begin{array}[]{*{20}{c}}0&1\\ 0&{-1}\end{array}}\right),\quad\hfill\\ \left|{{\varphi_{I}}}\right\rangle=\frac{1}{{\sqrt{2}}}\left({\left|0\right\rangle+\left|{{1_{q}}}\right\rangle}\right),\quad\left|{{\varphi_{F}}}\right\rangle=\left|0\right\rangle,\hfill\\ \end{gathered}

(18)

we need to implement two kinds of unitaries ${U^{\rm cl}_{\left[{i\neq n}\right]}}$ and ${U^{\rm cl}_{\left[{n}\right]}}$ (see Appendix B for details). Each application of ${U^{\rm cl}_{\left[{i\neq n}\right]}}$ followed by $M_{P}$ adds one site to the cluster state, and the last unitary ${U^{\rm cl}_{\left[{n}\right]}}$ followed by $M_{P}$ emits the last photon and disentangles the source from the photonic MPS. The quantum circuit for implementing ${U^{\rm cl}_{\left[{i\neq n}\right]}}$ is shown in Fig. 2(a). ${U^{\rm cl}_{\left[{n}\right]}}$ can be simply implemented with a two-photon Raman $\pi$ -pulse that swaps all population from $\left|q\right\rangle$ to $\left|l\right\rangle$ , which we discuss in detail in Section VII.

III.2 Quantum Optimal Control (QOC) approach

The quantum gate approach can be implemented in a relatively easy fashion since one only needs to sequentially implement the gates. However, the phase gate ${S_{T}}\left(\theta\right)$ is rather slow to implement, and it is less obvious how to extend this approach to universally control $\mathcal{H}_{\rm col}$ for the generation of photonic MPS with a larger bond and physical dimension $D$ and $d$ . QOC allows one to implement unitary operations in high-dimensional Hilbert spaces with a speed that potentially reaches the maximum speed allowed by quantum evolution [75, 76], and has already been applied to Rydberg atomic array experiments [77, 78]. Thus we develop a QOC approach to implement required unitaries in our protocol based on the Gradient Optimization of Analytic Controls (GOAT) algorithm [79].

In our MPS generation scheme we need to construct unitaries ${U_{\left[i\right]}}$ that behave as in Eq. 2. At each site $i$ of the MPS, the $D\times D$ matrices ${\{V_{\left[i\right]}^{j}\}_{j=0,...,d-1}}$ together form an isometry ${\hat{V}_{\left[i\right]}}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}_{D}\to{\cal H}_{{\rm{src}}}^{\left({D,d}\right)}$ that satisfies the isometry condition $\hat{V}_{\left[i\right]}^{\dagger}{\hat{V}_{\left[i\right]}}=I_{D}$ . ${\hat{V}_{\left[i\right]}}$ can be implemented through any unitary of the form

{U_{\left[i\right]}}=\left({\begin{array}[]{*{20}{c}}{{{\hat{V}}_{\left[i\right]}}}&{{B_{1}}}\\ O&{{B_{2}}}\end{array}}\right),

(19)

where the first row corresponds to ${\cal H}_{{\rm{src}}}^{\left({D,d}\right)}$ and the second to the rest of $\mathcal{H}_{\rm col}$ . Here, $O$ is a zero matrix, which physically means that ${U_{\left[i\right]}}$ does not cause the population to leak out of ${\cal H}_{{\rm{src}}}^{\left({D,d}\right)}$ . The parts $B_{1}$ and $B_{2}$ can be arbitrary as long as ${U_{\left[i\right]}}$ is unitary. This freedom of choice pose a challenge for certain gradient-based QOC algorithms like Krotov [80] and GRAPE [81], which requires backward evolution from ${U_{\left[i\right]}}$ . Thus we develop a QOC approach to efficiently synthesize ${U_{\left[i\right]}}$ by modifying the cost function of the GOAT algorithm such that it exactly captures the $\hat{V}_{\left[i\right]}^{j}$ and $O$ parts of $U_{\left[i\right]}$ (see the Appendix C for details). With this algorithm and the control Hamiltonian $H_{\rm col}\left(t\right)$ , we can obtain the pulse sequence for efficiently implementing desired ${U_{\left[i\right]}}$ . As a demonstration, the pulse sequence for $U_{\left[i\neq n\right]}^{\rm cl}$ is shown in Fig. 2(b). Here we only use the resonant laser couplings for $\mathrm{L_{g}},\mathrm{L_{q}}$ and $\mathrm{L_{l}}$ , which simplifies experimental realization.

Generation of high-dimensioanl MPS

We further demonstrate the ability of generating MPS of high bond and physical dimensions by considering the following generalized GHZ state [42]

|\mathrm{GHZ}(n,d)\rangle=\dfrac{1}{\sqrt{d}}\sum_{i=0}^{d-1}|i\rangle_{\rm ph}^{\otimes n}.

(20)

This state serves as a quantum resource state in the resource theory of genuinely multipartite entanglement under biseparability-preserving operations [42]. Defining $E\left({d,i}\right)$ as the $d$ -dimensional matrix with everywhere zero but only the $i$ -th diagonal element being one, $|\mathrm{GHZ}(n,d)\rangle$ can be written as an MPS of form Eq. 4 with both bond and physical dimension $D=d$ as

\begin{array}[]{*{20}{l}}{V_{\left[i\right]}^{j}={V^{j}}=E\left({d,j+1}\right),\quad j\in\left({0,...,d-1}\right),}\\ {\left|{{\varphi_{I}}}\right\rangle=\frac{1}{{\sqrt{d}}}\sum\limits_{i=0}^{d-1}{\left|{{i_{q}}}\right\rangle},\quad\left|{{\varphi_{F}}}\right\rangle=\left|0\right\rangle.}\end{array}

(21)

This state can be sequentially generated with the following two unitaries in Eq. 3

\begin{array}[]{l}U_{\left[{i\neq n}\right]}^{\rm gG}\mathrel{\mathop{\mathchar 58\relax}}\left|{{j_{q}}}\right\rangle\to\left|{{j_{q}}{j_{l}}}\right\rangle,\\ U_{\left[n\right]}^{\rm gG}\mathrel{\mathop{\mathchar 58\relax}}\left|{{j_{q}}}\right\rangle\to\left|{{j_{l}}}\right\rangle,\quad\forall j\in\left({0,...,d-1}\right).\end{array}

(22)

Similar as in the cluster state case, each application of $U_{\left[{i\neq n}\right]}^{\rm gG}$ followed by $M_{P}$ adds one site to the state, and $U_{\left[{n}\right]}^{\rm gG}$ followed by $M_{P}$ produces the last photon and disentangles the source from the photonic MPS. The QOC pulse sequence for generating $U_{\left[{i\neq n}\right]}^{\rm gG}$ in the case of $d=3$ is shown in Appendix C, and $U_{\left[{n}\right]}^{\rm gG}$ is implemented in the same way as $U_{\left[{n}\right]}^{\rm cl}$ by a two-photon Raman $\pi$ -pulse between $\left|q\right\rangle$ and $\left|l\right\rangle$ . The ability of generating $|\mathrm{GHZ}(n,d)\rangle$ clearly demonstrates that our scheme is particularly suitable to create high-dimensional entanglement.

III.3 Realistic modeling including imperfections

There will be various imperfections in real implementations of the above photon generation protocol. During unitary operations, the dominant imperfections are related to the Rydberg state $\left|r\right\rangle$ . For example, the Rydberg decay ( $\Gamma_{r}$ ) and dephasing ( $\Gamma_{\phi}$ ) lead to a typical lifetime of Rydberg states around $\sim 50\mu s$ [82], which is significantly shorter than the ground-state spin-wave coherence time that already can reach sub-second level [67]. Also, the finite Rydberg non-linearity $U$ provided by the van der Waals interaction [83] leads to a small probability to create double Rydberg excitations. The photon emission also has a finite retrieval efficiency $p_{\rm em}$ . Due to the above imperfections, operating our scheme will generate a $n$ -photon density matrix $\rho_{\rm ph}$ with a non-unit fidelity with respect to the targeted state ${{\cal F}_{{\rm{\rm ph}}}}=\left\langle{{\psi_{{\rm{MPS}}}}}\right|{\rho_{{\rm{\rm ph}}}}\left|{{\psi_{{\rm{MPS}}}}}\right\rangle$ . Due to the sequential nature of our scheme, ${{\cal F}_{\rm ph}}$ decays exponentially with the number of photonic qudits generated $n_{\rm ph}$ as ${{\cal F}_{\rm ph}}={e^{-\xi\cdot{n_{\rm ph}}}}$ . Here $\xi$ represents the error per photon.

Let us denote the maximal Rabi frequency during unitary implementations in our scheme as ${\Omega_{{\rm{max}}}}$ . We assume our scheme runs in the strong driving regime ${\Omega_{{\rm{max}}}}\gg{\Gamma_{r}}({\Gamma_{\phi}})$ and good blockade $U\gg{\Omega_{{\rm{max}}}}$ , which has already been reached by current experiments (for example see Refs. [82, 78]). Including Rydberg imperfections into the modeling, we find that the system density matrix $\rho\left(t\right)$ evolves as

\begin{array}[]{*{20}{l}}\begin{gathered}\dfrac{{d\rho\left(t\right)}}{{dt}}=-i\left[{H^{\prime}\left(t\right),\rho\left(t\right)}\right]\hfill\\ +\sum\limits_{i,\alpha=\left\{{g,l,q}\right\}}{\dfrac{{{\Gamma_{r}}}}{2}\left({2\sigma_{\alpha r}^{i}\rho\sigma_{r\alpha}^{i}-\sigma_{rr}^{i}\rho-\rho\sigma_{rr}^{i}}\right)}\hfill\\ \end{gathered}\\ {+\sum\limits_{i}{\dfrac{{{\Gamma_{\phi}}}}{2}\left({2\sigma_{rr}^{i}\rho\sigma_{rr}^{i}-\sigma_{rr}^{i}\rho-\rho\sigma_{rr}^{i}}\right)},}\end{array}

(23)

where the doubly excited Rydberg states are included in the Hamiltonian ${H^{\prime}\left(t\right)}$ . Since in our scheme we need to successively implement unitaries many times and possibly reach high excitations for high bond or physical dimension, the errors will accumulate over time and lead to involved dynamics. To track the long-time evolution of our atomic array with many atoms in a reliable way, we formulate an effective model in the collective basis that takes all these errors into account (see Appendix D for details). This effective model has a Hilbert space $\mathcal{H}_{\rm eff}$ that contains the Rydberg doubly excited states as well as the mixed states created by decoherence.

The finite photon retrieval efficiency $p_{\rm em}$ also reduces ${{\cal F}_{\rm ph}}$ . Its effect can be described by introducing the process map ${W_{P}}\mathrel{\mathop{\mathchar 58\relax}}{\mathcal{H}_{{\text{eff}}}}\to{\mathcal{H}_{{\text{eff}}}}\otimes{\mathcal{H}_{{\text{\rm ph}}}}$ , that maps the system density matrix to a joint density matrix of the system and a photonic qudit (see the construction of ${W_{P}}$ in Appendix E). Using ${W_{P}}$ and the solution of the master equation in $\mathcal{H}_{\rm eff}$ for the unitary operation process, we can compute the photonic state fidelity ${{\cal F}_{\rm ph}}$ with the help of the matrix product density operator (MPDO) [84] approach (see Appendix F for details).

We obtain the scaling of $\xi$ by numerical simulation of the protocol to generate cluster state using the QOC pulse sequence in Fig. 2(b). By obtaining ${{\cal F}_{\rm ph}}$ as a function of photon number $n_{\rm ph}$ , we can extract the error of ${{\cal F}_{\rm ph}}$ per photon $\xi=\xi\left({{\Gamma_{r}}/{\Omega_{\max}},{\Gamma_{\phi}}/{\Omega_{\max}},{\Omega_{\max}}/U,{p_{\rm em}}}\right)$ as shown in Fig. 3. We see that in the strong driving regime and with good blockade we have

\xi={\beta_{0}}+\frac{{{\beta_{r}}{\Gamma_{r}}+{\beta_{\phi}}{\Gamma_{\phi}}}}{{{\Omega_{\max}}}}+\frac{{{\beta_{U}}\cdot\Omega_{\max}^{2}}}{{{U^{2}}}}-{\beta_{{\text{em}}}}\cdot\log{p_{{\text{em}}}}.

(24)

The non-universal coefficients $\left\{{{\beta_{i}}}\right\}$ depend on the desired photonic state and the pulse shape. For example in Fig. 3, ${\beta_{0}}\approx 5\times{10^{-4}}$ represents the error due to the pulse imperfection, $\beta_{r}\approx 17.1$ , $\beta_{\phi}\approx 6.7$ and $\beta_{U}\approx 1.9$ correspond to the Rydberg decay, dephasing and double excitation error respectively, and $\beta_{\rm em}\approx 0.48$ corresponds to the photon retrieval error. Note that $\beta_{r}$ and $\beta_{\phi}$ can be further reduced by using a shorter pulse, at the expense of increasing $\beta_{U}$ . In the rest of this manuscript we assume our protocol operates in the regime of small error, in which Eq. 24 is valid. Note that if one implement unitaries via the quantum gate approach [cf. Section III.1], the Stark pulse will further leads to unwanted small residual population transfer [73, 85], which scales as $\Omega_{\rm ST}^{2}/2\Delta_{\rm ST}^{2}$ when ${\Delta_{\rm ST}}\gg{\Omega_{\rm ST}}$ . This error does not exist in Eq. 24 obtained from the QOC approach.

In actual experiments, the maximum Rabi frequency ${\Omega_{\max}}$ can be tuned on demand and is only bounded by the maximal laser power available, while ${\Gamma_{r}},{\Gamma_{\phi}},U$ are determined by the properties of the atomic array. Thus, by defining the total decoherence rate $\Gamma\equiv\beta_{r}{\Gamma_{r}}+\beta_{\phi}{\Gamma_{\phi}}$ and choosing ${\Omega_{\max}}={\left[{\Gamma{U^{2}}/\left({2{\beta_{U}}}\right)}\right]^{1/3}}$ , $\xi$ takes the optimized value

{\xi_{{\text{opt}}}}={\beta_{0}}+3\beta_{U}^{1/3}{\left({\dfrac{\Gamma}{{2\left|U\right|}}}\right)^{2/3}}-{\beta_{{\text{em}}}}\cdot\log{p_{{\text{em}}}}.

(25)

The scaling Eq. 24 and Eq. 25 applies to atomic arrays with arbitrary geometries as long as they lie within the Rydberg blockade radius. Different array geometries will affect the Rydberg non-linearity $U$ , the photon retrieval efficiency ${p_{\rm em}}$ , and the Holstein-Primakoff approximation. The analysis of array geometry dependence will be the focus of the next section.

IV Photon retrieval in atomic arrays

As we have shown in the previous section, we need to achieve high photon retrieval efficiency $p_{\rm em}$ because $\cal F_{\rm ph}$ decreases as $p_{{\text{em}}}^{\beta_{\rm em}\cdot{n_{{\text{\rm ph}}}}}$ . Defining the photon retrieval error as ${\varepsilon_{\rm em}}=1-{p_{\rm em}}$ , the maximal retrieval efficiency of an atomic ensemble generally scales inversely to its optical depth OD ${\varepsilon_{{\rm{em}}}}\sim{\rm{O}}{{\rm{D}}^{-1}}$ [86, 49]. In recent years, however, some theoretical works have suggested that the use of atomic arrays with optimized collective excitation profiles can lead to dramatic improvements in retrieval efficiency [87, 88, 89]. Thus, in this section, we will characterize how the photon retrieval efficiency scales for atomic arrays and the impact of imperfections during the photon emission process. See Appendix G for more details about this section.

The photon retrieval process is shown in Fig. 1(a), which consists of a Raman pulse of Rabi frequency ${\Omega_{el}}\left(t\right)$ and detuning $\Delta\left(t\right)$ that transfers the population on $\left|l\right\rangle$ to a fast-decaying excited state $\left|e\right\rangle$ , which only decays to the state $\left|g\right\rangle$ . During the emission process $\left|e\right\rangle\to\left|g\right\rangle$ , the emitted photon can be rescattered by atoms due to the dipole-dipole interactions between atoms [90]. Starting from a singly excited state

\left|{{1_{l}}}\right\rangle=S_{l,{{\bf{k}}_{rg}}-{{\bf{k}}_{rl}}}^{\dagger}\left|0\right\rangle=\sum\limits_{j}{{u_{j}}{e^{i\left({{{\bf{k}}_{rg}}-{{\bf{k}}_{rl}}}\right)\cdot{{\bf{r}}_{{j}}}}}\sigma_{lg}^{j}}{\left|g\right\rangle^{\otimes N}},

(26)

the decay dynamics is governed by the following non-Hermitian Hamiltonian ${H_{\rm em}}(t)={H_{el}}(t)+{H_{{\rm{DDI}}}}$ with

\begin{array}[]{l}{H_{el}}\left(t\right)=\frac{{{\Omega_{el}}\left(t\right)}}{2}\sum\limits_{i}{\left({\sigma_{el}^{i}{e^{i{{\bf{k}}_{el}}\cdot{{\bf{r}}_{i}}}}+{\rm H.c.}}\right)}-\Delta(t)\sum\limits_{i}{\sigma_{ee}^{i}},\\ {H_{{\rm{DDI}}}}=-{\mu_{0}}d_{eg}^{2}\omega_{eg}^{2}\sum\limits_{j,l}{{\bf{d}}_{j}^{*}}\cdot{{\bf{G}}_{0}}\left({{{\bf{r}}_{j}},{{\bf{r}}_{l}},{\omega_{eg}}}\right)\cdot{{\bf{d}}_{l}}\sigma_{j}^{eg}\sigma_{l}^{ge}.\end{array}

(27)

Here, ${d_{eg}}$ and ${{\bf{d}}_{j}}$ are the dipole matrix element and unit atomic polarization vector. The dyadic Green’s tensor is [91]

	$\displaystyle{{\mathbf{G}}_{0}}\left({{{\mathbf{r}}_{j}},{{\mathbf{r}}_{l}},{\omega_{eg}}}\right)=\frac{{{{\text{e}}^{{\text{i}}{k_{0}}R}}}}{{4\pi k_{0}^{2}{R^{3}}}}\left[\left({k_{0}^{2}{R^{2}}+{\text{i}}{k_{0}}R-1}\right){\mathbf{I}}\hfill\right.$		(28)
	$\displaystyle\left.+\left({-k_{0}^{2}{R^{2}}-3{\text{i}}{k_{0}}R+3}\right)\frac{{{\mathbf{R}}\otimes{\mathbf{R}}}}{{{R^{2}}}}\right],$		(28)

with ${\bf{R=r}}_{j}-{\bf{r}}_{l}$ and $k_{0}={\omega_{eg}}/c$ with $c$ being the speed of light. From Eq. 27 and Eq. 28 we can compute the evolution of the atomic array and the photon retrieval efficiency following the approach in Ref. [88]. The photon retrieval efficiency is defined as the probability of the photon to be emitted into some detection mode ${{\bf{E}}_{{\rm{det}}}}(\bf r)$ .

To achieve high retrieval efficiency, one needs to optimize the excitation profile of the collective excitation to best match the detection mode [87, 88, 89]. We study two types of excitation profiles. One is the optimal profile as studied in Ref. [88], which gives the minimal photon retrieval error $\varepsilon_{\rm em}^{\rm opt}$ . Another is the profile created by a Gaussian beam, which we call the Gaussian profile, with its photon retrieval error denoted as $\varepsilon_{\rm em}^{\rm Gauss}$ . To assess our photonic MPS generation protocol, we only use the Gaussian profile, as it is experimentally easier to implement.

For concreteness, we fix the distance between atoms to be $d_{0}=0.6\lambda_{eg}$ , but similar subwavelength interatomic distances give similar results. We consider two types of retrieval schemes. In the scheme shown in Fig. 4(a), the photon is emitted in two opposite directions, where the detection mode is a symmetric superposition of two vector Gaussian modes [92] along the photon emission directions. This retrieval setup has been shown to lead to favorable scaling of the retrieval efficiency in the case of a single layer two-dimensional array [88]. We assume that the photons collected from these two directions can be later recombined, for example using the setup in [30]. The geometry of our atomic array is characterized by its transverse size ${L_{v}}$ and longitudinal size ${L_{z}}$ .

The Fig. 4(b) shows the scaling of the photon retrieval error as a function of ${L_{v}}$ for the two-directional retrieval scheme. For the case of ${L_{z}}=1$ we reproduced the scaling $\varepsilon_{\rm em}^{\rm opt}\approx{\left({\log{L_{v}}}\right)^{2}}/L_{v}^{4}$ in Ref. [88]. When increasing $L_{z}$ , we obtain $\varepsilon_{\rm em}^{\rm opt}\approx c^{\rm opt}_{{L_{z}}}\cdot{\left({\log{L_{v}}}\right)^{2}}/L_{v}^{4}$ . For the Gaussian excitation profile, we also have $\varepsilon_{\rm em}^{\rm Gauss}\approx c^{\rm Gauss}_{{L_{z}}}\cdot{\left({\log{L_{v}}}\right)^{2}}/L_{v}^{4}$ when $L_{z}$ is small. This indicates that for the two-directional retrieval scheme with an array of a few layers in the transverse directions, the optimal excitation profile is close to a Gaussian profile.

In Fig. 4(c) we show the behavior of the photon retrieval error as a function of ${L_{z}}$ . When $L_{z}$ increases, more of the photonic field can leak out through the transverse directions that are perpendicular to the detection mode. Thus, in general, the photon retrieval efficiencies increase when ${L_{z}}$ increases. For the Gaussian profile, one obtains a power-law scaling, that $c_{{L_{z}}}^{{\text{Gauss}}}\approx L_{z}^{0.7}$ . In the optimal case, the behavior is more involved, since in this case the excitation profile adapts to the detection mode, which partially compensates for the leakage of the photonic field.

The other retrieval scheme we consider is a typical uni-directional retrieval shown in Fig. 4(d). In this case, the detection mode is a single vector Gaussian beam. In Fig. 4(e), we show the scaling of the photon retrieval errors $\varepsilon_{\rm em}^{\rm opt}$ and $\varepsilon_{\rm em}^{\rm Gauss}$ as a function of ${L_{v}}$ . Apart from the region where ${L_{v}}$ is small, $\varepsilon_{\rm em}^{\rm opt}$ and $\varepsilon_{\rm em}^{\rm Gauss}$ do not depend on ${L_{v}}$ , which is expected from the result of the atomic ensemble, since the optical depth does not involve the transverse size. In Fig. 4(f) we show the scaling of $\varepsilon_{\rm em}^{\rm opt}$ and $\varepsilon_{\rm em}^{\rm Gauss}$ as a function of ${L_{z}}$ . We get ${\varepsilon^{\rm opt}_{\rm em}}\approx{\rm{1.45}}L_{z}^{-2}$ and ${\varepsilon^{\rm Gauss}_{\rm em}}\approx{\rm{0.84}}L_{z}^{-1}$ when ${L_{v}}/L_{z}$ is not very small. Here the scaling of $\varepsilon_{\rm em}^{\rm Gauss}$ is already reaching the best scaling for the atomic ensemble case [86, 49]. We further find that in a wide range of geometries $L_{v}$ and $L_{z}$ , the $\varepsilon^{\rm Gauss}_{\rm em}$ of uni-directional scheme can be related to that of the two-directional scheme through

\varepsilon_{\rm em}^{\rm Gauss}\left(\textrm{1Dir}\right)\approx 0.84{L^{-1}_{z}}+\varepsilon_{\rm em}^{\rm Gauss}\left(\textrm{2Dir}\right),

(29)

which indicates that the error ${\varepsilon_{\rm em}^{\rm Gauss}}$ of the uni-directional retrieval scheme stems from both the photonic field leaked from the direction opposite to the detection mode, which scales as $0.84{L^{-1}_{z}}$ , and the field leakage to the transverse directions, which share the same scaling as the error ${\varepsilon_{\rm em}^{\rm Gauss}}$ for the two-directional retrieval scheme. This behavior predicts an increase of $\varepsilon_{\rm em}^{\rm Gauss}$ when ${L_{v}}/{L_{z}}\ll 1$ , where the leakage to transverse direction is the dominant error. Such behavior is observed in Fig. 4(f) for $L_{v}=6$ case, and qualitatively explains the increase of $\varepsilon_{\rm em}^{\rm opt}$ in that region as well.

The scaling $\varepsilon_{\rm em}^{\rm Gauss}\propto L_{z}^{-1}$ of the uni-directional retrieval error with Gaussian profile indicates that it shares the same physics as the photon retrieval process in disordered atomic ensembles where the optical depth ${\rm OD}\propto{L_{z}}$ . Thus, just like the atomic ensemble case, we can add a cavity to enhance the optical depth of the atomic array. Assuming we add a cavity of finesse $F$ , the photon retrieval error is reduced approximated $F$ -fold. Such a cavity-enhanced photon retrieval from an optical lattice has been demonstrated in Ref. [93, 67] with $F\approx 50$ .

Imperfections during the photon emission process include (i) array defects, (ii) thermal atomic motion, and (iii) the deviation of the Holstein-Primakoff approximation. Here we briefly describe their effects, and the details are shown in Appendix E and in Appendix G.

Atomic defects.

As shown in Ref. [88], one expects that the relative decrease of the retrieval efficiency should be proportional to the ratio of the detection mode hitting the empty sites. Numerical results show that the following relation for the retrieval efficiency $p_{\rm em}^{\rm def}$ with defects holds [88]

p_{\rm em}^{\rm def}={p_{\rm em}}\left({1-\alpha_{\rm def}\frac{{\sum\limits_{j\in{\rm{def}}}{{{\left|{{E_{j}}}\right|}^{2}}}}}{{\sum\limits_{l}{{{\left|{{E_{l}}}\right|}^{2}}}}}}\right),

(30)

where ${E_{j}}={{\bf{E}}_{{\rm{det}}}}\left({{{\bf{r}}_{j}}}\right)\cdot{\bf{d}}_{j}^{*}$ is the intensity of the detection mode at the position of the $j$ -th atom, and the coefficient $\alpha_{\rm def}$ depends on $d_{0}$ and the array geometry.

Atomic thermal motion.

The effect of atomic thermal motion is modeled as a random spatial disorder with a standard deviation of $\sigma_{\rm th}$ [94]. This leads to an additional photon retrieval error that scales as $\sigma_{{\rm{th}}}^{2}/{d_{0}^{2}}$ , the same as predicted in Refs. [88, 95].

Holstein-Primakoff approximation [62].

The error due to the approximation [Eq. 11] comes in when our photonic MPS generation protocol in Section II involves states with multiple collective excitations. This error lead to a lower bound on the photon retrieval efficiency for generating MPS with bond dimension $D$ and physical dimension $d$ in the regime of $D+d-2\ll N$ as

{p_{{\rm{em}}}}\to{p^{\prime}_{{\rm{em}}}}\geq{\left({1-\frac{{\left({D+d-2}\right)\left({D+d-3}\right)}}{2N}}\right)^{1/\left({d-1}\right)}}{p_{{\rm{em}}}}.

(31)

For the cluster state, we have $D=2$ and $d=2$ , thus ${p^{\prime}_{{\rm{em}}}}\geq\left({1-\frac{1}{N}}\right){p_{{\rm{em}}}}$ .

V Optimal performance of the scheme

Given that the error per photon ${{\xi_{{\rm{opt}}}}}$ scales as Eq. 25 and with the retrieval efficiencies studied in Section IV, we can optimize the array geometry ${L_{v}}$ and ${L_{z}}$ to find the minimal ${\xi_{{\text{opt}}}}$ .

The array geometry will affect both the retrieval efficiency $p_{\rm em}$ and the Rydberg nonlinear shift $U$ . The $U$ is determined by the van der Waals interaction ${V_{ij}}={C_{6}}/{\left|{{{\vec{r}}_{i}}-{{\vec{r}}_{j}}}\right|^{6}}$ between Rydberg atoms at position ${\vec{r}_{i}}={d_{0}}\vec{i}={d_{0}}\left({{i_{x}},{i_{y}},{i_{z}}}\right)$ and ${\vec{r}_{j}}$ as [96]

U=\sqrt{\frac{{N\left({N-1}\right)}}{{\sum\limits_{\vec{i}\neq\vec{j}}{{{\left|{\vec{i}-\vec{j}}\right|}^{12}}}}}}\frac{{{C_{6}}}}{{d_{0}^{6}}}\equiv{f_{{L_{v}},{L_{z}}}}\frac{{{C_{6}}}}{{d_{0}^{6}}},

(32)

where we denote the geometric factor as ${f_{{L_{v}},{L_{z}}}}$ . Substituting Eq. 32 and $N=L_{v}^{2}{L_{z}}$ into Eq. 25 and using the lower bound of $p_{\rm em}$ in Eq. 31, we get

{\xi_{{\rm{opt}}}}{\rm{}}={\beta_{0}}+{\left({\frac{{27{\beta_{U}}}}{{4f_{{L_{v}},{L_{z}}}^{2}}}}\right)^{1/3}}{\left({\frac{{\Gamma d_{0}^{6}}}{{\left|{{C_{6}}}\right|}}}\right)^{2/3}}{\rm{}}-{\beta_{{\rm{em}}}}\cdot\log{{p^{\prime}_{{\rm{em}}}}}.

(33)

Since in experiments the parameters ${\Gamma},{C_{6}},{d_{0}}$ depend on the particular atomic configuration chosen, in our numerical optimization we scan over $L_{v}$ and $L_{z}$ to find the minimal $\xi_{\rm opt}$ with respect to a given $\left|{{C_{6}}}\right|/\left({\Gamma d_{0}^{6}}\right)$ . To quantify the performance more intuitively, we define an entanglement length ${N_{\rm ph}}\equiv\log 2/\xi_{\rm opt}$ , which is the photon number that can be generated with ${\cal F}_{\rm ph}=1/2$ .

As a demonstration, the result of this optimization for generating cluster state with optimal control pulse in Fig. 2(b) is shown in Fig. 5. For the two-directional retrieval scheme, we obtain ${N_{{\rm{\rm ph}}}}\sim{[{|{{C_{6}}}|/({\Gamma d_{0}^{6}})}]^{0.23}}$ , shown in Fig. 5(a). The corresponding optimal array geometry $L_{z}^{\rm opt}$ and $L_{v}^{\rm opt}$ can also be estimated from a power law scaling $L_{z}^{\rm opt}\sim{[{|{{C_{6}}}|/({\Gamma d_{0}^{6}})}]^{0.05}}$ and $L_{v}^{\rm opt}\sim{[{|{{C_{6}}}|/({\Gamma d_{0}^{6}})}]^{0.1}}$ . We find such power law scaling behavior is universal for both uni-directional and two-directional retrieval schemes, with different exponents depending on the non-universal constants $\{{{\beta_{i}}}\}$ and the specific retrieval scheme. We further compare the performance of different schemes in Fig. 5(b), from which we see the uni-directional retrieval + cavity scheme and two-directional retrieval scheme gives rise to favourable scaling of $N_{\rm ph}$ . To make a quantitative estimation of $N_{\rm ph}$ , we use current state-of-the-art experimental parameters [82], in which the Rydberg state is $|{70S,J=1/2,{m_{J}}=-1/2}\rangle$ , with

\begin{gathered}{C_{6}}={-862}.{\text{69GHz}}\cdot\mu{{\text{m}}^{\text{6}}},\hfill\\ \Gamma_{r}=19.6\text{kHz},\quad{\Gamma_{\phi}}=21.3\text{kHz},\hfill\\ \end{gathered}

(34)

and assume we work on a optical lattice with ${d_{0}}=532\rm{nm}$ [97], which leads to $\left|{{C_{6}}}\right|/\left({\Gamma d_{0}^{6}}\right)=7.9\times{10^{7}}$ . With this parameter set, the uni-directional protocol in free space can generate a cluster state of ${N_{\rm ph}}\approx 9$ photons with fidelity ${{\cal F}_{\rm ph}}=1/2$ . If we add a ring cavity of finesse $F=50$ (similar to that used in Ref. [67]), this improves to ${N_{\rm ph}}\approx 75$ photons. While for the two-directional retrieval scheme, this parameter set yields ${N_{\rm ph}}\approx 47$ photons. We estimate the required maximal Rabi frequency ${\Omega_{\max}}\sim 2\pi\times\rm{15MHz}$ for the this set of parameters, already available in current experiments [98]. Depending on the available laser power, one could synthesize pulses that lead to smaller values of ${\beta_{r}}({{\beta_{\phi}}})/{\beta_{U}}$ compared to that in Fig. 2(b) to reduce ${\Omega_{\max}}$ .

These results indicate that using the uni-directional retrieval + cavity scheme or the two-directional retrieval scheme, our scheme can deterministically generate strongly-entangled photonic MPS with a large number of photons. We further qualitatively estimate that ${N_{\rm ph}}\sim{D^{-2}}{d^{-2}}$ in Appendix H. The entanglement length $N_{\rm ph}$ can be improved by reducing $\Gamma d_{0}^{6}/\left|{{C_{6}}}\right|$ , which can be done by choosing a Rydberg state with a larger principal number, making the interatomic distance $d_{0}$ smaller, or reducing the decoherences of the Rydberg excitation. Using a cavity with larger finesse $F$ also improves $N_{\rm ph}$ . We also point out that our analysis based on the atomic array can be applied to the case of a disordered atomic ensemble as well.

VI Free-space Multi-port Emission

An important ingredient of the quantum network is to create large-scale distributed quantum entanglement. In principle, thanks to the sequential nature of our protocol, we can generate temporally separated photons and send them to a photon switch within an optical fiber to distribute photons into different ports within the fiber. Here we go one step further to propose a setup and protocol that can directly distribute the sequentially generated photons into multi-ports in free space.

In our MPS generation protocol, the direction of photon emission can be controlled by the photon retrieval laser $\rm L_{e}$ , which imprints a plane-wave phase pattern with a specific momentum. Due to the conservation of momentum, the emitted photon tends to fly to a phase-matched direction [49]. Based on this effect, we propose the following setup depicted in Fig. 6(a) for the multi-port photon generation, where we have a one-layer, two-dimensional square array. The laser $\rm{L}_{g}$ is perpendicular to the plane to create collective Rydberg excitation, and the lasers $\rm{L}_{q}$ ( $\rm{L}_{l}$ ) with plane-wave profile with fixed momenta is used to transfer the collective exciation between $\left|r\right\rangle$ and $\left|q\right\rangle$ ( $\left|l\right\rangle$ ). The emission direction of photons is determined by the momentum of the laser $\rm L_{e}$ . The photon(s) will fly to two directions that have the same momentum component parallel to the plane, and opposite momentum component perpendicular to the plane, characterized by the angle $\theta$ . Thus we choose the detection mode to be a symmetric superposition of two vector Gaussian modes along the photon emission directions.

The photon retrieval error of the excitation with Gaussain profile $\varepsilon_{{\rm{em}}}^{{\rm{Gauss}}}$ as a function of $\theta$ is characterized in Fig. 6(b). As expected, $\varepsilon_{{\rm{em}}}^{{\rm{Gauss}}}$ gradually increases with larger angle $\theta$ . The emitted photonic field along each direction is well characterized by a Gaussian beam [49, 88] with wavelength $\lambda_{eg}$ , beam waist $w_{0}$ , and asymptotic beam angle $\theta_{0}=\lambda_{eg}/{\pi{w_{0}}}$ . In order to make the photons in different ports distinguishable, the angular distance between different ports should at least be $2{\theta_{0}}$ . The minimal version of the multi-port device consist of two ports with angles ${\theta_{1}}={\theta_{0}}$ and ${\theta_{2}}=-{\theta_{0}}$ . Thus we substitute the retrieval efficiency $p_{\rm em}$ at the angle $\theta={\theta_{0}}$ into Eq. 33 and optimize $L_{v}$ to get the maximal entanglement length $N_{\rm ph}$ . The scaling of $N_{\rm ph}$ for the two-port device is shown in Fig. 5(b), and we get ${N_{{\rm{\rm ph}}}}\sim{[{|{{C_{6}}}|/({\Gamma d_{0}^{6}})}]^{0.24}}$ for this case, with corresponding optimal array size $L_{v}^{{\rm{opt}}}\sim{[{|{{C_{6}}}|/({\Gamma d_{0}^{6}})}]^{-0.1}}$ (data not shown). With the current experimental parameters [Eq. 34], we can generate cluster state of $N_{\rm ph}\approx 28$ with this two-port device. By allowing photon emission with larger angle $\theta$ , one can increase the number of ports available at the expense of increasing photon retrieval error. We also point out that such multi-port photon retrieval works for three-dimensional atomic arrays as well [49], where the emission direction is further affected by the Bragg scattering along the $z$ direction in Fig. 6(a).

VII Experimental Consideration

The level scheme in Fig. 1(a) can be realized in various types of atoms. In Fig. 7 we illustrate a possible level scheme based on rubidium-87 ( ${}^{87}\textrm{Rb}$ ), with quantum number of relevant states marked on the figure. The Rydberg state could be $\left|r\right\rangle=\left|{70S,J=1/2,{m_{J}}=1/2}\right\rangle$ , for example.

The experiment can be initialized by preparing all atoms in the state ${\left|g\right\rangle^{\otimes N}}$ through optical pumping. Each state can be coupled to the Rydberg state with two-photon transitions (shown as arrows with different colors in Fig. 7). By choosing the frequencies and the polarizations of different lasers, one can address the transitions depicted in Fig. 7 individually. As a particular example, one can directly couple the state $\left|l\right\rangle$ and $\left|q\right\rangle$ by two-photon Raman transition using the lasers $\rm{L_{l}}$ and $\rm{L_{q}}$ shown in Fig. 7. This allows us to simply implement the $U_{\left[n\right]}^{\rm cl}$ and $U_{\left[n\right]}^{\rm gG}$ discussed in Section III.

In Fig. 7, the on-demand photon emission is realized by a two-photon Raman transition that couples $\left|{{l}}\right\rangle\to\left|{{e}}\right\rangle$ through the intermediate state $\left|{{g}}\right\rangle$ with a detuning ${\Delta}$ [49]. Due to the selection rules on $m_{F}$ , the photon decays from $\left|{{e}}\right\rangle$ only to $\left|{{g}}\right\rangle$ . In the case of generating photonic MPS with physical dimension $d=2$ , at most a single photon is emitted each time. In that case, one can alternatively use a laser to couple state $\left|r\right\rangle$ and $\left|e\right\rangle$ , go through the process $\left|l\right\rangle\to\left|r\right\rangle\to\left|e\right\rangle$ to emit photons. In this case it is further possible to operate the protocol without using the state $\left|l\right\rangle$ , which could simplify the experimental realization, at the expense of storing the collective excitation on $\left|r\right\rangle$ that has a larger decoherence rate.

VIII Summary and Outlook

Summing up, we have proposed a physical platform and a protocol to deterministically generate multiphoton entangled states in free space based on a Rydberg-blockaded atomic array. We have developed both a quantum gate approach and the quantum optimal control approach to universally control a large Hilbert space, which enables the generation of photonic MPS with large bond and physical dimensions, and allows one to experimentally explore the rich regime of one-dimensional spin systems. The strong nonlinearity and long coherence time of the Rydberg state together with the good photon retrieval efficiencies of atomic arrays, lead to favorable scaling of the entanglement length of the photonic MPS, which may enable the generation of large-scale photonic quantum entanglement beyond the state of the art. Furthermore, one can control the photon emission angle of the sequentially generated photons, which makes this platform a potential multi-port device that can efficiently distribute entangled photons in free-space.

This work can be extended in many ways. First, one can further reduce the decoherence errors by using long-lived states in alkaline earth atoms [99], optimize the pulse shape, and apply open system optimal control techniques [100, 79, 101]. Second, a more involved atomic level scheme may allow one to produce photonic MPS with polarization encoding [49, 102], and the photon directional retrieval enables momentum multiplexing [103, 49], which can be used to generate higher-dimensional photonic tensor network states in free space [102]. The ability to generate tensor network states with high bond dimension and physical dimension further call for identifying genuine quantum resource states in this class of states that can be efficiently created in the near-term devices. Finally, the principle to obtain a large Hilbert space with bosonic modes and the associated control techniques developed here can be extended to various other platforms, for example waveguide QED [6] and circuit QED systems [104].

Acknowledgements.

We thank Jun Rui and Tao Shi for insightful discussions. Z.-Y.W, D.M., and J.I.C. acknowledge funding from ERC Advanced Grant QENOCOBA under the EU Horizon 2020 program (Grant Agreement No. 742102) and within the D-A-CH Lead-Agency Agreement through project No. 414325145 (BEYOND C). AGT acknowledges support from CSIC Research Platform on Quantum Technologies PTI001 and from Spanish project PGC2018-094792-B-100 (MCIU/AEI/FEDER, EU).

Appendix A Generation of ancilla-photon superposed states

In this appendix, we show that an ancilla-emitter system is further suitable for generating superposed states between the ancilla and the multi-photon MPS [cf. Eq. 5]. This is a specific class of the superposed states considered in Ref. [60], where many applications of such states are discussed. For example, one can generate a state of the form [60]:

|\psi\rangle=\frac{1}{{\sqrt{2}}}|0{\rangle_{m}}|{\rm{GH}}{{\rm{Z}}_{\rm{1}}}\rangle_{\rm ph}+\frac{1}{{\sqrt{2}}}|1{\rangle_{m}}\left|{{\rm{GH}}{{\rm{Z}}_{2}}}\right\rangle_{\rm ph},

(35)

where $\left|{{\rm{GH}}{{\rm{Z}}_{1}}}\right\rangle_{\rm ph}=\frac{1}{{\sqrt{2}}}\left({\left|{0000}\right\rangle_{\rm ph}-\left|{1111}\right\rangle_{\rm ph}}\right)$ and $\left|{{\rm{GH}}{{\rm{Z}}_{2}}}\right\rangle=\frac{1}{{\sqrt{2}}}\left({\left|{0011}\right\rangle_{\rm ph}+\left|{1100}\right\rangle_{\rm ph}}\right)$ . Such a coherent superposition allows one to choose the entanglement structure of the photons by measuring the ancilla in different basis.

To create the states of the form Eq. 5, let us consider adding an $m$ -level ancilla to the MPS source discussed in Section II.1. So the system Hilbert space becomes ${{\cal H}_{m}}\otimes{{\cal H}_{D}}\otimes{{\cal H}_{d}}$ . As pointed out in Section II.1, the $m$ -level ancilla is not necessarily a physical system. For example, in the system described in the main text (cf. Section II.2), the role of the $m$ -level ancilla can be played by a part of the Hilbert space ${{\cal H}_{q}}$ . In this case, one needs $\dim{{\cal H}_{q}}=4$ to act as two 2-dimensional ancillas, to create the state of form Eq. 35.

Starting from an initial state

\left|{{\psi_{I}}}\right\rangle={\sum\limits_{i=0}^{m-1}{{\alpha_{i}}{{\left|i\right\rangle}_{m}}{{\left|{{\varphi_{I,i}}}\right\rangle}_{D}}}}{\left|0\right\rangle_{d}},

(36)

where ${\left\{{\left|{{\varphi_{I,i}}}\right\rangle}\right\}_{i=0,...,D-1}}$ correspond to the boundary condition of photonic MPS, in the $s$ -th round of photon generation we apply the unitary $\hat{U}\left[s\right]$ of following form

\hat{U}\left[s\right]=\sum\limits_{ijj^{\prime}kk^{\prime}}{U_{jk}^{j^{\prime}k^{\prime}}}\left[{i,s}\right]({\left|i\right\rangle_{m}}\langle i|)\otimes({\left|{j^{\prime}}\right\rangle_{D}}\langle j|)({\left|{k^{\prime}}\right\rangle_{d}}\langle k|),

(37)

where due to the unitarity of $\hat{U}\left[s\right]$ ,

\sum\limits_{p,q}{U_{pq}^{*j^{\prime}k^{\prime}}\left[{i^{\prime},s}\right]U_{jk}^{pq}\left[{i,s}\right]}={\delta_{ii^{\prime}}}{\delta_{jj^{\prime}}}{\delta_{kk^{\prime}}},\quad\forall s=1,...,n.

(38)

After the unitary operation, we trigger a photon emission [cf. Eq. 1]. After applying this sequence for $n$ rounds, the state becomes

\begin{array}[]{l}\left|{{\psi_{F}}}\right\rangle={M_{P}}\hat{U}\left[n\right]...{M_{P}}\hat{U}\left[1\right]\left|{{\psi_{I}}}\right\rangle\\ =\sum\limits_{i=0}^{m-1}{{\alpha_{i}}{{\left|i\right\rangle}_{m}}\left({U_{\left[{i,n}\right]}^{{k_{n}}}...U_{\left[{i,1}\right]}^{{k_{1}}}{{\left|{{\varphi_{I,i}}}\right\rangle}_{D}}}\right){{\left|0\right\rangle}_{d}}}{\left|{{k_{n}}...{k_{1}}}\right\rangle_{{\rm{\rm ph}}}},\end{array}

(39)

where ${\left({U_{\left[{i,s}\right]}^{{k_{1}}}}\right)_{jj^{\prime}}}\equiv U_{j0}^{j^{\prime}{k_{1}}}\left[{i,s}\right]$ . In the last step we can disentangle the $D$ -level ancilla with the photons and the $m$ -level ancilla [14]. For example, we can map the $D$ -level ancilla to a state $\left|{{\varphi_{F}}}\right\rangle_{D}$ at last. After tracing out the $D$ -level ancilla and the $d$ -level emitter, we get the desired state [cf. Eq. 5] of the $m$ -level ancilla and the photons, with

\left|{\psi_{{\rm{MPS}}}^{i}}\right\rangle_{\rm ph}=\sum\limits_{{k_{n}}...{k_{1}}}{{}_{D}\left\langle{{\varphi_{F}}}\right|U_{\left[{i,n}\right]}^{{k_{n}}}...U_{\left[{i,1}\right]}^{{k_{1}}}{{\left|{{\varphi_{I,i}}}\right\rangle}_{D}}{{\left|{{k_{n}}...{k_{1}}}\right\rangle}_{{\rm{ph}}}}}.

(40)

Here the unitarity of ${\left\{{\hat{U}\left[s\right]}\right\}_{s=1,...,n}}$ [cf. Eq. 38] lead to the canonical form of ${\left\{{\left|{\psi_{{\rm{MPS}}}^{i}}\right\rangle_{\rm ph}}\right\}_{i=0,...,\left({m-1}\right)}}$ .

The above process can also be viewed as a quantum random access memory [61], where the input quantum address is the initial state Eq. 36, and each (quantum) memory cell is a multi-photon MPS.

Appendix B The quantum gate approach

Here, we present the details of the quantum gate approach for generating $D=2,d=2$ photonic MPS. As discussed in Section III.1, the key step is to implement arbitrary unitaries in ${{\cal H}_{\rm mid}}={\rm span}({|0\rangle,|{{1_{r}}}\rangle,|{{1_{q}}}\rangle,|{{1_{r}}{1_{q}}}\rangle})$ . The resonant driving $\mathrm{L_{g}}$ can be decomposed into the $X$ and $Y$ component, with

\begin{array}[]{l}X\left(\theta\right)=\exp\left(-{i\left({{I_{q}}\otimes{\sigma_{x}}}\right)\theta/2}\right),\\ Y\left(\theta\right)=\exp\left(-{i\left({I_{q}\otimes{\sigma_{y}}}\right)\theta/2}\right),\end{array}

(41)

and we can construct $Z(\theta)=X({-\pi/2})Y({-\theta})X({\pi/2})$ . On the other hand, the resonant driving $\mathrm{L_{q}}$ can be decomposed into

\begin{array}[]{l}{X_{q}}\left(\theta\right)=\exp\left[{-i\theta\left({{a_{q}}\otimes\sigma_{r}^{\dagger}+a_{q}^{\dagger}\otimes{\sigma_{r}}}\right)/2}\right],\\ {Y_{q}}\left(\theta\right)=\exp\left[{\theta\left({a_{q}^{\dagger}\otimes{\sigma_{r}}-{a_{q}}\otimes\sigma_{r}^{\dagger}}\right)/2}\right],\end{array}

(42)

which couples the states $\left|{{1_{r}}}\right\rangle$ and $\left|{{1_{q}}}\right\rangle$ , but also couples $\left|{{1_{r}}{1_{q}}}\right\rangle$ to $\left|{{2_{q}}}\right\rangle$ , which is out of ${{\cal H}_{\rm mid}}$ . To avoid this leakage, we choose to construct a SWAP gate and the CNOT gate within $\mathcal{H}_{\rm mid}$ without population leakage to state $\left|{{2_{q}}}\right\rangle$ .

Let us apply ${\rm L}_{q}$ with Rabi frequency ${\Omega_{rq}}$ and a large detuning $\Delta_{\rm ST}\gg{\Omega_{rq}}$ . This induces different AC Stark shifts between $\left|{{1_{r}}}\right\rangle$ and $\left|{{1_{q}}}\right\rangle$ , and between $\left|{{1_{r}1_{q}}}\right\rangle$ and $\left|{{2_{q}}}\right\rangle$ :

\left|{{1_{r}}}\right\rangle\leftrightarrow\left|{{1_{q}}}\right\rangle\mathrel{\mathop{\mathchar 58\relax}}\frac{{\Omega_{rq}^{2}}}{{2\Delta}},\qquad\left|{{1_{r}}{1_{q}}}\right\rangle\leftrightarrow\left|{{2_{q}}}\right\rangle\mathrel{\mathop{\mathchar 58\relax}}\frac{{\Omega_{rq}^{2}}}{{\Delta}}.

(43)

By applying the pulse for a time $t=4{\Delta_{{\text{ST}}}}\theta/\Omega_{{\text{ST}}}^{2}$ , we obtain a phase gate ${S_{T}}\left(\theta\right)={\rm{diag}}\left({1,{e^{-i\theta}},{e^{i\theta}},{e^{-2i\theta}},{e^{2i\theta}}}\right)$ with the basis $\{{\left|0\right\rangle,\left|{{1_{r}}}\right\rangle,\left|{{1_{q}}}\right\rangle,\left|{{1_{r}}{1_{q}}}\right\rangle,\left|{{2_{q}}}\right\rangle}\}$ . Using ${S_{T}}\left(\theta\right)$ and resonant laser drivings, the SWAP and CNOT gates are constructed as:

	$\displaystyle{\rm{SWAP}}$	$\displaystyle={Y_{q}}\left({-\pi/4}\right)Z\left({\pi/2}\right){S_{T}}\left({\pi/2}\right)\hfill$		(44)
		$\displaystyle\cdot{X_{q}}\left({\pi/4}\right){S_{T}}\left({-\pi/2}\right)Z\left({-\pi/2}\right),\hfill$		(44)

{\rm{CNOT=}}X\left({5\pi/16}\right)Z\left({\pi/2}\right){S_{T}}\left({\pi/2}\right)X\left({3\pi/16}\right).

(45)

These two gates together with the single-qubit rotations produced by $\mathrm{L}_{\rm g}$ constitute a universal two-qubits gate set in $\mathcal{H}_{\rm mid}$ without leakage error.

We know each unitary in our protocol should act as Eq. 2. For the case of generating MPS with $D=2$ and $d=2$ , this is equivalent to implementing the isometry ${\hat{V}_{\left[i\right]}}\mathrel{\mathop{\mathchar 58\relax}}{{\cal H}_{D=2}}\to{\cal H}_{{\rm{src}}}^{\left({2,2}\right)}$ . One can synthesize a quantum circuit to implement desired isometries, for example using the UniversalQCompiler package [105].

As a demonstration, we construct the unitarie that generates the cluster state [Eq. 18]. Using a SVD decomposition approach [14], we find two isometries ${{\hat{V}}^{\rm cl}_{\left[{i\neq n}\right]}}$ and ${{\hat{V}}^{\rm cl}_{\left[{n}\right]}}$ needed to generate $n$ -photon cluster state

\begin{gathered}{{\hat{V}}^{\rm cl}_{\left[{i\neq n}\right]}}=\left({\begin{array}[]{*{20}{c}}{V_{\left[i\right]}^{0}}\\ {V_{\left[i\right]}^{1}}\end{array}}\right)=\frac{1}{{\sqrt{2}}}\left({\begin{array}[]{*{20}{c}}1&0\\ 1&0\\ 0&1\\ 0&{-1}\end{array}}\right),\qquad\hfill\\ {{\hat{V}}^{\rm cl}_{\left[n\right]}}=\left({\begin{array}[]{*{20}{c}}{V_{\left[n\right]}^{0}}\\ {V_{\left[n\right]}^{1}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}1&0\\ 0&0\\ 0&1\\ 0&0\end{array}}\right).\hfill\\ \end{gathered}

(46)

The corresponding unitaries ${{\hat{U}}^{\rm cl}_{\left[{i\neq n}\right]}}$ and ${{\hat{U}}^{\rm cl}_{\left[{n}\right]}}$ have the form of Eq. 19. The gate sequences for implementing ${{U}_{\left[i\neq n\right]}}$ and ${{U}_{\left[n\right]}}$ are shown in Fig. 2, where the gates in the second qubit are implemented with $X,Y,Z$ , while the first-qubit rotations are implemented as ${\rm{SWAP}}\to{\textrm{rotations on the second qubit}}\to{\rm{SWAP}}$ .

Appendix C The quantum optimal control (QOC) approach

In QOC we aim to find the time-dependence of ${\left\{{{\Omega_{r\beta}}\left(t\right)}\right\}_{\beta=g,l,q}}$ in $H_{\rm col}(t)$ [Eq. 13] that yields the desired $U_{\rm targ}\in\mathcal{H}_{\rm col}$ of the form Eq. 19, with the embedded isometry ${\hat{V}_{\rm targ}}$ . First, we assume that ${\left\{{{\Omega_{r\beta}}\left(t\right)}\right\}_{\beta=g,l,q}}$ are analytical functions parameterized by a set of parameters $\bar{\alpha}$ . The dynamics of our system under the control Hamiltonian $H_{\rm col}(\bar{\alpha},t)$ gives rise to the time-ordered $(\mathcal{T}$ ) evolution operator (with $\hbar=1$ )

U\left({\bar{\alpha},T}\right)={\cal T}\exp\left(-i{\int_{0}^{T}{H_{\rm col}(\bar{\alpha},t)dt}}\right).

(47)

In numerical calculation, we truncate the dimensions of ${{\cal H}_{l}}$ and ${{\cal H}_{q}}$ such that $\dim{U_{\rm targ}}={N_{h}}$ . The zero matrix $O$ thus is of dimension $\left({{N_{h}}-d\cdot D}\right)\times D$ . We can denote the unitary generated by our QOC pulse with the same block structure as in Eq. 19

U\left({\bar{\alpha},T}\right)=\left({\begin{array}[]{*{20}{c}}{V^{\prime}\left({\bar{\alpha},T}\right)}&{{{B^{\prime}_{1}}}\left({\bar{\alpha},T}\right)}\\ {O^{\prime}\left({\bar{\alpha},T}\right)}&{{{B^{\prime}_{2}}}\left({\bar{\alpha},T}\right)}\end{array}}\right).

(48)

We aim to get $V^{\prime}\left({\bar{\alpha},T}\right)={\hat{V}_{\rm targ}}$ and $O^{\prime}\left({\bar{\alpha},T}\right)=O$ . Let us define $F_{V}\left({\bar{\alpha}}\right)={\text{Tr}}\left({{{\hat{V}_{\rm targ}}^{\dagger}}V^{\prime}\left({\bar{\alpha},T}\right)}\right)$ and $F_{O}\left({\bar{\alpha}}\right)={\text{Tr}}\left({O{{\left({\bar{\alpha},T}\right)}^{\dagger}}O^{\prime}\left({\bar{\alpha},T}\right)}\right)$ . The deviation of $V^{\prime}\left({\bar{\alpha},T}\right)$ from $\hat{V}_{\rm targ}$ is quantified by $(1-\left|{{F_{V}}\left({\bar{\alpha}}\right)}\right|/D)$ and the deviation of $O^{\prime}\left({\bar{\alpha},T}\right)$ from $O$ can be quantified by $\frac{w}{D}\left|{{F_{O}}\left({\bar{\alpha}}\right)}\right|$ with a tunable penalty constant $w$ . Thus our optimization task is to minimize the cost function

g\left({\bar{\alpha}}\right)=1-\frac{1}{D}\left|{{F_{V}}\left({\bar{\alpha}}\right)}\right|+\frac{w}{D}\left|{{F_{O}}\left({\bar{\alpha}}\right)}\right|.

(49)

We can evaluate the gradient of $g\left({\bar{\alpha}}\right)$ as

\begin{array}[]{*{20}{l}}{{\partial_{\bar{\alpha}}}g\left({\bar{\alpha}}\right)=-\frac{1}{{D\left|{{F_{V}}}\right|}}{\text{Re}}\left[{F_{V}^{*}\cdot{\text{Tr}}\left({\hat{V}_{{\text{targ}}}^{\dagger}{\partial_{\bar{\alpha}}}V^{\prime}}\right)}\right]}\\ {+\frac{w}{{D\left|{{F_{O}}}\right|}}{\text{Re}}\left[{F_{O}^{*}\cdot{\text{Tr}}\left[{\left({{\partial_{\bar{\alpha}}}{{O^{\prime}}^{\dagger}}}\right)O^{\prime}+{{O^{\prime}}^{\dagger}}{\partial_{\bar{\alpha}}}O^{\prime}}\right]}\right],}\end{array}

(50)

where ${\partial_{\bar{\alpha}}}V^{\prime}\left({\bar{\alpha},T}\right)$ and ${\partial_{\bar{\alpha}}}O^{\prime}\left({\bar{\alpha},T}\right)$ are the corresponding blocks [cf. Eq. 19] of the matrix ${\partial_{\bar{\alpha}}}U\left({\bar{\alpha},T}\right)$ . The $U\left({\bar{\alpha},T}\right)$ and ${\partial_{\bar{\alpha}}}U\left({\bar{\alpha},T}\right)$ can be obtained using the GOAT algorithm [79] that numerically integrates the coupled system of the equation of motion

{\partial_{t}}\left({\begin{array}[]{*{20}{c}}U\\ {{\partial_{\bar{\alpha}}}U}\end{array}}\right)=-i\left({\begin{array}[]{*{20}{c}}H&0\\ {{\partial_{\bar{\alpha}}}H}&H\end{array}}\right)\left({\begin{array}[]{*{20}{c}}U\\ {{\partial_{\bar{\alpha}}}U}\end{array}}\right).

(51)

With the knowledge of $g\left({\bar{\alpha}}\right)$ and ${\partial_{\bar{\alpha}}}g\left({\bar{\alpha}}\right)$ , the control pulse can be optimized by typical gradient descent search methods. In our numerical optimization, we synthesize individual pulses for the real part and the imaginary part for each ${\left\{{{\Omega_{r\beta}}\left(t\right)}\right\}_{\beta=g,l,q}}$ , since this ensures a real value of each pulse amplitude. For each individual component, we choose a Fourier form ${f}\left({\bar{\alpha},t}\right)=\sum\nolimits_{j=1}^{{j_{\max}}}{{A_{j}}\sin\left({j\cdot{w}\cdot t}\right)}$ as our pulse amplitude, where the parameters are $\bar{\alpha}={\left\{{{A_{j}},{w}}\right\}_{j=1,...,{j_{\max}}}}$ . We also add two sigmoid functions to bound the amplitude and enforce a smooth start and end of the pulse

\begin{gathered}{S_{1}}\left({f\left({\bar{\alpha},t}\right)}\right)=b\left({1-\frac{2}{{1+\exp\left({{g_{1}}\cdot f\left({\bar{\alpha},t}\right)/b}\right)}}}\right),\hfill\\ {S_{2}}\left({t,T}\right)=1-\frac{2}{{1+{e^{{g_{2}}\left({T-t}\right)}}}},\hfill\\ \end{gathered}

(52)

where $b,g_{1}$ and $g_{2}$ are tuning parameters. The overall pulse form for each individual component is ${S_{1}}\left({b,{g_{f}},f\left({\bar{\alpha},t}\right)}\right)\cdot{S_{2}}\left({t,T}\right)$ .

As an example, the synthesized pulse for implementing $U_{\left[{i\neq n}\right]}^{\rm cl}$ (with the corresponding isometry ${\hat{V}}_{\left[{i\neq n}\right]}^{\rm cl}$ [Eq. 46]) for the cluster state generation is shown in Fig. 2, and the pulse for implementing $U_{\left[{i\neq n}\right]}^{\rm gG}$ for generating the generalized GHZ state Eq. 22 of physical dimension $d=3$ is shown in Fig. 8.

We can quantify the quality of the pulses in Fig. 2(b) and in Fig. 8 by computing the cost function without the $O$ matrix part as ${g_{V}}\left({{U_{{\rm{targ}}}}}\right)={\left.{g\left({\bar{\alpha}}\right)}\right|_{w=0}}$ , and getting $g_{V}\left({U_{\left[{i\neq n}\right]}^{\rm cl}}\right)={\rm{1.1}}\times{\rm{1}}{{\rm{0}}^{{\rm{-5}}}}$ , $g_{V}\left({U_{\left[{i\neq n}\right]}^{\rm gG}}\right)=2\times{10^{-3}}$ , which clearly demonstrate the effectiveness of this method.

Appendix D The effective modeling of long-time evolution for a Rydberg-blockaded atomic array

To compute the photonic state fidelity ${{\cal F}_{{\rm{\rm ph}}}}=\left\langle{{\psi_{{\rm{MPS}}}}}\right|{\rho_{{\rm{\rm ph}}}}\left|{{\psi_{{\rm{MPS}}}}}\right\rangle$ efficiently, it will be shown in Appendix F that we need to numerically obtain the full Liouvillian propagator ${W_{\cal L}}$ [Eq. 98] for the evolution process, which is a computationally demanding task, as its size grows with the Hilbert space dimension of the atomic array to the power of 4. Thus, directly solving the master equation Eq. 23 in the many-body Hilbert space ${\cal H}$ is only possible for array sizes of tens of atoms. To treat the case of larger atom numbers, we need to find an effective master equation that captures the dynamics of the system under the evolution of Eq. 23. For this purpose, we shall first understand the effects of various imperfections discussed in Section III.3.

Recall the basis elements of the collective Hilbert space ${{\cal H}_{\rm col}}$ [cf. Eq. 12]

\left|{{s_{r}}{m_{q}}{n_{l}}}\right\rangle=\frac{1}{{\sqrt{{m_{q}}!{n_{l}}!}}}{({\sigma_{r}^{\dagger}})^{{s_{R}}}}{({a_{q}^{\dagger}})^{{m_{q}}}}{({a_{l}^{\dagger}})^{{n_{l}}}}\left|0\right\rangle,

(53)

where $s_{r}=0,1$ and $m_{q},n_{l}$ are positive integers. First, we include the double Rydberg excitations due to finite van der Waals interaction between Rydberg atoms ${V_{ij}}={C_{6}}/r_{ij}^{6}$ . We model its effect with an effectively uniform level shift $U$ for Rydberg doubly excited states, which can be determined from Eq. 32). Thus, the Rydberg double-excitation can be taken into consideration by adding basis states with $s_{r}=2$ into ${\mathcal{H}_{r}}$ . Thus now the Rydberg collective excitations form a qutrit with the same properties as a truncated anharmonic oscillator with nonlinearity $U$ .

The Rydberg state $\left|{{r}}\right\rangle$ can dephase at a rate ${\Gamma_{\phi}}$ , and can decay to the three hyperfine ground states $\left|g\right\rangle,\left|q\right\rangle,\left|l\right\rangle$ with decay rate ${\Gamma_{r}}$ for each channel. Since the Rydberg decay transition has short wavelength (e.g. $\lambda_{\rm ryd}=294$ nm in Ref. [82]), for a normal ${}^{87}\mathrm{Rb}$ optical lattice with lattice constant ${d_{0}}=532$ nm [97] we have ${d_{0}}/{\lambda_{\rm ryd}}\approx 1.8$ , and therefore the dipole rescattering effect during the Rydberg decay process is very weak.

To get insight for constructing the effective master equation that captures the evolution under Eq. 23, first we analytically solve the decoherence dynamics under Eq. 23 with an initial state $\left|{{1_{r}}}\right\rangle=a_{r}^{\dagger}\left|0\right\rangle=\sum\nolimits_{i=1}^{N}{{u_{i}}{e^{i{{\bf{k}}_{rg}}\cdot{{\bf{r}}_{i}}}}\sigma_{rg}^{i}}\left|0\right\rangle$ by unravelling Eq. 23 into a non-Hermitian evolution $S\left({t,{t_{0}}}\right)\rho={e^{-i{H_{{\rm{nh}}}}t}}\otimes{e^{i{{\bar{H}}_{{\rm{nh}}}}t}}\rho$ and the quantum jumps $J$ [106] as

\begin{array}[]{l}\rho(t)=S\left({t,{t_{0}}}\right)\rho\left({{t_{0}}}\right)\\ +\sum\limits_{n=1}^{\infty}{\int_{{t_{0}}}^{t}{\rm{d}}}{t_{n}}\cdots\int_{{t_{0}}}^{{t_{2}}}{\rm{d}}{t_{1}}S\left({t,{t_{n}}}\right)J\cdots JS\left({{t_{1}},{t_{0}}}\right)\rho\left({{t_{0}}}\right),\end{array}

(54)

with the non-Hermitian Hamiltonian

{H_{{\rm{nh}}}}=H^{\prime}\left(t\right)-\frac{i}{2}\left({3{\Gamma_{r}}+{\Gamma_{\phi}}}\right)\sum\limits_{i=1}^{N}{\sigma_{rr}^{i}},

(55)

and the quantum jump given by

J\rho=\sum\limits_{i=1}^{N}{\left[{{\Gamma_{r}}\sum\limits_{\alpha=\left\{{g,q,l}\right\}}{\sigma_{\alpha r}^{i}\otimes\bar{\sigma}_{\alpha r}^{i}}+{\Gamma_{\phi}}\sigma_{rr}^{i}\otimes\bar{\sigma}_{rr}^{i}}\right]}\rho,

(56)

where $H^{\prime}\left(t\right)$ is same as in Eq. 23. For demonstration, we show the solution under purely spontaneous decay (choosing ${\Gamma_{\phi}}=0$ ) and under purely dephasing (choosing ${\Gamma_{r}}=0$ )

\begin{array}[]{l}\textrm{purely decay}\mathrel{\mathop{\mathchar 58\relax}}\qquad\rho\left(t\right)=\left({{e^{-3{\Gamma_{r}}t}}a_{r}^{\dagger}\otimes\bar{a}_{r}^{\dagger}+\frac{1}{3}\left({1-{e^{-3{\Gamma_{r}}t}}}\right)\left({I+\rho_{q,q}^{\left(1\right)}+\rho_{l,l}^{\left(1\right)}}\right)}\right)\left|{0\otimes\bar{0}}\right\rangle,\\ \textrm{purely dephasing}\mathrel{\mathop{\mathchar 58\relax}}\qquad\rho\left(t\right)=\left({{e^{-{\Gamma_{\phi}}t}}a_{r}^{\dagger}\otimes\bar{a}_{r}^{\dagger}+\left({1-{e^{-{\Gamma_{\phi}}t}}}\right)\rho_{r,r}^{\left(1\right)}}\right)\left|{0\otimes\bar{0}}\right\rangle,\end{array}

(57)

where we have defined the identity operator $I$ and the density matrix creation operator for singly-excited mixed state $\rho_{\alpha,\alpha^{\prime}}^{\left(1\right)}$ as

\rho_{\alpha,\alpha^{\prime}}^{\left(1\right)}=\sum\limits_{i=1}^{N}{{{\left|{{u_{i}}}\right|}^{2}}\sigma_{\alpha g}^{i}\otimes\bar{\sigma}_{\alpha^{\prime}g}^{i}},\qquad\alpha,\alpha^{\prime}=r,q,l.

(58)

In order to get a simplified description for $\rho_{\alpha,\alpha^{\prime}}^{\left(1\right)}$ , we introduce a set of collective creation operators as

a_{{M_{\alpha}}}^{\dagger}=\sum\limits_{i=1}^{N}{{u_{i}}\sigma_{\alpha g}^{i}},\qquad\alpha=r,l,q.

(59)

Here ${\{{a_{{M_{\alpha}}}^{\dagger}}\}_{\alpha=l,q}}$ are oscillator creation operators, and $a_{{M_{r}}}^{\dagger}$ is truncated up to the third Fock state since we include up to double Rydberg excitations. We introduce a projection operator $\mathcal{P}_{\rm diag}$ whose effect is to discard correlation terms with $i\neq j$ . Thus we can simplify the description of $\rho_{\alpha,\alpha^{\prime}}^{\left(1\right)}$ by rewriting it as a projection on the collective excitations as

\begin{array}[]{l}\rho_{\alpha,\alpha^{\prime}}^{\left(1\right)}={{\cal P}_{{\rm{diag}}}}\left[{\sum\limits_{i,j=1}^{N}{\left({{u_{i}}\sigma_{\alpha g}^{i}}\right)\otimes\left({{{\bar{u}}_{j}}\bar{\sigma}_{\alpha^{\prime}g}^{j}}\right)}}\right]\\ ={{\cal P}_{{\rm{diag}}}}\left[{a_{{M_{\alpha}}}^{\dagger}\otimes\bar{a}_{{M_{\alpha}}}^{\dagger}}\right].\end{array}

(60)

Similarly, we can analytically solve the dynamics starting from a doubly excited state $\left|{{2_{r}}}\right\rangle=\frac{1}{{\sqrt{2}}}{\left({a_{r}^{\dagger}}\right)^{2}}\left|0\right\rangle$ . The solution for the purely decay case ( ${\Gamma_{\phi}}=0$ ) is $\rho\left(t\right)={\rho^{\left(0\right)}}\left(t\right)+{\rho^{\left(1\right)}}\left(t\right)+{\rho^{\left(2\right)}}\left(t\right)$ , with

\begin{array}[]{l}{\rho^{\left(0\right)}}\left(t\right)={e^{-6{\Gamma_{r}}t}}\frac{1}{{\sqrt{2}}}{\left({a_{r}^{\dagger}}\right)^{2}}\otimes{\left({\bar{a}_{r}^{\dagger}}\right)^{2}}\left|{0\otimes\bar{0}}\right\rangle,\\ {\rho^{\left(1\right)}}\left(t\right)=\frac{2}{3}\left({{e^{-3{\Gamma_{r}}t}}-{e^{-6{\Gamma_{r}}t}}}\right)\left({a_{r}^{\dagger}\otimes\bar{a}_{r}^{\dagger}}\right)\cdot\left[{I+\rho_{q,q}^{\left(1\right)}+\rho_{l,l}^{\left(1\right)}}\right]\left|{0\otimes\bar{0}}\right\rangle,\\ {\rho^{\left(2\right)}}\left(t\right)=\frac{1}{9}{\left({1-{e^{-3{\Gamma_{r}}t}}}\right)^{2}}\left({{\rho^{\left(0\right)}}+2\rho_{q,q}^{\left(1\right)}+2\rho_{l,l}^{\left(1\right)}+2\rho_{ql,ql}^{\left(2\right)}+\rho_{qq,qq}^{\left(2\right)}+\rho_{ll,ll}^{\left(2\right)}}\right)\left|{0\otimes\bar{0}}\right\rangle,\end{array}

(61)

where we defined a general notation

\rho_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}^{\left({{n_{e}}}\right)}=\sum\limits_{{i_{1}}\neq...\neq{i_{{n_{e}}}}=1}^{N}{\prod\limits_{k=1}^{{n_{e}}}{{{\left|{{u_{{i_{k}}}}}\right|}^{2}}\sigma_{{\alpha_{k}}g}^{{i_{k}}}\otimes\bar{\sigma}_{\alpha^{\prime}_{k}g}^{{i_{k}}}}}.

(62)

Here $\left\{\alpha\right\}\equiv{\alpha_{1}}...{\alpha_{{n_{e}}}},\left\{{\alpha}^{\prime}\right\}\equiv\alpha_{1}^{\prime}...\alpha_{{n_{e}}}^{\prime}$ and ${\alpha_{i}},\alpha_{i}^{\prime}=r,l,q$ . Similar to Eq. 60, Eq. 62 can be rewritten using the $\mathcal{P}_{\rm diag}$ as

\rho_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}^{\left({{n_{e}}}\right)}={{\cal P}_{{\rm{diag}}}}\left[{\prod\limits_{k=1}^{{n_{e}}}{a_{{M_{{\alpha_{k}}}}}^{\dagger}\otimes\bar{a}_{{M_{\alpha^{\prime}_{k}}}}^{\dagger}}}\right].

(63)

Here the effect of $\mathcal{P}_{\rm diag}$ is to discard correlation terms where $\left\{{{i_{1}},...,{i_{n_{e}}}}\right\}\neq\left\{{{j_{1}},...,{j_{n_{e}}}}\right\}$ .

We can also solve the pure-dephasing ( ${\Gamma_{r}}=0$ ) dynamics starting from $\left|{{2_{r}}}\right\rangle$ as

$\displaystyle\rho\left(t\right)$	$\displaystyle=\left[{e^{-2{\Gamma_{\phi}}t}}\frac{1}{{\sqrt{2}}}{\left({a_{r}^{\dagger}}\right)^{2}}\otimes{\left({\bar{a}_{r}^{\dagger}}\right)^{2}}\hfill\right.$	(64)
	$\displaystyle+2t{e^{-2{\Gamma_{\phi}}t}}\left({a_{r}^{\dagger}\otimes\bar{a}_{r}^{\dagger}}\right)\cdot\rho_{r,r}^{\left(1\right)}\hfill$
	$\displaystyle\left.+\left({1-{e^{-2{\Gamma_{\phi}}t}}-2t{e^{-2{\Gamma_{\phi}}t}}}\right)\rho_{rr,rr}^{\left(2\right)}\right]\left\|{0\otimes\bar{0}}\right\rangle.\hfill$

Thus we see that during our sequential MPS generation protocol, the dephasing of $\left|r\right\rangle$ and the decay from $\left|r\right\rangle$ to $\left|l\right\rangle$ and $\left|q\right\rangle$ produce mixed excitations, and these mixed states with $n_{e}$ -excitations are represented by

{{\rho_{\rm mix}}={\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}\rho_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}^{\left({{n_{e}}}\right)}\left|{0\otimes\bar{0}}\right\rangle},

(65)

where ${\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}$ is a normalization factor.

Eq. (63) allows us to greatly simplify the description of $\rho_{\rm mix}$ since now it is characterized by its number of collective excitations. To carefully understand the evolution of ${\rho_{\rm mix}}$ under the master equation Eq. 23, we study its non-Hermitian evolution and the quantum jump one by one. Defining the non-Hermitian evolution operator as

{U_{\rm nh}}\left(t\right)={\cal T}\exp\left({-i\int_{0}^{t}{{H_{\rm nh}}\left(\tau\right)d\tau}}\right),

(66)

the non-Hermitian evolution of $\rho_{\rm mix}$ is

\begin{array}[]{*{20}{l}}{{U_{\rm nh}}\left(t\right)\otimes{{\bar{U}}_{nh}}\left(t\right){\rho_{{\text{mix}}}}={U_{\rm nh}}\left(t\right)\otimes{{\bar{U}}_{nh}}\left(t\right)}\\ \times{\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}{{{\cal P}_{{\text{diag}}}}\left[{\prod\limits_{k=1}^{{n_{e}}}{a_{{M_{{\alpha_{k}}}}}^{\dagger}\otimes\bar{a}_{{M_{{\alpha^{\prime}_{k}}}}}^{\dagger}}}\right]\left|{0\otimes\bar{0}}\right\rangle.}\end{array}

(67)

The effect of ${U_{\rm nh}}\left(t\right)$ can be understood on the single-operator level,

\begin{array}[]{*{20}{l}}{{U_{{\rm{nh}}}}\left(t\right)\sigma_{{\alpha_{1}}g}^{\dagger}...\sigma_{{\alpha_{n}}g}^{\dagger}\left|0\right\rangle=\sigma_{{\alpha_{1}}g}^{\dagger}\left(t\right)...\sigma_{{\alpha_{n}}g}^{\dagger}\left(t\right){U_{{\rm{nh}}}}\left(t\right)\left|0\right\rangle,}\\ {{\text{with}}\quad\sigma_{\alpha g}^{\dagger}\left(t\right)={U_{{\rm{nh}}}}\left(t\right)\sigma_{\alpha g}^{\dagger}U_{\rm nh}^{-1}\left(t\right).}\end{array}

(68)

When $H_{\rm nh}(t)$ does not include the driving $\mathrm{L_{g}}$ , we always have ${U_{\rm nh}}\left(t\right)\left|0\right\rangle=\left|0\right\rangle$ . Also notice that $H_{\rm nh}(t)$ is a collection of single-particle Hamiltonians ${H_{\rm nh}}(t)=\sum\nolimits_{i=1}^{N}{{H_{i}}}(t)$ , such that $\sigma_{\alpha^{\prime}g}^{\dagger}$ with $\alpha^{\prime}\neq\alpha$ does not get mixed into $\sigma_{\alpha g}^{\dagger}\left(t\right)$ . Thus the order of $\mathcal{P}_{\rm diag}$ and $U_{\rm nh}\left(t\right)$ in Eq. 67 can be exchanged as

\begin{array}[]{*{20}{l}}{{U_{\rm nh}}\left(t\right)\otimes{{\bar{U}}_{nh}}\left(t\right){{\cal P}_{{\rm{diag}}}}\left[{\prod\limits_{k=1}^{{n_{e}}}{a_{{M_{{\alpha_{k}}}}}^{\dagger}\otimes\bar{a}_{{M_{{\alpha^{\prime}_{k}}}}}^{\dagger}}}\right]\left|{0\otimes\bar{0}}\right\rangle}\\ {={{\cal P}_{{\rm{diag}}}}\left[{{U_{\rm nh}}\left(t\right)\otimes{{\bar{U}}_{nh}}\left(t\right)\prod\limits_{k=1}^{{n_{e}}}{a_{{M_{{\alpha_{k}}}}}^{\dagger}\otimes\bar{a}_{{M_{{\alpha^{\prime}_{k}}}}}^{\dagger}}}\right]\left|{0\otimes\bar{0}}\right\rangle.}\end{array}

(69)

When $\mathrm{L}_{g}$ is turned on, the non-Hermitian dynamics becomes more involved. We approach it from a different perspective, by looking at the purity $P_{r}$ of ${\rho_{\rm mix}}$ . One gets

P_{r}=\mathrm{Tr}\left[{{\rho^{2}_{\rm mix}}}\right]={\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}\sum\limits_{{i_{1}}\neq...\neq{i_{n_{e}}}}{{{\left|{{u_{{i_{1}}}}}\right|}^{4}}...{{\left|{{u_{{i_{{n_{e}}}}}}}\right|}^{4}}}.

(70)

In our scheme, we consider $\mathrm{L}_{g}$ to be a laser beam that shines on the whole array. Thus its beam profile $u_{i}$ in general scales as $u_{i}\sim O\left({1/\sqrt{N}}\right)$ , and ${\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}\sim O\left({1/{n_{e}}!}\right)$ . Thus ${{P_{r}}}\sim O\left(1/{{n_{e}}!{N^{{n_{e}}}}}\right)\to 0$ when $N\gg 1$ . On the other hand, we know that unitary evolution does not change $P_{r}$ . Thus, starting from $\rho_{\rm mix}$ , when the system is only driven by $\mathrm{L}_{g}$ , for an array with $N\gg 1$ the purity is always $P_{r}\approx 0$ . The consequence is that Rydberg mixed excitation(s) almost do not evolve back to $\left|g\right\rangle$ under $\mathrm{L}_{g}$ . Thus, in our effective model we make the approximation that ${\rho_{\rm mix}}$ does not change under $\mathrm{L}_{g}$ driving.

Similarly to the non-Hermitian evolution case, the action of quantum jump $J$ on ${\rho_{\rm mix}}$ can also be simplified by changing the order of $J$ and the projection $\mathcal{P}_{\rm diag}$ , since $J$ is also a collection of single atom jumps

J{\rho_{{\rm{mix}}}}={\cal N}_{\left\{\alpha\right\},\left\{{\alpha^{\prime}}\right\}}{{\cal P}_{{\rm{diag}}}}\left[{J\prod\limits_{k=1}^{{n_{e}}}{a_{{M_{{\alpha_{k}}}}}^{\dagger}\otimes\bar{a}_{{M_{\alpha^{\prime}{{}_{k}}}}}^{\dagger}}}\right]\left|{0\otimes\bar{0}}\right\rangle.

(71)

Thus, the evolution of ${\rho_{\rm mix}}$ can be obtained by a projection $\mathcal{P}_{\rm diag}$ on the Lindblad dynamics of collective excitations $\{{a_{M_{\alpha}}^{\dagger}}\}$ [cf. Eq. 59], as shown in Eq. 69 and Eq. 71. The Fock space of these collective excitations expands a Hilbert space ${{\cal H}_{\rm mix}}$ .

Together with the Lindblad evolution of collective excitations in $\mathcal{H}_{\rm col}$ , the Hilbert space of our whole effective model is ${{\cal H}_{\rm eff}}={{\cal P}_{\rm Ryd}}\left[{{\cal H}_{\rm col}}\otimes{{\cal H}_{\rm mix}}\right]$ . Here, ${{\cal P}_{\rm Ryd}}$ projects out the states with more than two Rydberg excitations. The unitary evolution preserves the population within ${{\cal H}_{\rm col}}$ and ${{\cal H}_{\rm mix}}$ , while the Rydberg decoherence processes couple these two spaces. This structure is illustrated in Fig. 9.

The system density matrix ${\rho_{\rm eff}}(t)\in{\mathcal{H}_{{\text{\rm eff}}}}$ evolves under the following Lindblad master equation

\begin{gathered}{{\dot{\rho}}_{\rm eff}}\left(t\right)=-i\left[{H_{\rm eff}\left(t\right),{\rho_{\rm eff}}\left(t\right)}\right]\hfill\\ +\sum\limits_{n}{\frac{{1}}{2}\left[{2{C_{n}}{\rho_{\rm eff}}\left(t\right)C_{n}^{\dagger}-{\rho_{\rm eff}}\left(t\right)C_{n}^{\dagger}{C_{n}}-C_{n}^{\dagger}{C_{n}}{\rho_{\rm eff}}\left(t\right)}\right]},\hfill\\ \end{gathered}

(72)

with Hamiltonian

{H_{{\rm{eff}}}}\left(t\right)={{\cal P}_{{\rm{Ryd}}}}\left[{H_{0}^{{\rm{eff}}}+H_{rg}^{{\rm{eff}}}\left(t\right)+H_{rq}^{{\rm{eff}}}\left(t\right)+H_{rl}^{{\rm{eff}}}\left(t\right)}\right].

(73)

Defining ${\omega_{{M_{\alpha}}}}={\omega_{\alpha}}$ with $\alpha\in\left\{{r,q,l}\right\}$ , we have

{H_{0}^{{\rm{eff}}}}=\sum\limits_{\alpha=r,l,q,{M_{r}},{M_{l}},{M_{q}}}{{\omega_{\alpha}}}a_{\alpha}^{\dagger}{a_{\alpha}}+\frac{U}{2}{\hat{n}_{r}}\cdot\left({{{\hat{n}}_{r}}-1}\right),

(74)

with ${\hat{n}_{r}}=a_{r}^{\dagger}{a_{r}}+a_{{M_{r}}}^{\dagger}{a_{{M_{r}}}}$ , and

{H^{{\rm{eff}}}_{rg}}\left(t\right)=\frac{{{\Omega_{rg}}\left(t\right)}}{2}\left({{e^{-i{\omega_{rg}}t}}a_{r}^{\dagger}+h(c).}\right),

(75)

	$\displaystyle{H^{{\rm{eff}}}_{r\alpha}}\left(t\right)$	$\displaystyle=\frac{{{\Omega_{r\alpha}}\left(t\right)}}{2}\left[{e^{-i{\omega_{r\alpha}}t}}({a_{r}}\otimes a_{\alpha}^{\dagger}\right.$		(76)
		$\displaystyle\left.+{a_{{M_{r}}}}\otimes a_{{M_{\alpha}}}^{\dagger})+{\rm H.c.}\right],\quad\alpha=l,q.$		(76)

Following the result of the analytical solution of Eq. 61, we can write down the jump operators in Eq. 72 for the Rydberg decay process as

\begin{array}[]{*{20}{l}}{{C_{R\to G}}=\sqrt{{\Gamma_{r}}}{a_{r}},\qquad{C_{{M_{r}}\to G}}=\sqrt{{\Gamma_{r}}}{a_{{M_{r}}}},}\\ {C_{R\to{M_{\alpha}}}}=\sqrt{{\Gamma_{r}}}{a_{r}}\sigma_{{M_{\alpha}}}^{\dagger},\qquad{C_{{M_{r}}\to{M_{\alpha}}}}=\sqrt{{\Gamma_{r}}}{a_{{M_{r}}}}\sigma_{{M_{\alpha}}}^{\dagger},\\ \mathrm{with}\quad\alpha=l,q,\end{array}

(77)

where the $\sigma_{{M_{\alpha}}}^{\dagger}$ for $\alpha=r,l,q$ have same structure as oscillator creation operators, but without the $\sqrt{n}$ dependence. The jump operators for the Rydberg dephasing process are

{C_{R\to{M_{r}}}}=\sqrt{{\Gamma_{\phi}}}{a_{r}}\sigma_{{M_{r}}}^{\dagger},\qquad{C_{{M_{r}}\to{M_{r}}}}=\sqrt{{\Gamma_{\phi}}}\sigma_{\phi},

(78)

where ${\sigma_{\phi}}={\rm diag}\left({0,1,\sqrt{2}}\right)$ to capture the evolution in Eq. 64.

Benchmark with exact simulation

We provide a benchmark of our effective modeling by comparing its prediction with that obtained by numerically solving the master equation Eq. 23. In order to numerically reach larger system size for solving Eq. 23, we consider a simplified atomic-level structure that only involves states $\left|g\right\rangle,\left|q\right\rangle,\left|r\right\rangle$ as shown in Fig. 10(a), and include all errors discussed in this section. The benchmark on this simplified level structure is enough to demonstrate the validity of our effective modeling since the state $\left|q\right\rangle$ and $\left|l\right\rangle$ experience the same errors.

We test our effective model on a cyclic Raman transition process between the state $\left|0\right\rangle$ and $\left|{{1_{q}}}\right\rangle$ , which the evolution starts from state $\left|0\right\rangle$ and going through following process

\left|0\right\rangle\to\left|{{1_{r}}}\right\rangle\to\left|{{1_{q}}}\right\rangle\to\left|{{1_{r}}}\right\rangle\to\left|0\right\rangle\to...

(79)

back and forth. Such a population transfer is very common in our protocol. We count each transfer process between $\left|0\right\rangle$ and $\left|{{1_{q}}}\right\rangle$ as one time of Raman transfer. The errors appear when the Rydberg state $\left|r\right\rangle$ gets populated during each Raman transfer, thus this cyclic process allows us to see how the errors accumulate in the long-time evolution.

Specifically, for this cyclic Raman transfer process we alternatively apply $\pi$ pulses with resonant laser coupling ${{\rm{L}}_{\rm{g}}}$ and ${{\rm{L}}_{\rm{q}}}$ of constant Rabi frequencies ${\Omega_{rg}}$ and ${\Omega_{rq}}$ , with coupling time ${t_{rg}}=\pi/{\Omega_{rg}}$ and ${t_{rq}}=\pi/{\Omega_{rq}}$ . For simplicity, here we consider the beam profile of ${{\rm{L}}_{\rm{g}}}$ as a plane wave, such that ${u_{i}}=1/\sqrt{N}$ in Eq. 8. In the numerical calculation, we truncate the many-body Hilbert space ${{\cal H}}$ up to two excitations, which is enough to simulate the case of the photonic MPS generation with $D=2$ and $d=2$ . In this case, we have

	$\displaystyle{\cal H}={\rm{span}}({\left\|g\right\rangle^{\otimes N}},{{\left\{{\left\|{{q_{i}}}\right\rangle,\left\|{{r_{i}}}\right\rangle}\right\}}_{i=1\sim N}},$		(80)
	$\displaystyle{{\left\{{\left\|{{r_{i}}{q_{j}}}\right\rangle,\left\|{{q_{i}}{r_{j}}}\right\rangle,\left\|{{r_{i}}{r_{j}}}\right\rangle,\left\|{{q_{i}}{q_{j}}}\right\rangle}\right\}}_{i<j=1\sim N}}),$		(80)

with $\dim{{\cal H}}=2{N^{2}}+1$ .

The basis of ${{\cal H}_{\rm eff}}$ for the simplified level structure before applying the projection ${{\cal P}_{{\rm{diag}}}}$ is:

\begin{array}[]{l}\left|{{s_{r}}{m_{q}},{{\tilde{s}}_{{M_{r}}}}{{\tilde{m}}_{{M_{q}}}}}\right\rangle=\\ \sqrt{\dfrac{{{{\cal N}_{\left\{\alpha\right\},\left\{\alpha\right\}}}}}{{{s_{r}}!{m_{q}}!}}}{\left({a_{r}^{\dagger}}\right)^{{s_{r}}}}{\left({a_{q}^{\dagger}}\right)^{{m_{q}}}}{\left({a_{{M_{r}}}^{\dagger}}\right)^{{{\tilde{s}}_{{M_{r}}}}}}{\left({a_{{M_{q}}}^{\dagger}}\right)^{{{\tilde{m}}_{{M_{q}}}}}}\left|0\right\rangle.\end{array}

(81)

In Fig. 10(b), we pictorially demonstrate the coupling between a part of the basis in ${\cal H}_{\rm eff}$ during this cyclic Raman transfer process.

For demonstration, in Fig. 11(a) we compare the system dynamics obtained by solving the effective master equation [Eq. 72] with the exact solution obtained by solving the exact master equation [Eq. 23] for the following choice of parameters:

\begin{array}[]{l}N=20,\qquad{\Gamma_{r}}/{\Omega_{\max}}=0.016,\\ {\Gamma_{\phi}}/{\Omega_{\max}}=0.016,\qquad U/{\Omega_{\max}}=5.\end{array}

(82)

We see a good agreement between the exact soluation and prediction of the effective model. To quantify the validity of the effective model, we compare these two dynamics on the full density matrix level. The exact master equation simulation provides us with the time-dependent density matrix ${\rho}\left(t\right)\in{\mathcal{H}}$ , while the effective model provides another time-dependent density matrix ${\rho_{\rm eff}}\left(t\right)\in{\mathcal{H}_{{\text{eff}}}}$ . To directly compare ${\rho_{\rm eff}}\left(t\right)$ and ${\rho}\left(t\right)$ , we reconstruct a density matrix in ${{\cal H}}$ from ${\rho_{\rm eff}}\left(t\right)$ by mapping the corresponding basis from ${{\cal{H}}_{\rm eff}}$ to ${{\cal H}}$ . For example

\begin{array}[]{l}\left|{{1_{r}}}\right\rangle\left\langle{{1_{r}}}\right|\to\sum\limits_{i,j=1}^{N}{{u_{i}}u_{j}^{*}{e^{i{{\bf{k}}_{rg}}\cdot\left({{{\bf{r}}_{i}}-{{\bf{r}}_{j}}}\right)}}}\left|{{r_{i}}}\right\rangle\left\langle{{r_{j}}}\right|,\\ \mathrm{and}\quad\left|{{1_{{M_{r}}}}}\right\rangle\left\langle{{1_{{M_{r}}}}}\right|\to\sum\limits_{i,j=1}^{N}{{{\left|{{u_{i}}}\right|}^{2}}}\left|{{r_{i}}}\right\rangle\left\langle{{r_{i}}}\right|.\end{array}

(83)

From this mapping we obtain a density matrix $\rho_{{\rm{eff}}}^{{\rm{mb}}}(t)\in{\cal H}$ predicted by our effective model. Then we can compute the infidelity $1-{F_{{\rm{eff}}}}\left(t\right)$ between $\rho_{{\rm{eff}}}^{{\rm{mb}}}(t)\in{\cal H}$ and ${\rho}(t)$ , where

{F_{{\rm{eff}}}}\left(t\right)={\left({{\rm{Tr}}\sqrt{\sqrt{\rho\left(t\right)}\rho_{{\rm{eff}}}^{{\rm{mb}}}\left(t\right)\sqrt{\rho\left(t\right)}}}\right)^{2}}.

(84)

The evolution of the infidelity for the dynamics in Fig. 11(a) is shown in Fig. 11(b). During the entire evolution that contains 8 Raman transfers, the infidelity stays very low. We further obtain the infidelity after 10 Raman transfers for various parameters, which is shown in Fig. 11(c-d). In all the regimes that our protocol runs (see the parameter range in Fig. 3), we see good agreement between the exact simulation and the prediction from our effective model. Thus this effective modeling reliably captures the dynamics of the atomic array over a long time.

Appendix E Multi-photon emission and the construction of the process map $W_{P}$ for photon emission process

With a finite retrieval efficiency $p_{\rm em}$ , the photon emission process of $\left|{{1_{l}}}\right\rangle$ is

\left|{{1_{l}}}\right\rangle\xrightarrow{{{M_{P}}}}\left|0\right\rangle\left({\sqrt{{p_{{\rm{em}}}}}{{\left|1\right\rangle}_{{\rm{\rm ph}}}}\left|0\right\rangle_{\epsilon_{{\text{\rm ph}}}}+\sqrt{1-{p_{{\rm{em}}}}}{{\left|0\right\rangle}_{{\rm{\rm ph}}}}\left|1\right\rangle_{\epsilon_{{\text{\rm ph}}}}}\right),

(85)

Where ${\epsilon_{{\text{\rm ph}}}}$ is a photonic environmental mode used to capture the erroneous jump. Under the Holstein-Primakoff approximation [Eq. 11], the multiphoton emission property for a state $\left|{{m_{q}}{n_{l}}}\right\rangle$ can be obtained directly from the single-photon emission property [49], which gives

\begin{gathered}\left|{{m_{q}}{n_{l}}}\right\rangle\xrightarrow{{{M_{P}}}}\hfill\\ \sum\limits_{j=0}^{{n_{l}}}{\sqrt{{\binom{n_{l}}{j}}p_{{\text{em}}}^{j}{{\left({1-{p_{{\text{em}}}}}\right)}^{{n_{l}}-j}}}}\left|{{m_{q}}}\right\rangle{\left|j\right\rangle_{{\text{\rm ph}}}}\left|{{n_{l}}-j}\right\rangle_{\epsilon_{{\text{\rm ph}}}}.\hfill\\ \end{gathered}

(86)

Error due to the Holstein-Primakoff approximation

In real experiments there will be an additional photon retrieval error due to the deviation of the collective excited states in many-body Hilbert space $\cal H$ from that obtained under the Holstein-Primakoff approximation [Eq. 12]. For simplicity, let us consider the collective excitation profile being close to a plane wave as ${u_{i}}\approx 1/\sqrt{N}$ . In this case the exact form of ${\left|{{m_{q}}{n_{l}}}\right\rangle_{\rm ext}}$ in $\cal H$ and ${\left|{{m_{q}}{n_{l}}}\right\rangle}$ are

\begin{gathered}\begin{aligned} {\left|{{m_{q}}{n_{l}}}\right\rangle_{{\text{ext}}}}=&\sqrt{\frac{{{N^{{m_{q}}+{n_{l}}}}\left({N-{m_{q}}-{n_{l}}}\right)!}}{{{m_{q}}!{n_{l}}!N!}}}\hfill\\ &\cdot{\left({S_{q,{{\mathbf{k}}_{rg}}-{{\mathbf{k}}_{rq}}}^{\dagger}}\right)^{{m_{q}}}}{\left({S_{l,{{\mathbf{k}}_{rg}}-{{\mathbf{k}}_{rl}}}^{\dagger}}\right)^{{n_{l}}}}\left|0\right\rangle,\end{aligned}\hfill\\ \left|{{m_{q}}{n_{l}}}\right\rangle=\frac{1}{{\sqrt{{m_{q}}!{n_{l}}!}}}{\left({a_{q}^{\dagger}}\right)^{{m_{q}}}}{\left({a_{l}^{\dagger}}\right)^{{n_{l}}}}\left|0\right\rangle.\hfill\\ \end{gathered}

(87)

Approximately expanding the operators $S$ and $a$ in Eq. 87 both using spin operators [cf. Eq. 8], we can compute the overlap of these two states in the regime of ${m_{q}}+{n_{l}}\ll N$ as

\left\langle{{m_{q}}{n_{l}}}\right.{\left|{{m_{q}}{n_{l}}}\right\rangle_{{\text{ext}}}}\approx\sqrt{1-\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}.

(88)

To estimate the effect that this difference has on photon retrieval, we decompose $\left|{{m_{q}}{n_{l}}}\right\rangle$ as

	$\displaystyle{\left\|{{m_{q}}{n_{l}}}\right\rangle_{{\text{ext}}}}$	$\displaystyle\approx\sqrt{1-\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}\left\|{{m_{q}}{n_{l}}}\right\rangle\hfill$		(89)
		$\displaystyle+\sqrt{\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}{\left\|\phi\right\rangle_{{\text{err}}}},\hfill$		(89)

where ${\left|\phi\right\rangle_{\rm err}}$ is a state that is orthogonal to ${\left|{{m_{q}}{n_{l}}}\right\rangle}$ . We overestimate the photon retrieval error by assuming that ${\left|\phi\right\rangle_{\rm err}}$ does not generate a directional photon that can be retrieved. This lead to an state-dependent retrieval efficiency

{p_{{\rm{em}}}}\left({\left|{{m_{q}}{n_{l}}}\right\rangle_{\rm ext}}\right)=\left({1-\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}\right)p_{{\rm{em}}}^{{n_{l}}}.

(90)

For simplicity, we represent this state dependence by a renormalized single photon emission rate $p^{\prime}_{\rm em}$ for ${\left|{{m_{q}}{n_{l}}}\right\rangle}_{\rm ext}$ that scales as

{p^{\prime}_{{\rm{em}}}}\left({{{\left|{{m_{q}}{n_{l}}}\right\rangle}_{{\rm{ext}}}}}\right)={\left({1-\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}\right)^{1/{n_{l}}}}{p_{{\rm{em}}}}.

(91)

We further overestimate the multi-excitation retrieval error by determining ${p^{\prime}_{{\rm{em}}}}$ from largest excitation number, which relates to the bond dimension $D$ and physical dimension $d$ of the photonic state that we want to generate as follows:

\max{m_{q}}=D-1,\quad\max{n_{l}}=d-1.

(92)

Thus we get an lower bound of the retrieval efficiency [cf. Eq. 31]

{p^{\prime}_{{\rm{em}}}}\geq{\left({1-\frac{{\left({D+d-2}\right)\left({D+d-3}\right)}}{2N}}\right)^{1/\left({d-1}\right)}}{p_{{\rm{em}}}}.

(93)

For example in the cluster state generation, the state with maximal possible excitation number is $\left|{{1_{q}}{1_{l}}}\right\rangle$ , with $D=2$ and $d=2$ . Thus for the cluster state generation, we get a reduced retrieval efficiency ${p^{\prime}_{{\rm{em}}}}\geq\left({1-\frac{1}{N}}\right){p_{{\rm{em}}}}$ .

The construction of $W_{P}$

Now we can construct the the process map ${W_{P}}\mathrel{\mathop{\mathchar 58\relax}}{\mathcal{H}_{{\text{eff}}}}\to{\mathcal{H}_{{\text{eff}}}}\otimes{\mathcal{H}_{{\text{\rm ph}}}}$ for the photon emission process. In the effective modeling [Appendix D], the general basis of $\mathcal{H}_{\rm eff}$ is

\left|{{s_{r}}{m_{q}}{n_{l}},{{\tilde{s}}_{{M_{r}}}}{{\tilde{m}}_{{M_{q}}}}{{\tilde{n}}_{{M_{l}}}}}\right\rangle=\sqrt{\frac{{{{\cal N}_{\left\{\alpha\right\},\left\{\alpha\right\}}}}}{{{s_{r}}!{m_{q}}!{n_{l}}!}}}{\left({a_{r}^{\dagger}}\right)^{{s_{r}}}}{\left({a_{q}^{\dagger}}\right)^{{m_{q}}}}{\left({a_{l}^{\dagger}}\right)^{{n_{l}}}}{\left({a_{{M_{r}}}^{\dagger}}\right)^{{{\tilde{s}}_{{M_{r}}}}}}{\left({a_{{M_{q}}}^{\dagger}}\right)^{{{\tilde{m}}_{{M_{q}}}}}}{\left({a_{{M_{l}}}^{\dagger}}\right)^{{{\tilde{n}}_{{M_{l}}}}}}\left|0\right\rangle,

(94)

with ${{s_{r}}+{{\tilde{s}}_{{M_{r}}}}\leqslant 2}$ since we include at most two Rydberg excitations. By introducing three additional environmental modes ${\epsilon_{r}},{\epsilon_{{M_{r}}}},{\epsilon_{{M_{l}}}}$ , we can denote the emission process ${M_{P}}\mathrel{\mathop{\mathchar 58\relax}}{\mathcal{H}_{{\text{eff}}}}\to{\mathcal{H}_{{\text{eff}}}}\otimes{\mathcal{H}_{{\text{\rm ph}}}}$ as

{M_{P}}\mathrel{\mathop{\mathchar 58\relax}}\left|{{s_{r}}{m_{q}}{n_{l}},{{\tilde{s}}_{{M_{r}}}}{{\tilde{m}}_{{M_{q}}}}{{\tilde{n}}_{{M_{l}}}}}\right\rangle\to\sum\limits_{j=0}^{{n_{l}}}{\sqrt{{{{\binom{n_{l}}{j}}\left({{{p}^{\prime}_{{\text{em}}}}}\right)^{j}}}{{\left({1-{{p}^{\prime}_{{\text{em}}}}}\right)}^{{n_{l}}-j}}}}\left|{{m_{q}},{{\tilde{m}}_{{M_{q}}}}}\right\rangle{\left|j\right\rangle_{{\text{\rm ph}}}}\left|{{n_{l}}-j}\right\rangle_{\epsilon_{{\text{\rm ph}}}}{\left|{{s_{r}}}\right\rangle_{{\epsilon_{r}}}}{\left|{{{\tilde{s}}_{{M_{r}}}}}\right\rangle_{{\epsilon_{{M_{r}}}}}}{\left|{{{\tilde{n}}_{{M_{l}}}}}\right\rangle_{{\epsilon_{{M_{l}}}}}}.

(95)

Here we assumed that all excitations on $\left|l\right\rangle$ and $\left|r\right\rangle$ completely decays away after the photon emission process. From Eq. 95 we can construct $W_{P}$ as

{W_{P}}={\rm Tr}{{}_{{\epsilon_{{\text{\rm ph}}}},{\epsilon_{r}},{\epsilon_{{M_{r}}}},{\epsilon_{{M_{l}}}}}}\left[{{M_{P}}\otimes{{\bar{M}}_{P}}}\right].

(96)

Appendix F Computing photonic state fidelity using the matrix product density operator (MPDO) approach

In this section we provide the formalism to compute the fidelity ${{\cal F}_{\rm ph}}={\left\langle{{\psi_{\rm MPS}}}\right|{\rho_{\rm ph}}\left|{{\psi_{\rm MPS}}}\right\rangle}$ for generating $\left|{{\psi_{\rm MPS}}}\right\rangle$ [Eq. 4] with Rydberg-blockaded atomic array [cf. Section II], where ${\rho_{\rm ph}}$ is the eventually created photonic density matrix in our protocol including all errors. Let us consider the effective modeling of the atomic array [cf. Appendix D] with the density matrix ${\rho_{{\text{eff}}}}(t)\in{\mathcal{H}_{{\text{eff}}}}$ . The ${\rho_{\rm ph}}$ is determined by the dynamics of ${\rho_{\rm eff}}(t)$ , which follows Eq. 72. Denote ${N_{{\text{eff}}}}\equiv\dim{\mathcal{H}_{{\text{eff}}}}$ , Eq. 72 can be written as a vectorized form

\frac{{d\vec{\rho}_{\rm eff}\left(t\right)}}{{dt}}=\mathcal{L}\left(t\right)\vec{\rho}_{\rm eff}\left(t\right),\quad{\vec{\rho}_{\rm eff}}(t)=\sum\limits_{a,b=0}^{{N_{\rm eff}-1}}{{\rho_{ab}}}(t)\left|{a\otimes\bar{b}}\right\rangle.

(97)

Here $\mathcal{L}\left(t\right)$ is the Liouville operator, and ${\bar{b}}$ represent the complex conjugate of $b$ . The solution of Eq. 97 is

{\vec{\rho}_{\rm eff}}\left(T\right)=\mathcal{T}\left\{{{e^{\int_{0}^{T}{\mathcal{L}\left(t\right)dt}}}}\right\}{\vec{\rho}_{\rm eff}}\left(0\right)={W_{\cal L}}{\vec{\rho}_{\rm eff}}\left(0\right),

(98)

where the Liouvillian propagator ${W_{\cal L}}$ is of dimension $\dim\left({{W_{\cal L}}}\right)=N_{\rm eff}^{4}$ .

The process map ${W_{P}}$ [Eq. 96] for the photon emission process can be written as

{W_{P}}=\sum\limits_{i,j=0}^{d_{\rm max}-1}{\sum\limits_{a,b,c,s=0}^{{N_{\rm eff}}-1}{W_{P,abcs}^{ij}}}\left|{c,\bar{s},i,\bar{j}}\right\rangle\left|{a,\bar{b}}\right\rangle,

(99)

which maps the system density matrix with vectorized basis $|{a\otimes\bar{b}}\rangle$ to an source–photon joint density matrix with vectorized basis $|{c,\bar{s},i,\bar{j}}\rangle=|{c\otimes\bar{s}}\rangle\otimes{|{i\otimes\bar{j}}\rangle_{{\text{\rm ph}}}}$ . Here $d_{\rm max}-1>d-1$ denote the maximal photon number that can be generated in the emission process, although the $d^{\prime}-1>d-1$ levels are not used in the computation of ${\cal F}_{\rm ph}$ .

Thus in total, each round of photon generation consisting of an evolution under the master equation [Eq. 72] followed by the photon emission $W_{P}$ , results in a map

	$\displaystyle{{\vec{\rho}}_{\left[k\right]}}$	$\displaystyle=\sum\limits_{{i_{k}},{j_{k}}=0}^{d_{\rm max}-1}{N_{\left[k\right]}^{{i_{k}},{j_{k}}}{{\vec{\rho}}_{\left[{k-1}\right]}}},\hfill$		(100)
	$\displaystyle N_{\left[k\right]}^{{i_{k}},{j_{k}}}$	$\displaystyle=W_{P}^{{i_{k}}{j_{k}}}{W_{{{\cal L}_{\left[k\right]}}}}.\hfill$		(100)

Starting from an initial state $\left|{{\varphi_{I}}}\right\rangle\in{\mathcal{H}_{{\text{eff}}}}$ and successively operate the protocol for $n$ rounds, we arrive at a source–photon density matrix

{\vec{\rho}_{\left[n\right]}}=\sum\limits_{\left\{{{i_{k}},{j_{k}}}\right\}=0}^{d_{\rm max}-1}{N_{\left[n\right]}^{{i_{n}},{j_{n}}}...N_{\left[1\right]}^{{i_{1}},{j_{1}}}\left|{{\varphi_{I}},{i_{n}}...{i_{1}}}\right\rangle\left|{{\varphi_{I}},{j_{n}}...{j_{1}}}\right\rangle}.

(101)

The $n$ -photon density matrix ${\rho_{\rm ph}}$ is obtained by tracing out the source

{\vec{\rho}_{\mathrm{\rm ph}}}=\sum\limits_{\left\{{{i_{k}},{j_{k}}}\right\}=0}^{d_{\rm max}-1}{{\text{Tr}}\left[{N_{\left[n\right]}^{{i_{n}},{j_{n}}}...N_{\left[1\right]}^{{i_{1}},{j_{1}}}\tilde{B}}\right]\left|{{i_{n}},...,{i_{1}}}\right\rangle_{\rm ph}\left|{{j_{n}},...,{j_{1}}}\right\rangle_{\rm ph}},

(102)

with $\tilde{B}=\sum\nolimits_{a=0}^{N_{\rm eff}-1}{\left|{{\varphi_{I}}}\right\rangle\langle a|\otimes\left|{{\varphi_{I}}}\right\rangle\langle a|}$ . Defining $B=\left|{{\varphi_{I}}}\right\rangle\left\langle{{\varphi_{F}}}\right|$ , the fidelity ${{\cal F}_{\rm ph}}$ can then be efficiently evaluated as

	$\displaystyle{{\cal F}_{{\rm{\rm ph}}}}$	$\displaystyle=\sum\limits_{\left\{{{i_{k}},{j_{k}}}\right\}=0}^{d-1}{{\rm{Tr}}\left[{N_{\left[n\right]}^{{i_{n}},{j_{n}}}...N_{\left[1\right]}^{{i_{1}},{j_{1}}}\tilde{B}}\right]}$		(103)
		$\displaystyle\times{\rm{Tr}}\left[{\left({V_{\left[n\right]}^{{j_{n}}}\otimes\bar{V}_{\left[n\right]}^{{i_{n}}}}\right)...\left({V_{\left[1\right]}^{{j_{1}}}\otimes\bar{V}_{\left[1\right]}^{{i_{1}}}}\right)\left({B\otimes B}\right)}\right].$		(103)

Appendix G photon retrieval process

In this section, we present the details for calculating the photon retrieval efficiency ${p_{\rm em}}$ . Our calculations in this section follow Ref. [88] closely.

As explained in Section IV, we consider the photon emission from the singly-excited initial state [Eq. 26], and the decay dynamics can be described by an ansatz $\left|{\psi_{\rm em}\left(t\right)}\right\rangle=\left({{l_{j}}\left(t\right)\sigma_{lg}^{j}+{e_{j}}\left(t\right)\sigma_{eg}^{j}}\right)\left|0\right\rangle$ which evolves under a non-Hermitian Hamiltonian ${H_{\rm em}}(t)={H_{el}}(t)+{H_{{\rm{DDI}}}}$ [Eq. 27]. In our calculation we set the polarizations of all atoms along $x$ direction to maximize the photon emission speed along $z$ direction, that ${{\bf{d}}_{j}}={{\bf{x}}}\forall j$ [cf. Fig. 4(a)]. We get the equation of motion as

\begin{array}[]{l}{{\dot{e}}_{j}}={\rm{i}}\Delta{e_{j}}-{\rm{i}}\dfrac{{{\Omega_{{el}}}(t)}}{2}{l_{j}}+{\rm{i}}{\Gamma_{\rm em}}\sum\limits_{l}{{M_{jl}}}{e_{l}},\\ {{\dot{l}}_{j}}=-{\rm{i}}\dfrac{{{\Omega_{{el}}}(t)}}{2}{e_{j}},\end{array}

(104)

where ${\Gamma_{{\rm{em}}}}={\mu_{0}}\omega_{eg}^{3}d_{eg}^{2}/3\pi\hbar c$ is the spontaneous emission rate from $\left|e\right\rangle$ to $\left|g\right\rangle$ and ${M_{jl}}=3\pi k_{0}^{-1}{\bf{d}}_{j}^{*}\cdot{{\bf{G}}_{0}}\left({{{\bf{r}}_{j}},{{\bf{r}}_{l}},{\omega_{eg}}}\right)\cdot{{\bf{d}}_{l}}$ . Eq. 104 can be further simplified by noticing that ${M_{jl}}$ is a symmetric complex matrix, thus can be transpose-diagonalized [107], and we will obtain a set of eigenvalues $\left\{{\lambda_{\xi}}\right\}$ and corresponding eigenmodes $\left\{{{{\mathbf{v}}_{\xi}}}\right\}$ that satisfies ${\bf{v}}_{\xi}^{T}\cdot{{\bf{v}}_{{\xi^{\prime}}}}={\delta_{\xi{\xi^{\prime}}}}$ and $\sum\nolimits_{\xi}{{{\bf{v}}_{\xi}}}\otimes{\bf{v}}_{\xi}^{T}={\bf{I}}$ .

By solving atomic state evolution [Eq. 104] and relating it to the emitted photonic field via an input-output relation [108, 109, 88], one can compute the photon retrieval efficiency ${p_{\rm em}}$ as the probability to go to a defined detection mode ${{\bf E}_{\det}}\left(r\right)$ . One can define a detection mode operator [88]

{\hat{E}_{{\rm{det}}}}={\rm{i}}{d_{eg}}\sqrt{\frac{{{k_{0}}}}{{2{\epsilon_{0}}{F_{{\rm{det}}}}}}}\sum\limits_{j}{{\bf E}_{{\rm{det}}}^{*}}\left({{{\bf{r}}_{j}}}\right)\cdot{{\bf d}_{j}}\sigma_{ge}^{j},

(105)

where the factor ${F_{{\rm{det}}}}=\int_{z={\rm{const}}}{{{\rm{d}}^{2}}}{\bf{r}}{\bf{E}}_{{\rm{det}}}^{*}({\bf{r}})\cdot{{\bf{E}}_{{\rm{det}}}}({\bf{r}})$ normalizes ${\hat{E}_{{\rm{det}}}}$ such that $\left\langle{\hat{E}_{{\rm{det}}}^{\dagger}{{\hat{E}}_{{\rm{det}}}}}\right\rangle$ represents the photon number per unit time emitted into the detection mode. The retrieval efficiency ${p_{\rm em}}$ can be represented as

{p_{\rm em}}=\int_{0}^{\infty}{\rm{d}}t\left\langle{\hat{E}_{{\rm{det}}}^{\dagger}(t){{\hat{E}}_{{\rm{det}}}}(t)}\right\rangle.

(106)

By solving the dynamics of $\left|{\psi_{\rm em}\left(t\right)}\right\rangle$ , one can show [88]

{p_{{\rm{em}}}}=\dfrac{{{S_{{\lambda_{eg}}}}}}{{4{F_{{\rm{det}}}}}}\sum\limits_{j,l}^{N}{{u_{j}}{K_{jl}}u_{l}^{*}},

(107)

where ${S_{{\lambda_{eg}}}}=(3/2\pi)\lambda_{eg}^{2}$ is the resonant atomic optical cross-section. The matrix $K$ has a form

{K_{jl}}={\rm{i}}\sum\limits_{\xi,{\xi^{\prime}}=1}^{N}{{v_{\xi,j}}v_{{\xi^{\prime}},l}^{*}\frac{{E_{\xi}^{*}{E_{{\xi^{\prime}}}}}}{{{\lambda_{\xi}}-\lambda_{{\xi^{\prime}}}^{*}}}},

(108)

where ${E_{\xi}}=\sum\nolimits_{m}{{v_{\xi,j}}}{E_{j}}$ and ${E_{j}}={{\bf{E}}_{{\rm{det}}}}\left({{{\bf{r}}_{j}}}\right)\cdot{\bf{d}}_{j}^{*}$ .

As introduced in the main text, we consider two retrieval schemes: a uni-directional retrieval scheme where the photon goes to a single direction, and a two-directional retrieval scheme where the photon goes to two opposite directions (shown in Fig. 4). For the uni-directional retrieval scheme, we set the detection mode as a single vector Gaussian mode [92]

\begin{gathered}E_{{\text{det}}}^{x}({\mathbf{r}})={E_{0}}\int_{0}^{1}{\text{d}}bb{{\text{e}}^{-{b^{2}}k_{0}^{2}w_{0}^{2}/4}}{{\text{e}}^{{\text{i}}{k_{0}}z\sqrt{1-{b^{2}}}}}{J_{0}}\left({b{k_{0}}\rho}\right),\hfill\\ E_{{\rm{det}}}^{z}({\bf{r}})=-{\rm{i}}{E_{0}}\frac{x}{\rho}\int_{0}^{1}{\rm{d}}b\frac{{{b^{2}}{J_{1}}\left({b{k_{0}}\rho}\right)}}{{\sqrt{1-{b^{2}}}}}\cdot{{\rm{e}}^{-{b^{2}}k_{0}^{2}w_{0}^{2}/4+{\rm{i}}{k_{0}}z\sqrt{1-{b^{2}}}}},\hfill\\ \end{gathered}

(109)

where ${J_{0}}$ and ${J_{1}}$ are the Bessel functions. For the two-directional retrieval scheme, we choose a symmetric superposition of two vector Gaussian modes from two opposite directions as our detection mode.

The smallest photon retrieval error $\varepsilon_{\rm em}^{\rm opt}$ corresponds to the largest eigenvalue of the matrix $K$ in Eq. 107, with the corresponding eigenvector as its excitation profile, denoted as the optimal excitation profile. We also study the case of initial atomic excitation profile as a Gaussian mode ${l_{j}}\left(0\right)\propto E_{\det}^{x}\left({{{\bf{r}}_{j}}}\right)$ , which in principle can be easily created with the laser $\rm{L_{g}}$ and $\rm{L_{l}}$ , where $\rm{L_{g}}$ has the Gaussian profile and $\rm{L_{l}}$ has plane-wave profile. The retrieval error in this case is denoted as $\varepsilon_{\rm em}^{\rm Gauss}$ , as shown in Section IV.

Effect of array defects and atomic thermal motion

To model the effect of atomic defects, we generate random atomic configurations where a fraction (up to 20%) of positions in the array are unoccupied. For each fraction, we generate 100 realizations of the configuration. As shown in Ref. [88], one expects that the relative decrease of the retrieval efficiency should be proportional to the ratio of the detection mode hitting the empty sites. The results for the uni-directional (two-directional) retrieval is shown in Fig. 12(a) (Fig. 12(b)), where we separate a series of the intervals of $\sum\nolimits_{j\in\rm def}{{{\left|{{E_{j}}}\right|}^{2}}/\sum\nolimits_{l}{{{\left|{{E_{l}}}\right|}^{2}}}}$ , and average the photon retrieval errors within each interval. From these data we obtain the scaling of the modification of $p_{\rm em}$ due to atomic defects as in Eq. 30. we also see that $\alpha_{\rm def}$ decreases with the increase of array size, showing that an array with a larger size is more robust against this error. These results are consistent with Ref. [88].

To model the effect of atomic thermal motion, we add random spatial disorder $\delta{\mathbf{r}_{j}}=\left({{\delta_{x,j}},{\delta_{y,j}},{\delta_{z,j}}}\right)$ where each component is randomly distributed with a standard deviation $\sigma_{\rm th}$ . The error of the effect of photon retrieval efficiency are shown in Fig. 12(c) and Fig. 12(d), with each point obtained by averaging over 50 realizations of atomic positions with the same ${\sigma_{\rm th}}/d_{0}$ . In those figures, we see that the effect of this random disorder is proportional of ${\sigma_{\rm th}^{2}}/{d_{0}^{2}}$ , which is the same as predicted in Refs. [88, 95].

Appendix H Scaling of the achievable photon number with the bond and physical dimension of MPS

As pointed out in the introduction, MPS with moderate high bond and physical dimensions already find many applications, and can well capture the ground states of 1D local gapped Hamiltonians [34, 35, 36, 37]. Even for one-dimensional spin models at a critical point, a moderately high $D$ is enough to capture the ground state when this chain only have moderate system size, as the main deviation to the thermodynamic limit result comes from the finite size effect, instead of a finite-entanglement effect provided by the limited $D$ [110]. In this section, we provide a simple estimation of the scaling of the maximum achievable photon number as a function of the bond dimension and the physical dimension of MPS.

The photon number $N_{\rm ph}$ is determined by the error rate per photon as Eq. 33. To produce an MPS with bond dimension $D$ and physical dimension $d$ , we need to implement a unitary on the Hilbert space of dimension $2\cdot\dim{\cal H}_{src}^{\left({D,d}\right)}=2Dd$ . There is numerical evidence suggesting that [111] the time cost ${T_{D,d}}$ of implementing a general unitary in such a Hilbert space using the quantum optimal control approach scales as ${T_{D,d}}\sim{\left({Dd}\right)^{2}}$ . Thus, it takes more time to implement unitaries when we want to produce photonic MPS with higher bond and physical dimensions, which lead to stronger decoherence as the coefficients ${\beta_{U}},{\beta_{r}},{\beta_{\phi}}$ in Eq. 33 are proportional to ${T_{D,d}}$ . Also the retrieval efficiency is bounded by Eq. 31. In the regime of $D+d\ll N,{\beta_{0}}\approx 0$ and ${p_{{\rm{em}}}}\approx 1$ , the worst case estimation yields the following scaling of error

	$\displaystyle\xi_{{\rm{opt}}}^{D,d}$	$\displaystyle\approx{T_{D,d}}{\left({\dfrac{{27{\beta_{U}}}}{{4f_{{L_{v}},{L_{z}}}^{2}}}}\right)^{1/3}}{\left({\dfrac{{\Gamma d_{0}^{6}}}{{\left\|{{C_{6}}}\right\|}}}\right)^{2/3}}$		(110)
		$\displaystyle-{\beta_{{\rm{em}}}}\dfrac{{(D+d-2)(D+d-3)}}{{2\left({d-1}\right)N}}.$		(110)

Thus the dominant part lead to a qualitative scaling of ${N_{{\rm{\rm ph}}}}\sim{D^{-2}}{d^{-2}}$ .

Note that there are ways to reduce the error. As discussed in Section V, one can use higher Rydberg levels to improve the non-linearity and add a cavity to reduce the photon retrieval error.

References

O’Brien et al. [2009] J. L. O’Brien, A. Furusawa, and J. Vučković, Photonic quantum technologies, Nature Photonics 3, 687 (2009).
Gisin and Thew [2007] N. Gisin and R. Thew, Quantum communication, Nature photonics 1, 165 (2007).
Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
Degen et al. [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Reviews of Modern Physics 89, 1 (2017).
Burnham and Weinberg [1970] D. C. Burnham and D. L. Weinberg, Observation of simultaneity in parametric production of optical photon pairs, Physical Review Letters 25, 84 (1970).
González-Tudela et al. [2015] A. González-Tudela, V. Paulisch, D. E. Chang, H. J. Kimble, and J. I. Cirac, Deterministic Generation of Arbitrary Photonic States Assisted by Dissipation, Physical Review Letters 115, 1 (2015).
González-Tudela et al. [2017] A. González-Tudela, V. Paulisch, H. J. Kimble, and J. I. Cirac, Efficient Multiphoton Generation in Waveguide Quantum Electrodynamics, Phys. Rev. Lett. 118, 213601 (2017).
Gheri et al. [1998] K. M. Gheri, C. Saavedra, P. Törmä, J. I. Cirac, and P. Zoller, Entanglement engineering of one-photon wave packets using a single-atom source, Physical Review A 58, R2627 (1998).
Saavedra et al. [2000] C. Saavedra, K. M. Gheri, P. Törmä, J. I. Cirac, and P. Zoller, Controlled source of entangled photonic qubits, Physical Review A 61, 62311 (2000).
Schön et al. [2007] C. Schön, K. Hammerer, M. M. Wolf, J. I. Cirac, and E. Solano, Sequential generation of matrix-product states in cavity QED, Physical Review A 75, 1 (2007).
Lindner and Rudolph [2009] N. H. Lindner and T. Rudolph, Proposal for pulsed On-demand sources of photonic cluster state strings, Physical Review Letters 103, 1 (2009).
Nielsen and Mølmer [2010] A. E. B. Nielsen and K. Mølmer, Deterministic multimode photonic device for quantum-information processing, Physical Review A 81, 1 (2010).
[13] K. Tiurev, M. H. Appel, P. L. Mirambell, M. B. Lauritzen, A. Tiranov, P. Lodahl, and A. S. Sørensen, High-fidelity multi-photon-entangled cluster state with solid-state quantum emitters in photonic nanostructures, arXiv:2007.09295 .
Schön et al. [2005] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf, Sequential generation of entangled multiqubit states, Physical Review Letters 95, 1 (2005).
[15] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix Product State Representations, arXiv:0608197 [quant-ph] .
Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
Orús [2019] R. Orús, Tensor networks for complex quantum systems, Nature Reviews Physics 1, 538 (2019).
Lindner et al. [2016] N. H. Lindner, D. Cogan, O. Kenneth, E. R. Schmidgall, D. Gershoni, I. Schwartz, Y. Don, and L. Gantz, Deterministic generation of a cluster state of entangled photons, Science 354, 434 (2016).
Eichler et al. [2015] C. Eichler, J. Mlynek, J. Butscher, P. Kurpiers, K. Hammerer, T. J. Osborne, and A. Wallraff, Exploring Interacting Quantum Many-Body Systems by Experimentally Creating Continuous Matrix Product States in Superconducting Circuits, Physical Review X 5, 1 (2015).
Besse et al. [2020] J.-C. Besse, K. Reuer, M. C. Collodo, A. Wulff, L. Wernli, A. Copetudo, D. Malz, P. Magnard, A. Akin, M. Gabureac, G. J. Norris, J. I. Cirac, A. Wallraff, and C. Eichler, Realizing a deterministic source of multipartite-entangled photonic qubits, Nature Communications 11, 1 (2020).
[22] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions, arXiv:0407066 [cond-mat] .
Verstraete et al. [2008] F. Verstraete, V. Murg, and J. I. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Advances in Physics 57, 143 (2008).
Pichler et al. [2017] H. Pichler, S. Choi, P. Zoller, and M. D. Lukin, Universal photonic quantum computation via time-delayed feedback, Proceedings of the National Academy of Sciences 114, 11362 (2017).
Dhand et al. [2018] I. Dhand, M. Engelkemeier, L. Sansoni, S. Barkhofen, C. Silberhorn, and M. B. Plenio, Proposal for Quantum Simulation via All-Optically-Generated Tensor Network States, Physical Review Letters 120 (2018).
Xu and Fan [2018] S. Xu and S. Fan, Generate tensor network state by sequential single-photon scattering in waveguide QED systems, APL Photonics 3, 10.1063/1.5044248 (2018).
Economou et al. [2010] S. E. Economou, N. Lindner, and T. Rudolph, Optically generated 2-dimensional photonic cluster state from coupled quantum dots, Physical Review Letters 105, 1 (2010).
Gimeno-Segovia et al. [2019] M. Gimeno-Segovia, T. Rudolph, and S. E. Economou, Deterministic Generation of Large-Scale Entangled Photonic Cluster State from Interacting Solid State Emitters, Physical Review Letters 123, 70501 (2019).
Russo et al. [2019] A. Russo, E. Barnes, and S. E. Economou, Generation of arbitrary all-photonic graph states from quantum emitters, New Journal of Physics 21, 55002 (2019).
Bekenstein et al. [2020] R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin, and M. D. Lukin, Quantum metasurfaces with atom arrays, Nature Physics 16, 676 (2020).
Lubasch et al. [2018] M. Lubasch, A. A. Valido, J. J. Renema, W. S. Kolthammer, D. Jaksch, M. S. Kim, I. Walmsley, and R. García-Patrón, Tensor network states in time-bin quantum optics, Physical Review A 97, 1 (2018).
Greenberger et al. [1989] D. M. Greenberger, M. A. Horne, and A. Zeilinger, Going beyond Bell’s theorem, in Bell’s theorem, quantum theory and conceptions of the universe (Springer, 1989) pp. 69–72.
Briegel and Raussendorf [2001] H. J. Briegel and R. Raussendorf, Persistent entanglement in arrays of interacting particles, Physical Review Letters 86, 910 (2001).
Huang [a] Y. Huang, Computing energy density in one dimension, (a), arXiv:1505.00772 .
Dalzell and Brandão [2019] A. M. Dalzell and F. G. Brandão, Locally accurate MPS approximations for ground states of one-dimensional gapped local Hamiltonians, Quantum 3, 10.22331/q-2019-09-23-187 (2019).
Huang [b] Y. Huang, Approximating local properties by tensor network states with constant bond dimension, (b), arXiv:1903.10048 .
[37] N. Schuch and F. Verstraete, Matrix product state approximations for infinite systems, arXiv:1711.06559 .
Jarzyna and Demkowicz-Dobrzański [2013] M. Jarzyna and R. Demkowicz-Dobrzański, Matrix product states for quantum metrology, Physical Review Letters 110, 1 (2013).
Chabuda et al. [2020] K. Chabuda, J. Dziarmaga, T. J. Osborne, and R. Demkowicz-Dobrzański, Tensor-network approach for quantum metrology in many-body quantum systems, Nature Communications 11, 1 (2020).
Erhard et al. [2020] M. Erhard, M. Krenn, and A. Zeilinger, Advances in high-dimensional quantum entanglement, Nature Reviews Physics 2, 365 (2020).
Michael et al. [2016] M. H. Michael, M. Silveri, R. T. Brierley, V. V. Albert, J. Salmilehto, L. Jiang, and S. M. Girvin, New class of quantum error-correcting codes for a bosonic mode, Physical Review X 6, 1 (2016).
Contreras-Tejada et al. [2019] P. Contreras-Tejada, C. Palazuelos, and J. I. De Vicente, Resource Theory of Entanglement with a Unique Multipartite Maximally Entangled State, Physical Review Letters 122, 120503 (2019).
Zhou et al. [2003] D. L. Zhou, B. Zeng, Z. Xu, and C. P. Sun, Quantum computation based on d-level cluster state, Physical Review A 68, 11 (2003).
Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Valence-bond states for quantum computation, Physical Review A 70, 1 (2004).
Clark [2006] S. Clark, Valence bond solid formalism for d-level one-way quantum computation, Journal of Physics A: Mathematical and General 39, 2701 (2006).
Wang et al. [2017] D. S. Wang, D. T. Stephen, and R. Raussendorf, Qudit quantum computation on matrix product states with global symmetry, Physical Review A 95, 1 (2017).
Wei [2018] T. C. Wei, Quantum spin models for measurement-based quantum computation, Advances in Physics: X 3, 547 (2018).
Saffman and Walker [2002] M. Saffman and T. G. Walker, Creating single-atom and single-photon sources from entangled atomic ensembles, Physical Review A 66, 4 (2002).
Porras and Cirac [2008] D. Porras and J. I. Cirac, Collective generation of quantum states of light by entangled atoms, Physical Review A 78, 1 (2008).
Pedersen and Mølmer [2009] L. H. Pedersen and K. Mølmer, Few qubit atom-light interfaces with collective encoding, Physical Review A 79, 1 (2009).
Dudin and Kuzmich [2012] Y. O. Dudin and A. Kuzmich, Strongly interacting Rydberg excitations of a cold atomic gas, Science 336, 887 (2012).
Li et al. [2013] L. Li, Y. O. Dudin, and A. Kuzmich, Entanglement between light and an optical atomic excitation, Nature 498, 466 (2013).
Miroshnychenko et al. [2013] Y. Miroshnychenko, U. V. Poulsen, and K. Mølmer, Directional emission of single photons from small atomic samples, Physical Review A 87, 1 (2013).
Feng et al. [2014] W. Feng, Y. Li, and S. Y. Zhu, Effect of atomic distribution on cooperative spontaneous emission, Physical Review A 89, 1 (2014).
Petrosyan and Mølmer [2018] D. Petrosyan and K. Mølmer, Deterministic Free-Space Source of Single Photons Using Rydberg Atoms, Physical Review Letters 121, 1 (2018).
Jaksch et al. [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Fast quantum gates for neutral atoms, Physical Review Letters 85, 2208 (2000).
Lukin et al. [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Physical Review Letters 87, 37901 (2001).
D’Alessandro [2007] D. D’Alessandro, Introduction to quantum control and dynamics (Chapman and Hall/CRC, 2007).
Werschnik and Gross [2007] J. Werschnik and E. K. Gross, Quantum optimal control theory, Journal of Physics B: Atomic, Molecular and Optical Physics 40 (2007).
[60] J. Miguel-Ramiro, A. Pirker, and W. Dür, Genuine quantum networks: superposed tasks and addressing, arXiv:2005.00020 .
Giovannetti et al. [2008] V. Giovannetti, S. Lloyd, and L. MacCone, Quantum random access memory, Physical Review Letters 100, 1 (2008).
Holstein and Primakoff [1940] T. Holstein and H. Primakoff, Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet, Phys. Rev. 58, 1098 (1940).
Yuan and Lloyd [2007] H. Yuan and S. Lloyd, Controllability of the coupled spin-1 2 harmonic oscillator system, Physical Review A 75, 52331 (2007).
Mischuck and Mølmer [2013] B. Mischuck and K. Mølmer, Qudit quantum computation in the Jaynes-Cummings model, Physical Review A 87 (2013).
[65] T. Hofmann and M. Keyl, Controlling a d-level atom in a cavity, arXiv:1712.07613 .
Dudin et al. [2013] Y. O. Dudin, L. Li, and A. Kuzmich, Light storage on the time scale of a minute, Physical Review A 87, 1 (2013).
Yang et al. [2016] S. J. Yang, X. J. Wang, X. H. Bao, and J. W. Pan, An efficient quantum light-matter interface with sub-second lifetime, Nature Photonics 10, 381 (2016).
Law and Eberly [1996] C. K. Law and J. H. Eberly, Arbitrary control of a quantum electromagnetic field, Physical Review Letters 76, 1055 (1996).
Santos [2005] M. F. Santos, Universal and deterministic manipulation of the quantum state of harmonic oscillators: A route to unitary gates for Fock state qubits, Physical Review Letters 95, 1 (2005).
Strauch [2012] F. W. Strauch, All-resonant control of superconducting resonators, Physical Review Letters 109, 1 (2012).
Krastanov et al. [2015] S. Krastanov, V. V. Albert, C. Shen, C. L. Zou, R. W. Heeres, B. Vlastakis, R. J. Schoelkopf, and L. Jiang, Universal control of an oscillator with dispersive coupling to a qubit, Physical Review A 92 (2015).
[72] T. Fösel, S. Krastanov, F. Marquardt, and L. Jiang, Efficient cavity control with SNAP gates, arXiv:2004.14256 .
Schmidt-Kaler et al. [2004] F. Schmidt-Kaler, H. Häffner, S. Gulde, M. Riebe, G. Lancaster, J. Eschner, C. Becher, and R. Blatt, Quantized AC-Stark shifts and their use for multiparticle entanglement and quantum gates, Europhysics Letters 65, 587 (2004).
Lee et al. [2019] M. Lee, K. Friebe, D. A. Fioretto, K. Schüppert, F. R. Ong, D. Plankensteiner, V. Torggler, H. Ritsch, R. Blatt, and T. E. Northup, Ion-Based Quantum Sensor for Optical Cavity Photon Numbers, Physical Review Letters 122, 153603 (2019).
Caneva et al. [2009] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Optimal control at the quantum speed limit, Physical Review Letters 103, 1 (2009).
Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: From Heisenberg’s uncertainty principle to optimal quantum control, Journal of Physics A: Mathematical and Theoretical 50 (2017).
Levine et al. [2019] H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, Others, and M. D. Lukin, Parallel implementation of high-fidelity multiqubit gates with neutral atoms, Physical Review Letters 123, 170503 (2019).
Omran et al. [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Generation and manipulation of Schrödinger cat states in Rydberg atom arrays, Science 365, 570 (2019).
Machnes et al. [2018] S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm, Tunable, Flexible, and Efficient Optimization of Control Pulses for Practical Qubits, Physical Review Letters 120 (2018).
Palao and Kosloff [2002] J. P. Palao and R. Kosloff, Quantum Computing by an Optimal Control Algorithm for Unitary Transformations, Physical Review Letters 89, 3 (2002).
Khaneja et al. [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
Levine et al. [2018] H. Levine, A. Keesling, A. Omran, H. Bernien, S. Schwartz, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, High-Fidelity Control and Entanglement of Rydberg-Atom Qubits, Physical Review Letters 121, 1 (2018).
Saffman et al. [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with Rydberg atoms, Reviews of Modern Physics 82, 2313 (2010).
Verstraete et al. [2004] F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Matrix product density operators: simulation of finite-temperature and dissipative systems, Physical review letters 93, 207204 (2004).
Maller et al. [2015] K. M. Maller, M. T. Lichtman, T. Xia, Y. Sun, M. J. Piotrowicz, A. W. Carr, L. Isenhower, and M. Saffman, Rydberg-blockade controlled- not gate and entanglement in a two-dimensional array of neutral-atom qubits, Physical Review A 92, 10.1103/PhysRevA.92.022336 (2015).
Gorshkov et al. [2007] A. V. Gorshkov, A. André, M. Fleischhauer, A. S. Sørensen, and M. D. Lukin, Universal approach to optimal photon storage in atomic media, Physical Review Letters 98, 1 (2007).
Asenjo-Garcia et al. [2017a] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays, Phys. Rev. X 7, 31024 (2017a).
Manzoni et al. [2018] M. T. Manzoni, J. V. Porto, A. V. Gorshkov, D. E. Chang, M. Moreno-Cardoner, A. Asenjo-Garcia, J. V. Porto, A. V. Gorshkov, and D. E. Chang, Optimization of photon storage fidelity in ordered atomic arrays, New journal of physics 20, 83048 (2018).
Grankin et al. [2018] A. Grankin, P. O. Guimond, D. V. Vasilyev, B. Vermersch, and P. Zoller, Free-space photonic quantum link and chiral quantum optics, Physical Review A 98 (2018).
Lehmberg [1970] R. H. Lehmberg, Radiation from an $N$ -Atom System. I. General Formalism, Physical Review A 2, 883 (1970).
Novotny and Hecht [2012] L. Novotny and B. Hecht, Principles of nano-optics (Cambridge university press, 2012).
Chen et al. [2002] C. G. Chen, P. T. Konkola, J. Ferrera, R. K. Heilmann, and M. L. Schattenburg, Analyses of vector Gaussian beam propagation and the validity of paraxial and spherical approximations, Journal of the Optical Society of America A 19, 404 (2002).
Bao et al. [2012] X. H. Bao, A. Reingruber, P. Dietrich, J. Rui, A. Dück, T. Strassel, L. Li, N. L. Liu, B. Zhao, and J. W. Pan, Efficient and long-lived quantum memory with cold atoms inside a ring cavity, Nature Physics 8, 517 (2012).
Duan et al. [2002] L. M. Duan, J. I. Cirac, and P. Zoller, Three-dimensional theory for interaction between atomic ensembles and free-space light, Physical Review A 66, 1 (2002).
Shahmoon et al. [2017] E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, Cooperative Resonances in Light Scattering from Two-Dimensional Atomic Arrays, Physical Review Letters 118 (2017).
Walker and Saffman [2008] T. G. Walker and M. Saffman, Consequences of Zeeman degeneracy for the van der Waals blockade between Rydberg atoms, Physical Review A 77, 1 (2008).
Zeiher et al. [2015] J. Zeiher, P. Schauß, S. Hild, T. Macrì, I. Bloch, C. Gross, P. Schauß, S. Hild, T. Macrì, I. Bloch, and C. Gross, Microscopic characterization of scalable coherent Rydberg superatoms, Physical Review X 5, 31015 (2015).
Madjarov et al. [2020] I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, High-fidelity entanglement and detection of alkaline-earth Rydberg atoms, Nature Physics 16, 857 (2020).
Dunning et al. [2016] F. B. Dunning, T. C. Killian, S. Yoshida, and J. Burgdörfer, Recent advances in Rydberg physics using alkaline-earth atoms, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 0 (2016).
Schulte-Herbrüggen et al. [2011] T. Schulte-Herbrüggen, A. Spörl, N. Khaneja, and S. J. Glaser, Optimal control for generating quantum gates in open dissipative systems, Journal of Physics B: Atomic, Molecular and Optical Physics 44 (2011).
Abdelhafez et al. [2019] M. Abdelhafez, D. I. Schuster, and J. Koch, Gradient-based optimal control of open quantum systems using quantum trajectories and automatic differentiation, Physical Review A 99, 52327 (2019).
Wei [2019] Z. Y. Wei, Directional superradiant generation of strongly-entangled photonic states, Master thesis, Ludwig Maximilians University of Munich (2019).
Cox et al. [2019] K. C. Cox, D. H. Meyer, Z. A. Castillo, F. K. Fatemi, and P. D. Kunz, Spin-Wave Multiplexed Atom-Cavity Electrodynamics, 100, 1 (2019), arXiv:1907.04921 .
Fink et al. [2008] J. M. Fink, M. Göppl, M. Baur, R. Bianchetti, P. J. Leek, A. Blais, and A. Wallraff, Climbing the Jaynes-Cummings ladder and observing its √n nonlinearity in a cavity QED system, Nature 454, 315 (2008).
[105] R. Iten, O. Reardon-Smith, L. Mondada, E. Redmond, R. S. Kohli, and R. Colbeck, Introduction to UniversalQCompiler, arXiv:1904.01072 .
Plenio and Knight [1998] M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Reviews of Modern Physics 70, 101 (1998).
Asenjo-Garcia et al. [2017b] A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, Atom-light interactions in quasi-one-dimensional nanostructures: A Green’s-function perspective, Physical Review A 95, 1 (2017b).
Dung et al. [2002] H. T. Dung, L. Knöll, and D.-G. Welsch, Resonant dipole-dipole interaction in the presence of dispersing and absorbing surroundings, Phys. Rev. A 66, 63810 (2002).
Caneva et al. [2015] T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I. Cirac, and D. E. Chang, Quantum dynamics of propagating photons with strong interactions: a generalized input–output formalism, New Journal of Physics 17, 113001 (2015).
Pirvu et al. [2012] B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, Matrix product states for critical spin chains: Finite-size versus finite-entanglement scaling, Physical Review B 86, 75117 (2012).
Lee et al. [2018] J. Lee, C. Arenz, H. Rabitz, and B. Russell, Dependence of the quantum speed limit on system size and control complexity, New Journal of Physics 20, 1 (2018).

	$\displaystyle{\left\|{{m_{q}}{n_{l}}}\right\rangle_{{\text{ext}}}}$	$\displaystyle\approx\sqrt{1-\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}\left\|{{m_{q}}{n_{l}}}\right\rangle\hfill$		(89)
		$\displaystyle+\sqrt{\frac{{\left({{m_{q}}+{n_{l}}}\right)\left({{m_{q}}+{n_{l}}-1}\right)}}{{2N}}}{\left\|\phi\right\rangle_{{\text{err}}}},\hfill$		(89)