Quantum computational supremacy in the sampling of bosonic random walkers on a one-dimensional lattice

We study the sampling complexity of a probability distribution associated with an ensemble of identical noninteracting bosons undergoing a quantum random walk on a one-dimensional lattice. With uniform nearest-neighbor hopping we show that one can efficiently sample the distribution for times logarithmic in the size of the system, while for longer times there is no known efficient sampling algorithm. With time-dependent hopping and optimal control, we design the time evolution to approximate an arbitrary Haar-random unitary map analogous to that designed for photons in a linear optical network. This approach highlights a route to generating quantum complexity by optimal control only of a single-body unitary matrix. We study this in the context of two potential experimental realizations: a spinor optical lattice of ultracold atoms and an quantum gas microscope.

What level of quantum complexity cannot be efficiently simulated on a classical computer? This is a central question in quantum information science that impacts fundamental physics from the nature of condensed matter [1][2][3] to the geometry of space time [4]. While it is widely believed that a universal fault-tolerant quantum computer achieves such complexity, implementation of such a device is still a distant prospect. Nonetheless, more modest devices designed for a limited task could supersede the power of a classical computer and achieve so-called "quantum supremacy" [5,6]. In particular, a quantum device can yield random outcomes sampled from a probability distribution such that no classical computer could efficiently simulate its statistics. In their seminal paper, Aaronson and Arkhipov [7] showed that quantum supremacy could arise from "sampling complexity" in the most unlikely of places: linear optics. Based on the hierarchy of computational complexity classes, we believe that the number distribution of identical photons scattering from a linear optical network cannot be efficiently simulated. In this case, the nonclassical behavior arises solely from the quantum statistics of identical particles, or equivalently, the entanglement between modes, since photons are noninteracting particles in linear optics. The study of Boson Sampling that achieves quantum supremacy is a goal pursued worldwide [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].
In this work we extend the physical paradigm for studying Boson Sampling to a continuous-time quantum random walk of noninteracting identical bosons on a 1D lattice [28][29][30][31][32][33]. Our motivation is two-fold. Firstly, we seek to understand the minimal complexity necessary to achieve quantum supremacy. A 1D gas of noninteracting bosons is perhaps the simplest nontrivial many-body quantum system. Secondly, this paradigm can be realized in physical platforms that might lead to more scalable implementations, e.g., ultracold bosonic atoms in an optical lattice. Recent experiments with atomic gas microscopes demonstrate the possibility is observing controlled transport of modest sized ensembles of identical * gopu90@unm.edu atoms and direct readout of the atom number distribution at the lattice sites [34][35][36][37][38][39]. In contrast to linear optics, imperfections such as photon loss, misalignment of a large interferometer, detection efficiency, and the challenge of producing a large input (∼ 50) of identical bosons might be mitigated, improving the prospect of near-term demonstration of quantum supremacy.
Recent work has explored the complexity of the dynamics of bosonic atoms hopping in optical lattices. Deshpande et. al [40] studied nearest-neighbor hopping in hypercubic lattices in d-dimensions. They obtained upper and lower bounds for the scaling of time t with the number of bosons for which the problem is efficiently simulatable or hard to classically simulate, respectively. Elben et al. [41,42] studied the generation of pseudorandom unitary maps (t-designs) using random quenches on atoms trapped in lattices, in order to measure Renyi entropies in a Bose-Hubbard model. Our work focuses solely on 1D non-interacting Bose gases, with the goal of understanding the minimal complexity necessary for achieving quantum supremacy. We show that even a time-independent translationally invariant system system has complexity for which there is no known efficient classical algorithm, and that with time-dependent global control of the hopping Hamiltonian, we can engineer an arbitrary single-body transition matrix, and thus achieve the same complexity as an arbitrary linear optical network.
In Boson Sampling one considers N identical noninteracting bosons in M modes evolving by the linear map, Given an input Fock state n in = n in 2 , . . . , n in M , the transition probability to a corresponding output Fock state with l n out l = l n in l = N is Where Λ n out |n in is a "submatrix" obtained by repeating n in l times the l th column of Λ, and n out l times its l th row. Perm(Λ) is the permanent of the matrix. We restrict our arXiv:1805.01858v2 [quant-ph] 14 May 2018 attention here to a fiducial input state with 1 boson in each of N < M input modes and 0 bosons in the remaining M − N modes. A Boson sampler then outputs a random sequence n out with probability P (n out ). The complexity of Boson Sampling from this distribution on a classical computer depends on the complexity of calculating the permanent of Λ and its submatrices [43].
In the optical Boson Sampling setting, U is the Smatrix for scattering incoming photons into the outgoing modes of a linear optical network [7]. There, an arbitrary unitary matrix Λ l l can be constructed through a sequence of phase shifters and beam splitters [44]. Here we consider a continuous time quantum random walk of identical bosons on a lattice, described by the second quantized Hamiltonian H = l ,l h l l a † l a l , where h l l is the generator of tunneling (hopping) from one lattice site l to another l . The unitary evolution of creation operators is exactly the same as for linear optics, except that the S-matrix is now replaced by the time evolution operator, U (t) = l ,l u l l (t)a † l a l , where u l l (t) = e −iht l l is the single particle time evolution operator in the lattice site basis. Setting Λ l l = u l l (T ) for the final time T , we see the isomorphism between the quantum random walk and the linear optical network.
If the Hamiltonian h l l describes nonlocal hopping on an arbitrary graph, then one can obtain an arbitrary Λ. We restrict our attention here to nearestneighbor hopping. In the simplest case of uniform hopping h l l = −J(δ l ,l+1 + δ l ,l−1 ); here we consider periodic boundary conditions. The Hamiltonian is trivially diagonalized in the Bloch basis with creation operators M a † l , yielding the band structure in the tight-binding model and resulting time evolution operator While the time evolution is trivial, the same cannot be said in general for sampling complexity. In particular, sampling is basis dependent, because the permanent is. Transforming back to the lattice site basis, we have The protocol for Boson Sampling in this system is then as follows. A single boson is prepared at each of N sites of an M -site periodic lattice, each in an identical vibrational level. After some a final time T we measure the number of bosons in each of the sites, which provides a sample from a probability distribution that depends on the permanent of Λ ll = u ll (T ). We seek to understand the computational complexity of calculating this permanent. The matrix Λ ll has a structure imposed by the restricted nature of the Hamiltonian. It is a circulant ma- trix whose elements depend only on (l − l ) mod M . In a translationally invariant system, the transition amplitude depends only on the ballistic distance traveled. A plot of the transition amplitudes in Fig. 1 shows the familiar wave function for a continuous-time quantum random walk [45]. Note that the elements of the matrix decays exponentially after a "band" of elements, and again increase towards edges due to the boundary conditions. We call this band B, which is proportional to T . To first approximation we set Λ ll = 0 for |Λ ll | < ||Λ|| for some 1, where |||Λ|| is a norm of the matrix. The transition matrix then takes the form of a doubly banded matrix.
In the case where M is sufficiently large compared to N , and all of the bosons are trapped in the central region of the lattice, a submatrix with dimension M sufficiently smaller than M is a banded Toeplitz matrix. Schwartz [46] devised an algorithm to calculate the permanent of such matrices which scales as If we include all modes M , assuming a lattice on a ring, the transition matrix Λ is a cyclically banded matrix [47]. Cifuentes et. al [48], showed that one can exactly determine the permanent in of such matrices in O(N 2 2B ). This algorithm does not use the fact that the matrix has a circulant structure. We conjecture in this case that the algorithm can be improved to a scaling of O(log N 2 B ).
We address the complexity of Boson Sampling in this geometry in two regimes. If time of evolution is constant with system size, T ∼ O(1), so is B. The algorithms discussed then scale as O(N ) for the Cifuentes algorithm [48] and O(log N ) in the Schwartz algorithm [46]. Thus, for constant time of evolution (or at most logarithmic in N ), one can devise an algorithm to efficiently sample from the resulting occupation number probability distribution of quantum random walkers. On the other hand, if time of evolution grows faster than T ∼ O(log N ), (say polylog(N ) ) the scaling of the algorithms in [46] and [48] becomes super-polynomial (exponential in case of T ∼ O(N )).
While the known algorithms for calculating the permanent or even approximating it to within a multiplicative error scale poorly, we cannot rule out the possibility of the existence of an efficient classical algorithm. Nonetheless, given the growth in the graph structure for these matrices [48], we conjecture that this is indeed a hard problem, and noninteracting quantum random walkers show quantum complexity even when restricted to a 1D lattice with uniform nearest-neighbor hopping. As further evidence, consider a simulation of the many-body system based on matrix product states (MPS) of bosonic modes [49,50]. We find that the bond dimension of the MPS has a similar scaling to the algorithms described above. For constant time the bond dimension is logarithmic in system size and entanglement growth is polynomial. For a final time growing linearly with systems size the MPS entanglement grows exponentially.
It is important to note that quantum supremacy does not require us to sample from the exact probability distribution P . It is sufficient that the distribution P one samples from is a multiplicative approximation to P , i.e., |P − P | < P . We know of no multiplicative approximation here. Other classical sampling algorithms, such as Monte-Carlo sampling, will only allow one to calculate the permanent to within an additive approximation.
Since uniform 1D hopping yields banded Toeplitz transition matrices and no definitive hardness result, a more direct path to observe quantum supremacy is to connect to the arguments of Aaronson and Arkipov [7]. For that purpose we seek to generate a Haar random unitary transition matrix in an experimental realistic model with only nearest neighbor hopping. By employing the tools of optimal control [51,52] we can engineer any desired target unitary map on the system. This is feasible in a scalable way as we need only design the one-body unitary map u ll (T ), an M × M matrix, where M is the number of modes included in the model. Quantum supremacy is then achieved in the many-body system solely due to the quantum statistics of identical particles. This is in contrast to the similar efforts to observe quantum supremacy in quantum circuits with N qubits [53,54]. There, one seeks to implement sufficiently random 2 N × 2 N matrices, a task one cannot implement scalably via quantum control.
To achieve full control, we generalize here to spinoroptical lattices, a platform long studied for implementations of controlled transport [55][56][57][58][59][60][61][62][63]. Specifically, consider two optical lattices with circular σ ± polarization, whose nodes can be displaced by a parameter θ, e.g., in a lin-angle-lin lattice [55,56,58], or through direct synthesis of two independent standing waves [63]. Two chosen internal atomic states, labeled |↑ and |↓ can be trapped separately in each state-dependent lattice. For concreteness we consider two magnetic sublevels in the electronic ground state separated in energy by the hyperfine splitting, e.g., the stretched states, |↑ = |F = 4, m F = 4 , |↑ = |F = 3, m F = 3 , in cesium [58]. Microwaves   FIG. 2. Coherent transport in a spinor optical lattice. A microwave field drives spin flips between two different hyperfine magnetic sublevels |↑ and |↓ , correlated with transport of the atomic wavepacket in the lowest Wannier states for optical lattices with σ± polarization. The transition matrix element is determined by the wavepacket overlap (sketched below). The amplitude of hopping left(right) Ω L(R) is controlled by the microwave amplitude and the time-dependent angle θ, which determines the relative displacement of the two standing waves; the phase of the hoping matrix element is set by the time-dependent microwave phase φ. The overall translational symmetry is broken by the addition of a quadratic chemical potential. that drive spin flips between these internal states are then correlated with transport between the sublattices. We envision deep lattices along 1D so that tunneling of atoms within a given standing wave is highly suppressed and quantum walks occur only due to microwaveinduced hopping in the tight-binding approximation (see Fig. 2). We neglect here any interaction between the atoms, a good approximation if the transverse confinement is weak, or by tuning the interaction to zero via a Feshbach resonance [64].
Our goal is to generate an arbitrary unitary transport map Λ ll . In order to break the translational symmetry we add an additional quadratic light-shift potential. With this addition we show below that the system is controllable when the relative phase of the standing waves and phase of the microwave are time-dependent control parameters. We take as our control Hamiltonian H 0 is the time-independent "drift Hamiltonian," taken here to be a quadratic potential with l ↓ = (N + 1)/2, l ↓ = (N + 1)/2 + θ/2. H L and H R govern spindependent transport hopping to the left and right, driven by microwaves, with lattice labeling convention shown in Fig. 2. The hopping amplitudes are set by the global microwave Rabi frequency Ω 0 and the spatial overlap of atomic wavepackets in the tight-binding approximation. For deep lattices, where η = 2πx 0 /λ is the Lamb-Dicke parameter, x 0 the width of the vibrational ground state Gaussian wavepacket, and λ is the wavelength of the lattice [58]. Additionally, the phase of the hopping amplitude is set by the global phase of the microwave, φ(t). We do not include any controls that address individual lattice sites.
Our aim is to show that the above Hamiltonian is controllable in the one-body Hilbert space of a d-dimensional lattice. We do so by employing the tools of geometric control [65][66][67]. One must prove that the set {H 0 , H R , H L } form the generators the Lie algebra of the group SU (d)the unitaries that describe all one-body transition matrices (up to an overall phase). A detailed proof is given in the supplementary material, where we show that all linear combinations of multiple commutators of that generating set span the d 2 −1 generalized Gell-Mann matrices, a basis for su(d).
Given that the system is controllable in the sense described above, we use numerical optimal control to find time-dependent waveform that implements a desired unitary map. Generally, the goal is to optimize the fidelity between the unitary map arising from the timedependent Schödinger equation and the target transition matrix, F[θ, φ] = |Tr Λ † u[θ, φ, T ] | 2 /d 2 over the functions θ(t) and φ(t). We employ a variant of the well known GRAPE algorithm [68], building on the methods that have been employed in the Jessen group to design implement arbitrary SU(16) maps on the internal hyperfine levels [52]. We parameterize the control waveforms as piecewise constant functions so that the Hamiltonian is a function of the control vector b = {θ(t k ), φ(t k )|k = 1, 2, . . . , T /(δt)}, where δt is the duration of each step, so that the number of steps is T /δt. Using GRAPE we optimize F(b) yielding the control vector b.
The numerical simulation is done by first numerically generating a target Haar random unitary. The parameters are chosen as follows, η = 0.4, Ω L ( π 2 ) = Ω R ( π 2 ) = 1, V 0 = 1, δt = 2π. The total number of time steps is 2N 2 , so the total time is 2d 2 δt. The dimension of parameter space is 4d 2 . The target unitary, analytic form of the derivative of fidelity function, and an initial guess for control vector b are given as inputs to the GRAPE algorithm. The algorithm does a gradient descent by varying control vector and find one that minimizes infidelity. Numerical results show infidelity as low as 10 −11 .
For ease of numerics and controllability, we have employed periodic boundary conditions, but this does not impose physical constraints. We initialize the system with N atoms in d lattice sites such that d > N 2 . The optimization algorithm is restricted to a target unitary M × M Λ matrix in an subspace of dimension M = N 2 . If we trap the atoms in fixed N lattice sites, the optimization will be a partial isometry with N orthonormal columns in d-dimension. A general control waveform to achieve this target will require dN time steps. If d is chosen to be sufficiently larger than M = N 2 , then the amplitudes will not reach the boundary in any intermediate step, and therefore will not have any boundary effects.
Optimal control can never yield a fidelity exactly equal to one, and thus a perfect waveform to achieve a target Haar-random unitary is never possible. In addition, experimental errors will always occur in any implementation of a control waveform. However, as long as the final probability distribution is close to the ideal one in total variation distance, we can still guarantee its hardness [7]. Here, as long as the error that arises as the infidelity of generated unitary in our protocol is kept constant with respect to N , we can guarantee the hardness of output distribution [69].
In this letter we have shown that even the simplest trivially solvable many-body system, noninteracting Bosons in a 1D lattice with uniform nearest neighbor hopping, could exhibit quantum supremacy in sampling complexity. Moreover, if we extend this 1D system to include time-dependent control, one can generate an arbitrary a one-body Haar-random unitary transition matrix. We studied this for the case noninteracting atoms undergoing microwave-driven quantum walks in a spinor optical lattice. This system can achieve the same quantum supremacy as Boson Sampling of photons in an arbitrary linear optical network. While we have used quantum optimal control to reach a desired target Haar-random unitary, the complexity of control required to achieve quantum supremacy remains an open question. In particular, a random control waveform will yield a pseudorandom transition matrix that after a suitably long mixing time [70], and we expect this time to be much shorter than that required to reach a particular target Haar random unitary target. This could lead to a more efficient and scalable approach to Boson Sampling. Finally, of interest is the nature of sampling complexity when interactions between bosons are included. This connects Boson Sampling to the more general problem of quantum simulations of many-body systems, such as the well studied Bose-Hubbard model [71].
We thank Tyler Keating, Christopher Jackson, Rolando Somma, and Saleh Rahimi-Keshari for useful discussions. This work was supported in part by National Science Foundation grants PHY-1521016 and PHY-1521431.
After these manipulations, we have isolated the following terms in the Lie algebrâ Since we already have the termÂ N (φ) in our set, we obtain the termĈ =Â 1 (φ) +Â N −1 (φ), from the commutator.
Ĉ ,Ẑ 2 (θ = 0) ⇒ (2 − N )Â 1 (φ) + NÂ N −1 (φ) Combining the two equations above, we get the termŝ A 1 (φ), andÂ N −1 (φ) separately. It further follows that, If we repeat the above steps we obtain all terms of the formÂ (φ) ⊗ σ z = | + 1 |e −iφ + | + 1|e iφ ⊗ σ z . As we have already shown that this is sufficient to generate the Gell-Mann matrices of a system of spinor optical lattices of dimension 2N . tion. Fig. 3 shows the average infidelity achieved for random target unitaries for different dimensions, which can be as small as 10 −11 for low dimensions. Fig. 4 shows a sample control waveform for a random target unitary in dimension d = 2N = 10, which appear random, but is in fact of a form that drives the system to implement the desire unitary transition matrix Λ.