Boson sampling with ultracold atoms in a programmable optical lattice

Sampling from a quantum distribution can be exponentially hard for classical computers and yet could be performed efficiently by a noisy intermediate-scale quantum device. A prime example of a distribution that is hard to sample is given by the output states of a linear interferometer traversed by $N$ identical boson particles. Here, we propose a scheme to implement such a boson sampling machine with ultracold atoms in a polarization-synthesized optical lattice. We experimentally demonstrate the basic building block of such a machine by revealing the Hong-Ou-Mandel interference of two bosonic atoms in a four-mode interferometer. To estimate the sampling rate for large $N$, we develop a theoretical model based on a master equation that accounts for particle losses, but not include technical errors. Our results show that atomic samplers have the potential to achieve quantum advantage over today's best supercomputers with $N \gtrsim 40$.

A quantum sampling machine deals with the task of drawing from the probability distribution of the outcomes that are produced by measuring a quantum system in a highly entangled state.In essence, the idea is to use the randomness inherent to a measured quantum system to construct a hard-to-simulate sampling machine.Compared to other problems (e.g., decision problems), quantum sampling has the advantage that its computa-tional complexity can be ascertained for many quantum distributions [32] by relying only on a few widely held assumptions (e.g., no collapse of the polynomial hierarchy).Knowing the computational complexity of the specific problem allows us to gain important insights into the conditions (e.g., size of the Hilbert space) and class of quantum states [33][34][35][36] required to achieve a quantum advantage over classical machines [37,38].Quantum sampling also appeals for practical reasons because its computational hardness is generally robust to small experimental errors [39,40].Such a natural tolerance to errors makes quantum sampling a particularly suitable task to be performed with NISQ devices.Based on these motivations, several proposals have been put forward, where one draws samples from the state generated by a quantum circuit such as: constant-depth quantum circuits [41][42][43], instantaneous quantum polynomial-time circuits [44,45], random quantum circuits [39,[46][47][48][49], and linear quantum circuits of indistinguishable bosons [50], better known as boson samplers.
Boson sampling [51] refers to the problem of sampling from the probability distribution of the outcomes generated by N identical, noninteracting bosons that have travelled through an M -mode interferometer, with the initial and final N -particle states being of the form of Fock states.The state of each particle evolves according to the same M ×M unitary matrix U , mapping the initial into the final M modes of the linear interferometer.The probability of detecting a particular outcome comprising N bosons is proportional to the absolute square of the permanent of a N ×N submatrix of U [52,53].In spite of its compact analytical expression, the permanent (and likewise its absolute square) is very hard to compute, for it requires a time exponential in N [54].In fact, even its approximation to a multiplicative factor has been shown [55] to fall into the #P-hard complexity class [56].From a physics point of view, it is worth emphasizing that the hardness of this problem is due to the quantum statistics of indistinguishable bosons and not to interactions between particles [57,58].
A number of experiments have been reported demonstrating boson sampling in photonic quantum circuits [59][60][61][62][63][64][65][66][67][68][69][70], with the current record being N = 20 photons in M = 60 modes [71].Because of losses, however, it is hard to reach in the near future a much higher number of photons in a deterministic manner.This limitation along with the development [72] of more efficient classical algorithms for simulating boson sampling have prompted the study of variants of the problem that better cope with losses, such as lossy boson sampling [73,74], scattershot boson sampling [74][75][76][77] and Gaussian boson sampling [74,78].This latter in particular, which uses squeezed light instead of single photons, has demonstrated [79][80][81][82] a huge increase in the number of photons detected at the interferometer output, on the order of 100, leading to claims of a quantum advantage.
The quantum advantage of Gaussian boson sampling machines has been questioned, though, as it has been shown that classical sampling algorithms are able to efficiently draw samples from a sufficiently close distribution [83][84][85][86].For random quantum circuits [48,49], likewise, effective representations of the qubits' entangled state have been found [87,88] using tensor networks, which result in a tremendous speed-up of classical simulations, since only a tiny fraction of the Hilbert space is actually used when the gate fidelity is below a certain threshold.Such a race between quantum hardware and ever more efficient classical algorithms is indeed expected to continue in the coming years, promising new insights into what makes quantum systems advantageous from a computational complexity perspective.Remarkably, it was shown [36] based on fine-grained complexity arguments that a boson sampling quantum machine of the original type [50] with N = 100 and M = 500 achieves quantum advantage with respect to any (i.e., known and unknown) classical simulation algorithms.These numbers are large, yet not beyond the reach of scalable NISQ devices such as ultracold atoms in optical lattices.
In this article, we propose to use ultracold atoms in state-dependent optical lattices as a scalable architecture for boson sampling with hundreds of bosons.We also report on the experimental realization of the basic building block of the proposed boson-sampling machine, demonstrating the Hong-Ou-Mandel interference between two atoms trapped in state-dependent optical lattices.In our scheme, atoms cooled into their motional ground state play the role of identical bosons, while the lattice site as well as two internal atomic states serve as the bosonic modes.Distant modes, associated with different lattice sites, are brought together by state-dependent shift operations, which are realized with polarizationsynthesized optical lattices [89,90].Modes that are spatially overlapped are coupled in pairs, by a combination of microwave and site-resolved optical pulses, realizing the analog of phase-programmable photonic quantum circuits [91,92].
It should be mentioned that based on a similar motivation to establish a quantum advantage, other NISQ proposals alternative to photonic boson samplers have been put forward in the past years relying on trapped ions [93][94][95], superconducting circuits [96], and neutral atoms with microwave assisted tunneling [97].Furthermore, after the submission of this article, the first demonstration of boson sampling using ultracold atoms has been presented [98], involving 180 atoms spread across 1000 sites in a tunnel-coupled optical lattice.This recent experiment beautifully demonstrates the potential of ultracold atoms to realize a large-scale boson sampling device, once arbitrarily programmable quantum circuits can be created.Such arbitrary quantum circuits could be produced either using programmable tweezer arrays [99] or with state-dependent potentials, as discussed in this article.

II. BOSON SAMPLING WITH ATOMS IN OPTICAL LATTICES
A boson sampling quantum machine is in essence an M -port quantum circuit traversed simultaneously by N identical bosons that do not interact with each other.As there are no interactions between the particles, such a quantum circuit behaves as a linear interferometer, mapping each input mode into a superposition of the output modes, Here, a † i is the operator creating a boson in the i-th mode, and U ji is the matrix element of a unitary transformation U , which is randomly chosen from the uniform distribution (i.e., Haar measure) over all M × M unitaries.The randomness of U ensures that no particular feature can be exploited to efficiently simulate the boson sampler machine with a classical computer.
By detecting the occupation of the output modes, the machine thus directly samples from the probability distribution Here, n i denotes the number of bosons in the i-th output mode, and |ψ 0 ⟩ represents the initial state with N identical bosons, each occupying a particular input mode.According to best-known algorithms [72], sampling from P cannot be performed efficiently with classical computers, as it is bound to computing the permanent of N × N matrices, which requires a computation time of order O(N 2 2 N ) [100].
Importantly, to be hard to simulate by a classical computer, a boson sampler must have a number of modes much larger than that of particles, M ≫ N ≫ 1 [50].The gold standard satisfying this condition is given by the scaling law M = N 2 because it ensures that detecting two or more particles in any of the output modes has a small probability [101].In fact, only when the output modes are singly occupied, i.e., for the so-called collision-free outcomes, does the conjectured hardness of boson sampling hold [50].We therefore assume such a quadratic scaling in this paper.It should, however, be emphasized that this scaling has so far only been experimentally realized with a relatively small number of particles, N < 10, by photonics devices [71].
In the remainder of this section, we develop a concept how to implement such a boson sampling quantum machine with ultracold atoms in state-dependent optical lattices.We start with the key idea of how to construct arbitrary quantum circuits, and then discuss initialization and detection.
A. Arbitrary quantum circuits using polarization-synthesized optical lattices Such state-dependent shift operations can be performed using polarization-synthesized optical lattices [89].This piece of technology relies on the synthesis of polarization states of light to create two movable, fully independent periodic potentials, which selectively trap atoms in either the |+⟩ or |−⟩ internal state.Here, V 0 ± represents the trap depth and x ± (t) the position of the respective periodic potentials, x is the coordinate along the lattice direction, and λ L the wavelength of the lattice laser.The underlying concept of state-dependent optical potentials is suited to atomic species such as Rb and Cs [102,103].In this work, we will consider specifically the case of 133 Cs, where the internal modes are the hyperfine states |+⟩ = |F = 4, m F = 3⟩ and |−⟩ = |F = 3, m F = 3⟩, and λ L is set to the value of 870 nm, for which the trapping potential of right and left circularly polarized light selectively trap the two internal states.
The lattice potentials must be chosen sufficiently deep, with V 0 ± being of order of a few hundred recoil energies, to prevent atoms from tunneling to the neighboring sites.In such a deep-lattice regime, one can shift the atoms to the adjacent sites in a state dependent manner by simply varying the relative position, x + (t) − x − (t), as a function of time t.We have experimentally demonstrated [104] that repeated state-dependent shift operations preserve the coherence between the two internal states.Furthermore, we have shown in a recent work [105] that shifting the atoms by one lattice site can be rapid, with the minimum duration being bounded by the trap period at around 10 µs.
Crucially, a quantum circuit such as the one in Fig. 1, where the modes are locally coupled in alternating pairs, allow one to realize any arbitrary M ×M unitary transformation U of the input into the output modes [106,107].For a generic matrix U , a minimum number M (M −1)/2 of local operations T is required, arranged in a circuit of M -step depth [107].Such an operation T defines the basic unit of the quantum circuit, coupling together the modes |±⟩ associated with a given lattice site, where ϕ and θ are parameters depending on the particular site and time step.A programmable quantum device of this kind is said to be completely controllable [108].
We propose to implement T through composite pulses, where two types of elementary operations are stacked together: local phase imprints A(φ) and global Hadamard pulses H.In fact, for any T , one can find suitable angles θ and ϕ yielding the following decomposition: where A(φ) = exp[−iσ z φ/2] imprints onto the two modes |±⟩ a relative phase depending on the lattice site, whereas H = exp(−iσ x π/4) acts on all sites identically, realizing the equivalent of a beam splitter (here, σ i represent the Pauli matrices).Note that the local common-mode phase shift by ϕ/2, which appears in Eq. ( 4), can be avoided by conveniently adapting the algorithm by Clements et al. [107] to use T ′ instead of T as the basic unit.
The decomposition of T in Eq. ( 4) reveals a direct analogy to phase-programmable photonic circuits [91,92].Their structure reveals however an important difference: for ultracold atoms, a single spatial dimension suffices to wire the circuital modes, whereas at least two dimensions are necessary for photonic devices [70].The advantage of ultracold atoms simply arises from the fact that massive particles can be held in a specific location by a trapping potential.
The global Hadamard gates are readily implemented by means of π/2 microwave pulses, which act homogeneously on all lattice sites and require a time of order of 1 µs.The local operator A(φ) can be realized by exploiting the differential Stark shift that is produced by an array of laser beams focused on the target lattice sites through a high-numerical-aperture objective lens [109][110][111].Exploiting the vector polarizability of alkali atoms [102], one can imprint a purely differential phase shift φ onto the atoms by means of a circularly polarized light field.For Cs atoms, this condition is fulfilled when the wavelength of the addressing light field is chosen at λ A = 880 nm.Such local pulses also require that the addressing beam has a nonzero component along the quantization axis.This additional condition can be met by tilting the quantization axis with respect to the lattice direction (see Appendix A).
The differential phase shift φ is directly controlled by the product of the laser intensity and pulse duration.We estimate that the addressing pulses A(φ) can be realized in about 1 µs using approximately 10 µW of laser power per addressed lattice site.These pulses have a small impact on the coherence time of the atoms, because the probability that an atom scatters a photon off the addressing beam is approximately 4 × 10 −5 .
It is worth noting that the composite pulse scheme proposed to implement T offers an important advantage over other schemes that use local resonant pulses to couple the two hyperfine states, |±⟩.The reason for this is the difference in sensitivity to the crosstalk caused by the light field at the sites adjacent to the target lattice site.Differential Stark shift pulses, as in the proposed scheme, depend on the intensity of the addressing laser beams in their Hamiltonian, while resonant pulses directly depend on the respective electric fields.This different sensitivity implies that for the same intensity I leaking to the neighboring sites, crosstalk errors are smaller in the proposed scheme: the crosstalk infidelity is proportional to I 2 in the proposed scheme, in contrast to I in the resonant pulse schemes.

B. Initializing an array of identical atoms
Atom sorting techniques have been demonstrated [90,[112][113][114][115][116] where movable optical potentials are used to deterministically fill a predefined array of optical traps with one atom each.Recently, it has become possible to achieve densely-packed arrays with more than 1000 atoms [117,118] Once loaded into the desired sites of an optical lattice, atoms can be cooled to their motional ground state using sideband cooling techniques, making them indistinguishable in their mechanical degrees of freedom.Ground state probabilities above 90 % have already been achieved for atoms trapped in optical tweezers [115,[119][120][121], whereas higher values above 99 % are expected [122] for more tightly confined atoms in a threedimensional optical lattice.

C. Detection of individual atoms
The final state is measured by recording a fluorescence image of the atoms [111,123].Using a high-resolution objective lens, the positions of the individual atoms in the optical lattices can be reconstructed with a high fidelity, exceeding 99 %.Standard fluorescence-imaging techniques, however, only give information about the occupation of modes belonging to distinct lattice sites.To also resolve the occupation of modes associated with the same site, a state-sensitive detection scheme resolving the two hyperfine states, |±⟩, is needed.For this purpose, a long-distance state-dependent transport operation can be used to realize an optical Stern-Gerlach detection, mapping the two internal states to different lattice sites [124,125].Alternatively, one can use a magnetic Stern-Gerlach detection scheme in a multilayer optical lattice [126].
A particle number resolving detection is more demanding.Standard fluorescence imaging is not suitable because of light-induced collisions, which only allow the parity of the occupation number to be measured [127].One approach is to distribute the atoms to multiple sites prior to fluorescence imaging [128], similar to spatial multiplexing in photonic devices [63].An improved version of this approach consists in using a pinning lattice to detect the atoms [129,130].Alternatively, one can exploit interaction blockade to induce occupation-dependent tunneling to distinct sites of an optical lattice [126].

III. SCALING LAW OF AN ATOM-BOSON-SAMPLING MACHINE
To appreciate the quantum advantage of an atom boson sampler over classical machines, we study in this section the sampling rate, R, as a function of the number of bosons.

A. Ideal boson sampler
In the ideal case of no atom losses, the sampling rate is simply given by the repetition rate R rep at which N bosons are made to interfere with each other in a quantum circuit of M modes.Three terms determines this rate, where t init is the time for preparing the initial state of N atoms, t exec is the time for executing the quantum circuit, and t det is the time for detecting the atoms in the M output modes.
To estimate the execution time t exec , we assume that local Stark shift pulses can be carried out in parallel by optically addressing all sites simultaneously [16,131].The time for executing the quantum circuit is thus proportional to the total number of steps, which in turn is equal to the number of modes M (see Sec. II A).Considering that the number of modes is given by M = N 2 , as previously reasoned, we find that the execution time has a quadratic dependence on the number of bosons, t exec = N 2 t step , where t step is the time to perform the single step.
The initialization time t init is proportional to N for linear atom sorting [112,113], and to log(N ) for parallel atom sorting [90].For simplicity, we assume that this time is fixed at 500 ms, since the initialization of 100 atoms can be efficiently performed in less than this time [27].Likewise, we consider the detection time to be fixed, t det = 100 ms, since both positions and spin can be efficiently detected for all atoms in a single operation relying on fluorescence imaging [124,125].Scaling of boson sampling machines.Predicted sampling rate versus the number of bosons N is shown for atomic (blue), photonic (purple), and classical (grey) boson sampling machines.The atomic sampling is computed based on Eq. ( 6), where an extended model (see appendix B) is used for Pstep to account for both particle pairs and triplets.Our atomic sampling model disregards decoherence by technical noise sources; see text.The photonic and classical curves correspond to Eqs. (C1) and (C2), respectively.Under stateof-the-art conditions, atomic NISQ devices have the potential to overtake classical algorithms at N ≈ 37. Pioneering photonic experiments are marked with symbols [59], [60], and [61], while the best reported boson sampling experiment with [71].Note that the photonic experimental points refer to implementations of the original boson sampling problem, and thus do not include Gaussian boson sampling experiments, for the reasons presented in Sec.I. Our proofof-principle experiment, described in Sec.IV, is marked with symbol.The numerical results obtained with the master equation approach (see Appendix D) are indicated with symbols, with the vertical bars indicating the 1-σ statistical uncertainty.
Furthermore, if we post-select only those events with all atoms populating a different output mode (i.e., the socalled collision-free events), because these are the events hard to simulate with a classical machine, the sampling rate is reduced by a constant factor 1/e in the limiting case of large Thus, we conclude that under ideal conditions, an atom boson sampler can draw events from the boson distribution efficiently, since its computation time 1/R ideal scales with N 2 for sufficiently large N (i.e., polynomial time complexity).In contrast, classical computer simulations require an exponentially longer time to perform the same task, which scales with O(N 2 2 N ) [72] (see Sec. II).

B. NISQ boson sampler
In a realistic scenario typical of NISQ devices, the sampling rate is significantly degraded by state preparation errors, atom losses while executing the quantum circuit, and detection inefficiency.Along the lines of Refs.[68,72], we estimate the sampling rate as: where η det is the detection efficiency, η init the cooling efficiency into the motional ground state, and P surv is the probability that all N atoms survive.As reasoned in Sec.I, we only consider the case where no particle is lost.Equation ( 6) immediately reveals that the computation time 1/R NISQ scales exponentially with the number of particles.It is therefore important to carefully evaluate the expression in Eq. ( 6) to determine the conditions when an atom boson sampler offers a quantum advantage over classical computer simulations.
As described in Sec.II, cooling and detection of ultracold atoms can be done efficiently, with reported values of η init and η det above 90 % [115,[119][120][121] and 99 % [124], respectively.The survival probability P surv in Eq. ( 6) depends on two main loss mechanisms, which we discuss below.
The first mechanism responsible for the loss of atoms is collisions between one of the N trapped atoms and a particle from the background gas at room temperature, causing the atom involved in the collision to be ejected from the trap.For N atoms, the probability that no collision with the background gas occurs during a single step is given by the exponential formula P BG step = exp (−N t step /τ BG ), where τ BG is the background-gaslimited mean lifetime of a single atom.
The second mechanism leading to the loss of atoms is given by inelastic collisions of atoms occupying the same lattice site.Inelastic collisions cause the atoms to change their hyperfine state [132][133][134], and to acquire kinetic energy, thus leaving the motional ground state where they are initially prepared.Such inelastic collisions typically result in the loss of atoms from the trap because of the large energy separation between hyperfine states (several MHz between adjacent m F states, 9.2 GHz between different F states, expressed in frequency units).Threebody collisions are neglected because the probability of three (and more) particles occupying the same lattice site compared to that of pairs is negligible in the limit of a large number of modes (see Appendix B).To account for the two-atom lossy collisions, we introduce the survival probability of a pair of atoms located in the same lattice site, given by exp (−t/τ TB ), where τ TB is the mean lifetime limited by two-body collisions.Note that for simplicity we use the same constant τ TB without differentiating between the three possible spin configurations of the two bosons occupying the same site.To estimate the number of pairs of atoms that can collide onsite, we make the conservative assumption that all states of N bosons have equal probability of being occupied at every time step.This is a conservative assumption, which overestimates the probability of inelastic collisions at the initial steps since the atoms are first prepared in different sites.Based on these assumptions, the probability of finding k sites occupied by exactly a pair of atoms can be estimated as P k = e −3/2 (3/2) k /k! (see Appendix B).For example, P 0 is the probability of the collisionless subspace, where all atoms occupy distinct lattice sites.Thus, the survival probability per time step, limited by two-body lossy collisions, is obtained by the evaluating the following sum, where the expression on the right-hand side holds in the limit of large N (see Appendix B).For weak two-body losses, t ≪ τ TB , the survival probability P TB step decays as exp [−3t step /(2τ TB )], while in the limiting case of strong losses the survival probability approaches P 0 .
Combining the two loss mechanisms, the survival probability per step is simply given by P step = P BG step P TB step .Because the execution of the quantum circuit requires M = N 2 steps, the total survival probability of N atoms (where N ≫ 1) is thus A comparison of the two terms in Eq. ( 8) shows that under typical conditions t step ≪ τ TB ≪ τ BG two-body collisions are the dominant loss mechanism for N < N threshold = 3τ BG /(2τ TB ).As will be argued below, realistic experiments are expected to operate with a number of atoms below this threshold.
To evaluate R NISQ in Eq. ( 6), we consider two different scenarios, which are based on conservative and stateof-the-art assumptions, respectively.In the conservative scenario, we assume a step duration t step = 170 µs and a mean lifetime limited by two-body collisions τ TB = 40 ms, while in the state-of-the-art scenario, we consider t step = 33 µs and τ TB = 400 ms.For the initialization and detection times, we assume t init = 500 ms and t det = 100 ms, with efficiencies of η init = 0.90 and η det = 0.99 for the conservative scenario, and η init = η det = 0.99 for the state-of-the-art scenario.For the mean lifetime limited by background gas collisions, we take τ BG = 360 s in both scenarios.The threshold value N threshold is larger than 1000 atoms in both scenarios, implying that for a realistic number of atoms the dominant loss mechanism is inelastic two-body collisions rather than collisions with the background gas.
In Fig. 2, we show the sampling rate R NISQ as a function of the number of particles N , computed for an atom boson sampler with conservative and state-of-the-art assumptions (blue curves).To identify the regime of potential quantum advantage, we present in the same figure the sampling rate R classical of best algorithms [72] simulating a boson sampler using a standard laptop and the Tianhe-2 supercomputer (gray curves).For the sake of comparison, we also report in the figure the sampling rate expected for NISQ photonic devices (purple curves) for a conservative and state-of-the-art scenario; see Appendix C details.
To validate our model of the sampling rate R NISQ in Eq. ( 6), we have carried out exact numerical simulations based on a master equation approach (Appendix D) for N up to 4. The result of the numerical simulations (blue squares in Fig. 2) shows a very good agreement with the curves from the model.
In the scenario of N = 50, which is relevant for quantum advantage, we find that the sampling rate of the fastest atomic sampler proposed in this work (t step = 33 µs) is dominated by initialization and detection times.We expect a different situation in tunnel-coupled optical lattices [98], where a tunneling event (i.e., the equivalent of a time step) takes about t step ≈ 1 ms.For these samplers, the execution time t exec = N 2 t step is of order of a few seconds for N = 50 particles and will likely be the factor limiting the sampling rate.
We should emphasize that our model takes into account interactions between atoms only through their dissipative effects, which are accounted for in Eq. ( 6) in terms of atom losses.The model does not consider interaction-induced phase shifts.Modeling coherent interactions between atoms requires a different approach [135], and goes beyond the scope of this work.We expect that the interaction-induced phase shifts will make this boson sampling problem even harder to simulate on classical computers compared to its linear counterpart, as was originally proposed [55].
Finally, we would like to point out that, both in our theoretical model of Eq. ( 6) and in the numerical simulations, we have not considered technical errors affecting the operations during the quantum circuit.Future investigations on the scaling of atomic sampling devices will have to address not only atom losses, but also the loss of coherence.The specific technical errors will depend on the specific implementation, and a detailed analysis thereof goes beyond the scope of this work.We refer the reader to the appendix of Ref. [136] for a discussion of the main technical decoherence mechanisms, and to [89] for a characterization of the noise sources affecting polarization-synthesized optical lattices.The present work rather focuses on quantifying the effect of fundamental errors that exist across all neutral atom platforms.

IV. EXPERIMENTAL DEMONSTRATION OF THE HONG-OU-MANDEL INTERFERENCE
We have performed a proof-of-principle experiment with two atoms in a four-mode interferometer, which demonstrates the Hong-Ou-Mandel effect with atoms, as schematically shown in the inset of Fig. 1.Such an experiment establishes the basic building block of the envisaged boson sampling machine.
The Hong-Ou-Mandel effect with atoms has been previously demonstrated experimentally using movable optical tweezers [137], an optical lattice superimposed to a box potential [138], and a free-fall atom interferometer [139].Compared to these setups [140] and to related pro-posals based on microwave-induced tunneling [97], our setup is distinguished by the way modes are coupled, where the atoms are moved with state-dependent shift operations [89] instead of having them tunnel through an optical potential barrier.Our approach enables faster operations on the scale of few microseconds instead of milliseconds.
The setup used for our experimental demonstration is schematically depicted in Fig. 3(a).We start with a handful of Cs atoms, which are sparsely loaded in a onedimensional polarization-synthesized optical lattice [89].Using the atom sorting technique presented in Ref. [90], a pair of atoms is then selected and repositioned to a relative distance of twenty lattice sites with a success rate of about 99 %, mainly limited by an incorrect detection of the initial distance between the two atoms [123].
To make the two atoms identical, we cool them to the ground state of the lattice site potential in which they are respectively trapped.For this purpose, we use resolved sideband cooling, where the sideband transitions are driven by microwave radiation [141] for the direction along the lattice axis and by two Raman lasers [119] for the transverse directions.The optical lattice provides a tight confinement (ν ∥ ≈ 120 kHz sideband) along its longitudinal direction, while a hollow-tube potential collinear with the lattice axis also provides a tight confinement (ν ⊥ ≈ 20 kHz sideband) in the transverse directions.These trap frequencies are much larger than the recoil frequency (≈ 2 kHz), thus ensuring that the Lamb-Dicke condition necessary for ground state cooling is fulfilled.We alternate between microwave and Raman FIG. 3. Experimental setup probing the Hong-Ou-Mandel interference with two atoms.(a) Atoms are trapped in a polarization-synthesized optical lattice, formed by two optical standing waves of left (L) and right (R) circular polarization, which can be shifted by varying the optical phases ϕL and ϕR (adapted after [89]).Two Raman laser beams perform transverse sideband cooling.A spiral phase plate (SPP) generates a hollow-tube laser beam (intensity profile in the inset), enhancing the transverse confinement of the atoms.Cooling along the longitudinal and transverse directions is evidenced by the suppressed blue detuned sidebands in the (b) microwave and Raman (c) sideband spectra, respectively.
sideband cooling 3 times in order to cool the atoms in both the longitudinal and transverse directions.By allowing a slight ellipticity of the transverse potential, we lift the degeneracy of the transverse motional states, allowing Raman sideband cooling to be effective along both transverse directions.At the end of the cooling process, the atoms are polarized in state |+⟩ = |F = 4, m F = 4⟩ with a probability of 99 %.Note that in this section |+⟩ refers to a different Zeeman state than the one considered in the rest of this work (cf.Sec.II A).With the m F state chosen for the experimental demonstration, it must be taken into account that the state-dependent potentials depend on both left and right circular polarization components of the trapping light field [105].
Figures 3(b) and 3(c) report a typical microwave and Raman sideband spectrum recorded after the cooling procedure, demonstrating a pronounced suppression of the cooling (blue detuned) sideband with respect to the heating (red detuned) sideband.From the ratio of the sideband amplitudes [142], we derive a ground state probability of ≈ 95 % for the longitudinal direction and of ≈ 84 % for each transverse direction.Thus, the overall probability of occupying the motional ground state can be estimated as P 3D ≈ 95 % × 84 % × 84 % ≈ 67(9) %.
After sideband cooling, a magnetic field gradient along the lattice direction (11.6 G/cm) is ramped up in 10 ms and maintained until detecting the atoms by fluorescence imaging.The magnetic field gradient induces a position-dependent Zeeman shift (≈ 1.2 kHz per lattice site), which is used to selectively transfer one of the two atoms to the |−⟩ state.We perform such a selective spin flip by addressing the target atom with a microwave narrow-bandwidth π pulse (Gaussian shape, 7 kHz rms width).Spin-flip errors are clearly visible in the final fluorescence image, allowing them to be removed by postselection.The atom thus selected is then adiabatically transported in 1 ms to the site of the second atom by shifting the V − (x) lattice potential.With the two atoms occupying the same site, we apply a fast microwave π/2 pulse (square shape, 4.8 µs duration).This pulse acts much like the beam splitter of a Hong-Ou-Mandel optical interferometer, erasing the which-way information of the two impinging particles.Last, we shift the V + (x) lattice potential to map the internal states, |+⟩ and |−⟩, to two different locations 10 lattice sites apart, where the atoms are detected by position-resolved fluorescence imaging.
Two identical atoms are expected to bunch together in the same lattice site with unit probability because of the Hong-Ou-Mandel interference (bosonic quantum statistics).In practice, however, the two atoms can differ from each other because of their motional states.For two atoms in orthogonal states (i.e., fully distinguishable particles), the outcomes resemble those obtained from the toss of two independent coins, yielding a bunching probability of 1/2 (classical statistics).For partially distinguishable atoms like ours, the probability to bunch is determined by the so-called quantum purity of the state, γ ≈ P 2 3D , which represents the probability of the two atoms to be indistinguishable: Thus, the Hong-Ou-Mandel interference is established if we can show that the bunching probability of the two atoms fulfills P bunch > 1/2.
In the experiment, we distinguish three outcomes corresponding to the detection of zero, one, and two atoms in the final fluorescence image.Figure 4 shows the experimentally recorded probability for each of them.When both atoms are positively detected, the atoms are found at different locations, 10 sites apart (see fluorescence image in Fig. 4).Importantly, such an outcome can only occur when the two atoms have not bunched together.Its probability is thus directly related to P bunch through the expression: P 2 = S 2 (1 − P bunch ), where S is the single-atom survival probability.From the measurement of P 2 = 19(3) % and the independent characterization of the single-atom survival probability S = 84(1) % (see Appendix E), we therefore obtain P bunch = 73(4) %, which exhibits a 5-σ deviation from the reference value 1/2.
Such a value of P bunch establishes that the Hong-Ou-Mandel interference of the two atoms occurs with a probability γ = 2P bunch − 1 = 45(8) %, in good agreement with the value expected from the independent measurement of the ground state fraction, γ ≈ P 2 3D ≈ 45(12) %.We expect that γ can be significantly improved in the future with more efficient ground-state cooling of the atoms in a three-dimensional optical lattice.An analysis of all experimental outcomes, including those with zero and one atom detected, yields a value of P bunch = 73(6) % that is statistically consistent with the value we have derived from P 2 only (see Appendix E).
In the future, it will be interesting to test not only the degree of indistinguishability between atoms, but also their degree of entanglement.For that, one could study cross-correlations with a two-atom interferometer as presented in Refs.[143][144][145][146].

V. CONCLUSIONS
In this work, we have presented a scheme for the realization of programmable NISQ circuits with neutral atoms in state-dependent optical lattices.Quantum circuits based on the proposed scheme can be easily reprogrammed and scaled up to hundreds of modes.Both repogrammability and scalability are key to realize large random unitaries a prerequisite for any large-scale boson sampling machine.Furthermore, we have experimentally demonstrated the basic building block of an atom boson sampling machine by executing a quantum circuit with four modes and two indistinguishable atoms.We observed the atoms bunching in pairs, thus revealing their bosonic nature.The Hong-Ou-Mandel interference signal of the atom bunching is found to deviate from the outcome predicted for distinguishable (i.e., classical) particles by 5 σ.The degree of indistinguishability of the atoms is determined by the probability of occupying the motional ground state in the potential well of an optical lattice site.We independently measured the motional ground state occupancy of the atoms and showed that it is in good agreement with the ground state occupancy inferred from the observed bunching probability.
We have discussed in detail how to wire quantum circuits using a one-dimensional state-dependent optical lattice.Our analysis of NISQ devices has shown that controlling more than M = 500 lattice sites will be required to reach a quantum advantage over best supercomputers.Controlling such a large number of lattice sites may become difficult to realize in a one-dimensional geometry.However, our scheme can be readily extended to twodimensional state-dependent optical lattices [147], leading to a more compact and less resource-intensive platform.
For future studies, it will be interesting to investigate the role of controlled coherent collisions among atoms [148], which can be exploited to imprint collisional phases onto the quantum state when two or more particles meet at the same site [149].The inclusion of such nonlinearities is shown to augment the amount of correlations in the output distribution [135].For this reason, there is a potential that simulating a nonlinear boson sampler with classical computers will be an even harder task than the original linear problem.
Finally, we emphasize that beyond the boson sampling application, the proposed scheme can be used to implement reprogrammable parametrized quantum circuits a key component for quantum machine learning [28,29].
After the completion of our manuscript, the first demonstration of an atomic boson sampler has been put forward [98], involving 180 atoms spread across 1000 optical tweezer sites.

APPENDIX A: QUANTIZATION AXIS
The propagation direction of the addressing beam must be perpendicular to the optical lattice itself in order to allow it to be tightly focused on the individual sites of the lattice itself.The only way to satisfy this geometric condition and the condition stated in Sec.II A that the quantization axis must have a nonzero component along the addressing beam is to have the quantization axis form an angle α ̸ = 0 with the direction of the optical lattice.
However, tilting the quantization axis by α with respect to the optical lattice laser beams causes a wobbling of the optical lattice potential depth, V 0 ± , during the state-dependent shift operations.Because of this wobbling, the lattice depth is reduced by a factor | cos(α)| during transport.The extent of this effect is very small for small angles.For example, for α = 15 • the lattice depth is reduced by only ≈ 3 %.This wobbling can be easily taken into account and compensated for in the design of transport operations by optimal control [105].

APPENDIX B: PAIR DISTRIBUTION AND EXTENDED MODEL
In the limit of large N and after averaging over random quantum circuits, the state of the N bosons is described by a uniform statistical mixture ρ u , where all possible configurations of the N bosons over the M modes are equally weighted [101].Using this result, we derive a general expression P (k 2 , k 3 ) to account for the probability of finding exactly k 2 sites occupied by a pair and k 3 sites occupied by a triplet, whereas the other sites are singly occupied.
To calculate P (k 2 , k 3 ), we first determine the number of configurations containing k 3 sites occupied by a triplet of atoms.For that, we consider that there are M/2 k3 possible combinations to arrange k 3 triplets in M/2 sites, where we assume for simplicity that the overall number of modes, M , is even.Given the two spin states, for each site occupied by a triplet, we have 2 3 = 4 possible spin configurations, where the double brackets denote the multiset coefficient.Thus, the previous number of combinations must be multiplied by 4 k3 to obtain the total number of configurations for the k 3 triplets.Next, we consider that there are M/2−k3 k2 combinations to arrange k 2 pairs in the remaining M/2 − k3 sites.This number must be multiplied by 3 k2 to account for the 2 2 = 3 different spin configurations for each pair.Now, there are M/2−k2−k3 N −2k2−3k3 different combinations to place the re- < l a t e x i t s h a 1 _ b a s e 6 4 = " P 3 h z X l 0 X p x 3 5 2 P e u u J k M 0 f w B 8 7 n D z p 6 j s o = < / l a t e x i t > P (k 2 , 0) < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 F f M Z i Z y a p t R 5 j w 3 r C q C V k 9 4 0 2 g = " > A A A B 6 n i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h d 0 o 6 k U I e v E k E c 0 D k i X M T n q T I b O z y 8 y s E J Z 8 g h c P i n j 1 i 7 z 5 N 0 4 e B 0 0 s a C i q u u n R Z i h 2 i l 8 t b s x S y O U h g m q d c t z E + N n V B n O B I 7 y 7 V R j Q t m A 9 r B l q a Q R a j + b n D o i x 1 b p k j B W t q Q h E / X 3 R E Y j r Y d R Y D s j a v p 6 3 h u L / 3 m t 1 I S X f s Z l k h q U b L o o T A U x M R n / T b p c I T N i a A l l i t t b C e t T R Z m x 6 e R t C N 7 8 y 4 u k X i 5 5 5 6 X y / V m x c j 2 L I w e H c A Q n 4 M E F V O A W q l A D B j 1 4 h l d 4 c 4 T z 4 r w 7 H 9 P W J W c 2 c w B / 4 H z + A J 1 N j V 0 = < / l a t e x i t > N = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 f i 9  For increasing N , P (0, k3 > 0) tends to zero.
maining N − 2k 2 − 3k 3 particles in the rest of the sites.
For each singly occupied site, there are only two possible spin configurations.Finally, to obtain the probability P (k 2 , k 3 ), all the configurations should be divided by the overall number of possible bosonic configurations, which is given by the multiset coefficient M N .Hence, the probability of having exactly k 2 sites occupied by pairs and k 3 sites occupied by triplets is given by In the limit of large N , the expression in Eq. (B1) yields a vanishing probability to have k 3 > 0 triplets, lim N →∞ P (k 2 , k 3 > 0) = 0. Instead, when k 3 = 0, Eq. (B1) converges to a Poisson distribution for the occurrence of k 2 pairs, lim with average value λ = 3/(2c).Here, the constant factor c denotes the ratio c = M/N 2 , which in the main text has simply been assumed equal to 1. Figure 5 shows the distributions P (k 2 , 0) (left panel) and P (0, k 3 ) (right panel) for an increasing number of particles N = 3, 9 and 27, assuming c = 1.The graphs show that the probability associated with the collision-free subspace, P (0, 0), tends to the finite value exp[−3/(2c)].We also notice that for a relatively small number of particles N ≲ 10, the probability of finding a triplet, P (0, 1), has a nonnegligible value between 0.05 and 0.1.
It is straightforward to generalize Eq. (B1) to higher occupations (quartets, quintets, etc).However, the probability of having more than 3 particles in a site is negligible even for small number of particles and can thus be safely neglected.
The extended distribution P (k 2 , k 3 ) can be used to estimate the effect of two-body losses on the uniform state ρ u .Two-body losses are assumed to reduce the survival probability of a state containing k 2 pairs by a factor e −k2t/τTB .A particle triplet decays three times faster than a pair because there are 3 possible combinations for the colliding pairs.This can also be understood through the action of the particle decay operator V (t) (defined in Appendix D) on a state with k 3 particle triplets, where V (t)|k 3 ⟩⟨k 3 |V † (t) = exp (−3k 3 t/τ TB )|k 3 ⟩⟨k 3 |.The total survival probability is thus obtained evaluating the sum The sum in Eq. (B3) can be simplified in the limit of large N by ignoring particle triplets.Using the approximate expression in Eq. (B2) for P (k 2 , 0), we obtain B4) which corresponds to Eq. ( 7) in the main text when c = 1.The expression on the right-hand side of Eq. (B4) is obtained by computing the limit as N tends to infinity.
Note that in Fig. 2 and in Appendix D, the extended model of Eq. (B3) is used instead of Eq. ( 7) to describe the two-body losses, because it also applies in the limit of a small number of particles.

APPENDIX C: SCALING LAW FOR PHOTONIC AND CLASSICAL BOSON SAMPLERS
Along the lines of Ref. [68], the sampling rate of a photonic boson sampler is given by where R 0 /N is the rate in which N indistinguishable photons are created and η = η f η d c is the success probability of a single photon going through state preparation, circuit transmission, and successful detection.The success probability η is a product of a fixed probability η f , which does not scale with the number of modes in the circuit and depends only on the state preparation and detection, and the circuit transmission probability η M c , which accounts for the chance that a photon is absorbed at each beam splitter.For a square circuit [107] with M = N 2 , the sampling rate R ph is mainly determined by the term η N M c for sufficiently large N .Our conservative values (lower edge of the purple area in Fig. 2) are based on Ref. [71].The authors use a single-photon source working at R 0 = 76 MHz and report a boson sampling rate of R ph = 295 Hz for 5-photons in a 60×60 optical circuit.The transmission probability through the entire circuit is reported as 98.7% and correspondingly η c = 0.987 1/60 , which in turn allows us to extract the fixed preparation and detection probability of η f = 0.14.For the state-of-the-art estimate (upper edge of the purple area in Fig. 2), we assume an overall increased fixed preparation and detection probability of η f = 0.65 [69].
For a classical boson sampler, the time required by a recently developed algorithm [72] based on Metropolised independent sampling to generate a valid sample scales as where ã relates to the speed of the classical computer.This value has been reported to be ã = 3 × 10 −15 s in the case of the Tianhe-2 supercomputer [100].For the case of an ordinary computer, we choose this value to be ã = 3 × 10 −9 s.

APPENDIX D: BENCHMARK WITH EXACT SIMULATIONS
We derive a master equation model in order to benchmark the sampling rate formula in Eq. ( 6).In our model, an N -particle initial state is given by |ψ 0 ⟩, where ⟨ψ 0 |ψ 0 ⟩ = 1 and where n s is the number operator counting the number of atoms in the lattice site s.In Eq. (D2), the first and second terms account for single-body and twobody decay processes, respectively.The effect of V (t step ) is to leave the state in the N -particle subspace (i.e., the collision-free subspace in which we are interested), but with a reduced norm.When the operator is applied to a state with k pairs, it reduces the norm by a factor exp (−N t step /2τ BG ) exp (−kt step /2τ TB ).The survival probability of the boson sampler is thus characterized by the norm squared, ⟨ψ|ψ⟩.
We use the master equation model given in Eq. (D1) to carry out numerical simulations for up to N = 4 bosons, M = N 2 modes, and M steps.The exact numerical simulations allow us to validate the survival probability model in Eq. (B3), which was used to estimate the sampling rate of an atom boson sampler in Sec.III.We consider in particular the survival probability from time step j−1 to j, which is given by p j = ⟨ψ j |ψ j ⟩/⟨ψ j−1 |ψ j−1 ⟩, where |ψ j ⟩ = j i=1 U i V (t step )|ψ 0 ⟩.In Fig. 6(a), we plot p j versus j and for different values of τ TB .For clarity, we omit the effect of single-particle losses as their effect is trivial, and for simplicity we choose as initial state a pure state given by a uniform superposition of all bosonic states of N particles.For τ TB = t exec , we find an almost exact correspondence between P step (black dashed line) and p j (blue solid curve) for all steps j, suggesting that, for weak losses, Eq. (B3) provides an accurate description of the overall survival probability.For lower values of τ TB , such as τ TB = t exec /5 or τ TB = t exec /20, this correspondence only holds for the initial state, while for the rest of the steps, P TB step < p j , showing that Eq. (B3) represents a lower bound on the survival probability for each step.In Fig. 6(b), we evaluate the total survival probability ⟨ψ M |ψ M ⟩ as a function of τ TB (blue solid curve), and compare it with the result from the simple model, (P TB step ) M .We find that the simple model gives a correct description of two-body losses when τ TB ≳ t exec and a lower bound estimate when τ TB ≲ t exec .Although numerical benchmarks have been performed for up to N = 4, to our knowledge, it is reasonable to assume that the observed results are generalizable to a much larger number of particles.

APPENDIX E: MONTE CARLO ANALYSIS OF THE HONG-OU-MANDEL EXPERIMENT
In order to extract P bunch from the experimental results, it is important to consider that atoms bunched in the same site cannot be directly detected by fluorescence imaging without losing them.In fact, independent measurements performed in our apparatus on bunched atoms show that with a probability P LIC,0 = 71(5) % no atom is left, while with a probability P LIC,1 = 1 − P LIC,0 just a single atom is detected in the final fluorescence image as a consequence of light-induced collisions [150].Hence, identical atoms lead to the detection of either zero or one atom, while distinguishable atoms are both detected if after the Hong-Ou-Mandel π/2 pulse they occupy different internal states.Examples of the three possible outcomes with 0, 1 and 2 atoms detected are shown in Fig. 4(a).
To measure the single-atom survival probability S, we performed independent measurements similar to the Hong-Ou-Mandel interference experiment outlined in the main text, omitting the microwave π/2 pulse.Without the π/2 pulse, P bunch =0, independent of the quantum purity of the state, γ.From this calibration experiment, we can directly determine S = 84(1) %.The value of S is limited by losses at the beginning of the transverse cooling process, which could be avoided in the future with an improved experimental procedure.
We employ a Monte Carlo analysis to more rigorously analyze the statistical outcomes recorded for the Hong-Ou-Mandel experimental sequence that is described in Sec.IV.The Monte Carlo simulation uses the predetermined light-induced collision probability P LIC,0 and the single-atom survival probability S. Further input parameters are the bunching probability P bunch , the addressing probability of the narrow bandwidth MW pulse, and the probability to successfully reconstruct the position of the atoms.Using a nonlinear least squares regression, we fit the generated Monte Carlo events to the measured data to extract the underlying experimental parameters, yielding P bunch = 73(6) %.

FIG. 1 .
FIG. 1. Illustration of a quantum circuit based on ultracold atoms in state-dependent optical lattices.Each site accommodates two modes, |±⟩, of the quantum circuit represented by two internal states of the atom.A representative initial state is shown, where every second mode of the circuit is occupied.State-dependent shift operations connect distant modes, while local operations couple the internal modes.Inset: Hong-Ou-Mandel interference of two indistinguishable atoms, where a combination of local phase-imprinting pulses and microwave pulses realize the equivalent of a generalized photonic beam splitter.

Figure 1
Figure 1 illustrates how to "wire" an arbitrary quantum circuit using ultracold atoms in state-dependent optical lattices.The idea here is to use the lattice sites along with two internal states of the atom, |+⟩ and |−⟩, to represent the modes of the quantum circuit, so that M/2 lattice sites accommodate M modes.Modes associated with FIG. 2.Scaling of boson sampling machines.Predicted sampling rate versus the number of bosons N is shown for atomic (blue), photonic (purple), and classical (grey) boson sampling machines.The atomic sampling is computed based on Eq.(6), where an extended model (see appendix B) is used for Pstep to account for both particle pairs and triplets.Our atomic sampling model disregards decoherence by technical noise sources; see text.The photonic and classical curves correspond to Eqs. (C1) and (C2), respectively.Under stateof-the-art conditions, atomic NISQ devices have the potential to overtake classical algorithms at N ≈ 37. Pioneering photonic experiments are marked with symbols[59],[60], and[61], while the best reported boson sampling experiment with[71].Note that the photonic experimental points refer to implementations of the original boson sampling problem, and thus do not include Gaussian boson sampling experiments, for the reasons presented in Sec.I. Our proofof-principle experiment, described in Sec.IV, is marked with symbol.The numerical results obtained with the master equation approach (see Appendix D) are indicated with symbols, with the vertical bars indicating the 1-σ statistical uncertainty.

FIG. 4 .
FIG. 4. Detection of the Hong-Ou-Mandel interference.Right: measured probability for the detection of zero (P0), one (P1), or two atoms (P2).The dashed lines show the case of distinguishable particles as a reference (see Appendix E).The error bars represent the 68 % Clopper-Pearson confidence intervals.Left: representative fluorescence images of the two atoms after the Hong-Ou-Mandel sequence, with the internal states |−⟩ and |+⟩ mapped to the sites 0 and 10.When the atoms are bunched in the same site, inelastic lightinduced collisions result in either one or zero atoms being detected.Super-resolution microscopy allows resolving the lattice sites despite the diffraction-limited optical resolution of ≈ 4 sites [123].

6 <
l a t e x i t s h a 1 _ b a s e 6 4 = " a l v I 4 i 5 K 1 9 W + 2 a 5 s O 8 A H i w h h p 0 0

9 <
e m x g / o 8 p w J n C U b 6 c a E 8 o G t I c t S y W N U P v Z 5 N Q R O b Z K l 4 S x s i U N m a i / J z I a a T 2 M A t s Z U d P X 8 9 5 Y / M 9 r p S a 8 8 D M u k 9 S g Z N N F Y S q I i c n 4 b 9 L lC p k R Q 0 s o U 9 z e S l i f K s q M T S d v Q / D m X 1 4 k 9 X L J O y u V 7 0 + L l e t Z H D k 4 h C M 4 A Q / O o Q K 3 U I U a M O j B M 7 z C m y O c F + f d + Z i 2 L j m z m Q P 4 A + f z B 6 Z l j W M = < / l a t e x i t > N = l a t e x i t s h a 1 _ b a s e 6 4 = " p L D q u h T p F H 5 L n z F z Y v e s h r B Q C G w = " > A A A B 6 3 i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 9 k t Y r 0 I R S + e p I L 9 g H Y p 2 T T b h i b Z J c k K Z e l f 8 O J B E a / + I W / + G 7 P t H r T 1 w c D j v R l m 5 g U x Z 9 q 4 7 r d T W F v f 2 N w q b p d 2 d v f 2 D 8 q H R 2 0 d J Y r Q F o l 4 p L o B 1 p Q z S V u G G U 6 7 s a J Y B J x 2 g s l t 5 n e e q N I s k o 9 m G l N f 4 J F k I S P Y Z N L 9 d a 0 + K F f c q j s H W i V e T i q Q o z k o f / W H E U k E l Y Z w r H X P c 2 P j p 1 g Z R j i d l f q J p j E m E z y i P U s l F l T 7 6 f z W G T q z y h C F k b I l D Z q r v y d S L L S e i s B 2 C m z G e t n L x P + 8 X m L C K z 9 l M k 4 M l W S x K E w 4 M h H K H k d D p i g x f G o J J o r Z W x E Z Y 4 W J s f G U b A j e 8 s u r p F 2 r e p f V 2 s N F p X G T x 1 G E E z i F c / C g D g 2 4 g y a 0 g M A Y n u E V 3 h z h v D j vz s e i t e D k M 8 f w B 8 7 n D x T Z j Z 0 = < / l a t e x i t > N = 27 < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 F f M Z i Z y a p t R 5 j w 3 r C q C V k 9 4 0 2 g = " > A A A B 6 n i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h d 0 o 6 k U I e v E k E c 0 D k i X M T n q T I b O z y 8 y s E J Z 8 g h c P i n j 1 i 7 z 5 N 0 4 e B 0 0 s a

M
m=1 n m |ψ 0 ⟩ = N |ψ 0 ⟩.Here, n m is the number operator acting on mode m.The evolution of such a quantum state is given by|ψ⟩ = U M V (t step )U M −1 . . .U 1 V (t step )|ψ 0 ⟩,(D1)where U j are the coherent operations at step j, and V (t) = exp (−Ht) represents particle decay operator acting for a time t.The Hermitian operator H is

FIG. 6 .
FIG. 6. Numerical simulation benchmark of the survival probability.Blue lines represent the exact simulation for N = 4 and M = 16, where each point is the average of 30 realizations with different random unitaries.(a) Success probability pj from time step j − 1 to j for different decay rates τTB = texec, τTB = texec/5 and τTB = texec/20.The black dashed lines are obtained using the simple model in Eq. (B3).(b) Total survival probability ⟨ψM |ψM ⟩ versus τ −1 TB .The black dashed line represents the simple model, (P TB step ) M .