Simulating the Dynamics of Single Photons in BosonSampling Devices with Matrix Product States

TL;DR

This paper demonstrates how matrix product states can be used to simulate and analyze the quantum dynamics of photons in BosonSampling devices, revealing entanglement growth and extending to dissipative and fermionic systems.

Contribution

It introduces a novel application of matrix product states to simulate photon dynamics in BosonSampling, including dissipative and fermionic circuits, enhancing understanding of quantum interference.

Findings

Matrix product states effectively simulate photon evolution in BosonSampling.

The method reveals entanglement entropy growth during photon dynamics.

Applicable to dissipative and fermionic quantum circuits.

Abstract

BosonSampling is a well-defined scheme for demonstrating quantum supremacy with photons in near term. Although relying only on multi-photon interference in nonadaptive linear-optical networks, it is hard to simulate classically. Here we study BosonSampling using matrix product states, a powerful method from quantum many-body physics. This method stores the instantaneous quantum state during the evolution of photons in the optical quantum circuit, which allows us to reveal the dynamical features of single photons in BosonSampling devices, such as entanglement entropy growth. We show the flexiblility of this method by also applying it to dissipative optical quantum circuits, as well as circuits with fermionic particles. Our work shows that matrix product states is a powerful platform to simulate optical quantum circuits. And it is readily extended to study quantum dynamics in multi-particle quantum walks beyond BosonSampling.

Figures

Equations

$$|\psi ⟩=\sum _{n_1,n_2,\dots ,n_M}{c}_{n_1,n_2,\dots ,n_M}|n_1,n_2,\dots ,n_M⟩,$$

$$c_{n_1,n_2,\dots ,n_M}=\sum _{a_1,\dots ,a_{M+1}}{B}_{a_1,a_2}^{n_1}{B}_{a_2,a_3}^{n_2}\dots B_{a_M,a_{M+1}}^{n_M}.$$

$$\sum _{a_{i+1},n_i}{B}_{a_i,a_{i+1}}^{n_i}\mathrm{conjugate}(B_{a_i^{\prime },a_{i+1}}^{n_i})={\delta }_{a_i,a_i^{\prime }},$$

$$B_{a_i,a_{i+1}}^{n_i^{\prime }}=\sum _{n_i}{R}_{n_i}^{n_i^{\prime }}{B}_{a_i,a_{i+1}}^{n_i}.$$

$${\Phi }_{a_i,a_{i+2}}^{n_i,n_{i+1}}=\sum _{a_i,a_{i+1}}{S}_{a_i,a_i}{B}_{a_i,a_{i+1}}^{n_i}{B}_{a_{i+1},a_{i+2}}^{n_{i+1}},$$

$${\Phi }_{a_i,a_{i+2}}^{n_i^{\prime },n_{i+1}^{\prime }}=\sum _{n_i,n_{i+1}}{\Phi }_{a_i,a_{i+2}}^{n_i,n_{i+1}}{O}_{n_i,n_{i+1}}^{n_i^{\prime },n_{i+1}^{\prime }},$$

$$\sum _{a_{i+1}^{\prime }}{U}_{a_i,a_{i+1}^{\prime }}^{n_i^{\prime }}{S}_{a_{i+1}^{\prime },a_{i+1}^{\prime }}{V}_{a_{i+1}^{\prime },a_{i+2}}^{n_{i+1}^{\prime }}=\mathrm{svd}({\Phi }_{(a_{i+2},n_{i+1}^{\prime })}^{(a_i,n_i^{\prime })}),$$

$$\sum _{a_i,n_i}{U}_{a_i,a_{i+1}}^{n_i}\mathrm{conjugate}(U_{a_i,a_{i+1}^{\prime }}^{n_i})={\delta }_{a_{i+1},a_{i+1}^{\prime }}.$$

$$B_{a_{i+1}^{\prime },a_{i+2}}^{n_{i+1}^{\prime }}$$

$$B_{a_i,a_{i+1}^{\prime }}^{n_i^{\prime }}$$

$$\sum _{i=1}^M{n}_i=N.$$

$$|\varphi ⟩=|1{⟩}^N\otimes |0{⟩}^{M-N}.$$

$$U_{\mathrm{ps}}^b(\varphi )|n⟩=e^{-2\pi \mathrm{i}n\varphi }|n⟩.$$

$$\displaystyle U_{{\rm bs}}^{b}(\theta)\left(\begin{array}[]{c}\hat{a}^{\dagger}\\ \hat{b}^{\dagger}\end{array}\right)=\left(\begin{array}[]{cc}\cos(2\pi\theta)&-\sin(2\pi\theta)\\ \sin(2\pi\theta)&\cos(2\pi\theta)\end{array}\right)\left(\begin{array}[]{c}\hat{a}^{\dagger}\\ \hat{b}^{\dagger}\end{array}\right)$$

$$U_{\mathrm{bs}}^b(\theta )|m,n⟩=$$

$$=$$

$$\frac{{\left(\sin (2\pi \theta ){\hat{a}}^{†}+\cos (2\pi \theta ){\hat{b}}^{†}\right)}^n}{\sqrt{n!}}|0,0⟩$$

$$U_{\mathrm{bs}}^b(\theta )|1,1⟩=$$

$$-\sqrt{2}\cos (2\pi \theta )\sin (2\pi \theta )|0,2⟩$$

$$+\left({\cos }^2(2\pi \theta )-{\sin }^2(2\pi \theta )\right)|1,1⟩.$$

$$U_{\mathrm{bs}}^f(\theta )|1,1⟩=|1,1⟩,$$

$$\rho ={\sigma }^N\otimes |0⟩⟨0{|}^{M-N},$$

$$|\rho ⟩⟩=|\sigma ⟩{⟩}^N\otimes |00⟩{⟩}^{M-N}.$$

$$U^o=\mathrm{kron}(U,\mathrm{conjugate}(U)).$$

$$\sigma =\bar{n}|1⟩⟨1|+(1-\bar{n})|0⟩⟨0|,$$

$$S=-\mathrm{tr}\left[{\rho }_A\log ({\rho }_A)\right],{\rho }_A={\mathrm{tr}}_B\rho .$$

$$S^o=-{\mathrm{tr}}_A\left[R\log (R)\right],R=⟨⟨\rho ||\rho ⟩{⟩}^{-1}{\mathrm{tr}}_B|\rho ⟩⟩⟨⟨\rho |.$$

$${ϵ}_j^b$$

$${ϵ}_j^f$$

Full text

I Introduction

Multi-photon interference lies at the heart of optical quantum information processing pan2012 ; sihui2012 , including optical quantum computing kok2007 , quantum cryptography lo2012 , quantum simulation Guzik2012 , and quantum metrology Giovannetti2011 . Based on some basic interference phenomena, such as the Hong-Ou-Mandel (HOM) effect hong1987 , it is possible for us to make use of quantum advantages to realize some complex processing tasks beyond the current classical limit Lund2017 , by increasing the number of single photons and the complexity of the optical quantum circuit.

In particular, BosonSampling Aaronson2013 , an intermediate model of linear optical quantum computing, utilizes only passive multi-port optical interferometer for the evolution of input single-photon sources, and then samples the output events using single-photon detectors, which is a well-defined problem that can be naturally and efficiently solved on a specialized photonic quantum simulator. However, it has been shown to be classically intractable because the amplitude associated with each term in the output quantum state of BosonSampling is related to the calculation of a matrix permanent Aaronson2013 , which is known to be a $\#{P}$-hard problem Valiant1979 . Thus, BosonSampling is considered as a promising candidate to experimentally demonstrate the $quantum$ $supremacy$ in the near future Harrow2017 . So far, a number of elegant BosonSampling experiments has been achieved with linear optics on a small scale Broome2013 ; Tillmann2013 ; Crespi2013 ; Spring2013 ; Spagnolo2014 ; Carolan2014 ; Carolan2015 ; he2017 ; wang2017 ; Loredo2017 ; Loredo2017 ; wang2018 . In addition, BosonSampling can be regarded as a special multi-particle quantum walk, and thus investigating its nonclassical correlations during the quantum dynamics is also an interesting task Urbina2014 ; Tamma2015 ; Laibacher2018 ; Shchesnovich2016 ; Walschaers2016 ; Huang2017 .

Previous works which use classical algorithms to solve the sampling task for BosonSampling have mostly focused on pushing the boundaries of quantum supremacy, namely to explore what is the largest scale that is reachable by the best classical algorithm and the best supercomputer in the world. In Neville2017 , it is reported that classical simulation of outputting one sample of a $30$-photon BosonSampling takes half an hour on a laptop, using the Metropolised independence sampling algorithm (MIS). In WuYang2018 , the permanent of a matrix with size $48\times 48$ was calculated using Ryser’s algorithm Ryser1963 , which takes $196608$ cores and an execution time of around $4500$ seconds on Tianhe-2 supercomputer. Recent work also show that BosonSampling does not necessarily require to compute a huge number of permanents to sample from the output distribution CliffordClifford2017 .

In this work, we focus on simulating the evolution of photons in the optical quantum circuit for BosonSampling (or, more broadly, the dynamical behavior of multi-particle quantum walks), rather than simulating the sampling tasks like those in the previous works. By simulating the behavior of photon in complex linear optical quantum circuits, it is possible to measure intermediate observables such as correlations, entanglement entropy, to have a better understanding of the underlying dynamical process, which is important for the study of multi-particle quantum walk. However, it is difficult to store the entire information about the evolution of multi-particle in a optical quantum circuit with $M$-mode and $N$-input photons. Firstly, considering the size of the Hilbert space it can easily be seen that, the number of parameters \left({\begin{array}[]{*{20}{c}}{M+N-1}\\ {M-1}\end{array}}\right) of the output quantum state at each depth of the optical quantum circuit grows exponentially with the number of photons and modes in the system, thus one would soon run out of memory if one represents such a quantum state exactly on a classical computer. Secondly, for some special cases, such as BosonSampling, calculating a single amplitude of the output quantum state requires $O(2^{n}n^{2})$ runtime according to Ryser’s algorithm Ryser1963 .

We use matrix product states (MPS) method with a global U(1) symmetry (particle number conservation symmetry) to directly simulate the evolution of photons in the optical quantum circuit for BosonSampling. MPS has been a very successful numerical tool in solving one-dimensional quantum many-body problems Perez-GarciaCirac2007 ; Schollwock2011 , for both the unitary systems White1992 ; White1993 ; DelayVidal2004 and the dissipative systems GuoPoletti2018 ; BernierPoletti2018 ; XuPoletti2018 . Briefly speaking, MPS works by compactly rewriting a one-dimensional quantum manybody state as a product of many low dimensional tensors. For certain types of quantum systems which obey the so-called area law Hastings2007 , the sizes of the low dimensional tensors are bounded, which makes MPS to be a very memory efficient representation of the quantum state. The time complexity of MPS is also low in that the local unitary operations can be performed on local tensors only, without affecting other tensors. Moreover, measuring observables, especially local observables with MPS can be performed extremely efficiently Schollwock2011 .

For quantum systems with particle number conserving symmetry, MPS algorithms will usually be accelerated if the quantum symmetry is explicitly takeing into account, namely using a U(1) symmetric MPS SinghVidal2010 ; GuoPoletti2019a . We note that a non-symmetric matrix product states method has been applied to study time-bin quantum optics in Lubasch2018 , without using the number conserving property of the underlying optical quantum circuit. In particular, we highlight that we have developed a generic platform to study optical quantum circuits based on the U(1) symmetric MPS. It can be used to simulate different particle types including both the bosons and the fermions. Simulating an open quantum system can be done almost in parallel with the unitary case. The changes of the entanglement entropy as well as other physical observables during each depth of the evolution can be measured efficiently with this method.

The paper is organized as follows: In Sec. II, we briefly introduce the MPS algorithms and how it is used to study optical quantum circuits. In Sec. III, we show the entanglement entropy of the quantum state in a bosonic circuit is much larger and grows much faster than that in a fermionic circuit, and how the numerical complexity of simulating BosonSampling using MPS method is affected by losses. Finally we summarize in Sec. IV.

II method

In this section, we briefly show how to use MPS to simulate many-particle dynamics in various cases, namely bosons and fermions in an ideal environment as well as in a lossy environment.

A quantum many-body state living on a bosonic or fermionic optical lattice of size $M$ ($M$ modes) can be written as

[TABLE]

where $n_{i}$ labels the local Fock state on the $i$-th site of the lattice, $n_{i}$ is unbounded for bosons and $n_{i}=0,1$ for fermions. The coefficient $c_{n_{1},\dots,n_{M}}$ is an $M$ dimensional tensor. Assuming the size of the local Hilbert space is $l$, then the tensor $c$ would contain $l^{M}$ complex numbers if quantum symmetry is not taken into account. MPS rewrites $c$ as a product of $3$ dimensional tensors:

[TABLE]

Here we have used $a_{i}$ to denote the auxiliary degree of freedom, with $a_{1}=a_{M+1}=1$ added for notational convenience. The largest size of $a_{i}$ is referred to as the bond dimension $D$, which essentially characterizes the complexity of MPS based algorithms.

MPS keeps invariant when inserting a pair of matrices $U$, $U^{-1}$ between a pair of sites $i$ and $i+1$ (such a pair of neighbouring sites is also referred to as a bond $i$). To fix this gauge degree of freedom, one can prepare the MPS in canonical forms Schollwock2011 . In our simulation we keep MPS in the right canonical form, which means that each tensor $B_{a_{i},a_{i+1}}^{n_{i}}$ satisfies the condition

[TABLE]

where ${\rm conjugate}(a)$ means taking the element-wise conjugate of the tensor $a$.

Local unitary transformations on the quantum state in Eq.(1) can be easily translated into operations on MPS. Here we exemplify this by showing how to apply one-body and nearest-neighbour two-body operations to MPS. The application of one-body operations $R_{n_{i}}^{n_{i}^{\prime}}$ on site $i$ is simply

[TABLE]

We can see that one-body operation will not break the right canonical condition in Eq.(3) as long as $R_{n_{i}}^{n_{i}^{\prime}}$ is unitary.

For a two-body operation $O_{n_{i},n_{i+1}}^{n_{i}^{\prime},n_{i+1}^{\prime}}$ acting on bond $i$, we first multiply $B_{a_{i},a_{i+1}}^{n_{i}}$, $B_{a_{i+1},a_{i+2}}^{n_{i+1}}$ and the singular vector $S_{a_{i},a_{i}}$ between sites $i-1$ and $i$

[TABLE]

then we apply $O_{n_{i},n_{i+1}}^{n_{i}^{\prime},n_{i+1}^{\prime}}$ to $\Phi$ and get

[TABLE]

after that we do a singular value decomposition (SVD) to $\Phi$

[TABLE]

where we have grouped the indexes $(a_{i+2},n_{i+1}^{\prime})$ and $(a_{i},n_{i}^{\prime})$ so that SVD is performed on a two-dimensional tensor. From the definition of SVD, the tensor $V$ satisfies the right canonical condition in Eq.(3) and the tensor $U$ satisfies the left canonical condition

[TABLE]

With these results, we can update $B_{a_{i+1}^{\prime},a_{i+2}}^{n_{i+1}^{\prime}}$ and $B_{a_{i},a_{i+1}^{\prime}}^{n_{i}^{\prime}}$

[TABLE]

It is straightforward to verify that the new tensors $B_{a_{i+1}^{\prime},a_{i+2}}^{n_{i+1}^{\prime}}$ and $B_{a_{i},a_{i+1}^{\prime}}^{n_{i}^{\prime}}$ produced in Eqs.(9, 10) still satisfy the right canonical condition in Eq.(3). The division by very small singular values in Eq.(10) could make this algorithm numerically unstable, which could be overcome by the algorithm in Hastings2009 with negligible additional effort. We note that since we are using a U(1) symmetric MPS, all the tensors mentioned above are assumed to be generic tensors which are protected by a global U(1) symmetry. Such symmetric tensors can be simply taken as a list of dense tensors which are labeled by good quantum numbers, whose corresponding tensor operations mimic those of multi-dimensional dense tensors. The U(1) symmetric MPS leverages the fact that the total quantum number of the underlying quantum system is conserved, thus could often greatly reduce the memory requirements of MPS algorithms. For the implementation of generic tensor operations on U(1) symmetry protected tensors, one can refer to, for example SinghVidal2010 . The usage of U(1) symmetric MPS allows us to directly consider all the local degrees of freedom for bosonic particles, without the need to choose different truncations of the local Hilbert space until the result converges, as is done in Lubasch2018 .

In BosonSampling, only phase shifters and beam splitters are used. A phase shifter is a one-body operation, and a beam splitter is a two-body operation. Moreover, as in Fig. 1, the beam splitter acts only on nearest-neighbouring sites. Therefore the above operations would suffice for BosonSampling. We also note that the successive application of phase shifters and beam splitters to MPS would be analogous to the time-dependent matrix product states (t-MPS) algorithm Vidal2003 ; DaleyVidal2004 ; WhiteFeiguin2005 ; Schollwock2011 .

II.1 Bosonic circuit

One important property of a optical bosonic circuit made of phase shifters and beam splitters is that the total number of particles in the system is conserved. Namely if the input quantum state has a fixed total number of particles $N$, then during the evolution, $N$ is preserved. This means that we only need to consider the coefficients $c_{n_{1},n_{2},\dots,n_{M}}$ which satisfy

[TABLE]

Without loss of generality, the initial state $|\phi\rangle$ is chosen so that the first $N$ modes have one photon as inputs, while for the rest there is no input photon. This state could be written as

[TABLE]

To this end, we note that in the literature the phase shifter and beam splitter are often written as linear mappings in the operator space, namely mapping from one quantum operator to another. To use them with MPS, one need to rewrite them as quantum operators which map one quantum state to another. A phase shifter with a rotation angle $\phi$ can straightforwardly be written as

[TABLE]

For a beam splitter with an angle $\theta$, it is usually written in the operator space as

[TABLE]

where $\hat{a}^{\dagger}$ and $\hat{b}^{\dagger}$ are bosonic creation operators. Written as an operator on a quantum state $|m,n\rangle$, it would be

[TABLE]

For example, for an input state $|1,1\rangle$, Eq.(II.1) will be

[TABLE]

To ensure that the optical quantum circuit corresponds to a Haar random unitary matrix, we use the approach in RussellLaing2017 , namely each $\phi$ is generated according to a uniform distribution, while each $\theta$ is generated according to a specific distribution which is dependent on the depth and the position.

II.2 Fermionic circuit

In comparison with BosonSampling, sampling fermions (we call it FermionSampling in this paper) is related to calculating determinants of matrices, for which there exists efficient classical algorithms. Nevertheless, it would still be insightful to see whether FermionSampling is easy to simulate or not with MPS based algorithms. The simulation of FermionSampling is done by mapping the fermionic circuit to a one-dimensional spin chain by Jordan-Wigner transformation JordanWigner1928 ; LiebMattis1961 . As a result, the beam splitter for the fermionic operator $U_{{\rm bs}}^{f}$, in the zero- and one-particle sections, would remain the same as Eq.(II.1), while in the two-particle section, it is trivially

[TABLE]

due to the fermionic anti-commutation relation.

II.3 Lossy bosonic circuit

We also consider bosonic circuit in a lossy environment Garcia-PatronShchesnovich2017 ; OszmaniecBrod2018 ; RenemaGarcia-Patron2019 . In particular, we consider the simple case studied in Garcia-PatronShchesnovich2017 , where it is assumed that photons are injected into an optical circuit with uniform-loss channels of transmission ($i$.$e$. affecting all modes equally), and only the loss in the circuit is considered. In this simple and special case, the lossy linear optical circuit can be modeled by a lossless circuit, but with input to be a mixed state instead of a pure state. In this case, the quantum state is a mixture of different total number of particles. Nevertheless, by using the method in GuoPoletti2019a , we could still benefit from a U(1) symmetric MPS. Such an input state can be written as

[TABLE]

where each $\sigma$ is a local density operator. Here $N$ means the largest possible total number of particles instead. To deal with this kind of input state with MPS, one could reshape $\rho$ into a long vector which we denote as $|\rho\rangle\rangle$

[TABLE]

Here $|\sigma\rangle\rangle$ and $|00\rangle\rangle$ denote vectors resulting from squashing of the local density operators $\sigma$ and $|0\rangle\langle 0|$ in Eq.(24). The local Hilbert space would then be enlarged to $(N+1)\times(N+1)$. Correspondingly, any unitary gate $U$ operating on the corresponding unitary system, will become

[TABLE]

Here ${\rm kron}(a,b)$ means the kronecker product of two tensors $a$ and $b$.

We note the usage of U(1) symmetric MPS here allows us to choose much larger bond dimension $D$. In our simulations of the lossy cases, we have made use of this properties, which allows us to study an open system with $D$ larger than $10000$. In comparison, with a non-symmetric MPS, this case is almost intractable with a personal computer since the size of the local Hilbert space $l$ will already be $49$ for $N=6$ (A two-body operation would require to do an SVD on a dense matrix of size $Dl\times Dl$).

For simplicity, in this work we focus on a special type of initial local mixed state $\sigma$, which is

[TABLE]

$\bar{n}$ is the average filling of bosons which satisfies $0\leq\bar{n}\leq 1$. However, more general intial mixed state, as well as more general losses could also be studied with our setup.

III result

An important advantage of MPS is that one could explicitly monitor the growth of entanglement entropy $S$ during the dynamics, which is defined as

[TABLE]

Here the subscripts $A$ and $B$ mean two subsystems resulting from a bipartition of the system, and the density operator $\rho=|\phi\rangle\langle\phi|$ with $|\phi\rangle$ being the wave function of the system. In Fig. 2, we randomly generate $10^{3}$ optical quantum circuits corresponding to Haar unitary matrices for bosons and fermions with $M=20$ and $N=10$, and evolve them to a maximum depth $d_{\max}=M=20$. For all the simulations done in this work, we have chosen a maximum bond dimension $D=15000$ to ensure that it would neven be reached, and an SVD truncation threshold $t=10^{-8}$. In Fig. 2(a) we plot the distribution of entropy $S$ of the final quantum state as a function of the bond index $i$ corresponding to the bipartition position, for both the bosonic and the fermionic cases. We can see that in the bosonic case, the distribution is unbalanced (not reflection-symmetric around the central bond) which is a remnant of the unbalanced initial state, while in the fermionic case the distribution is balanced irrespective of the initial state. In in Fig. 2(b), we plot the growth of entropy $S$ as a function of the evolution depth $d$. The entanglement entropy $S$ of the quantum system strongly affects our ability to simulate it with MPS algorithms Schollwock2011 . From Fig. 2(a) and Fig. 2(b), we could see clearly that the entanglement entropy $S$ of the quantum state in a bosonic circuit is much larger and grow much faster than that in a fermionic circuit, which means that simulating the multi-particle evolution in a bosonic circuit is more difficult than that in a fermionic circuit using MPS. However, we must note that although the growth of entanglement entropy is the key factor can directly affect the efficiency of MPS Schollwock2011 , it is not equivalent to the classical hardness of BosonSampling, as there are cases with entanglement at the output that are classically simulatable Tichy2014 .

In fact, it is very simple and natural to measure $S$ during the dynamical process in the optical quantum circuit in Fig. 1. The singular values are stored alongside with the MPS whenever a two-body gate operation is performed, thus allowing to compute the entropy with almost no additional effort. We can simulate the evolution of photons in a BosonSampling circuit with $M=26$, $N=13$ and $d_{\max}=26$ in only about $5500$s on a personal computer with $2$ cores of 2.8 GHz frequency and $16$ Gb memory, where we have measured the entropy $S$ of each depth. We note that if we measure the entropy by first reconstructing the quantum state at each depth, we need to spend a huge amount of time to calculate all of the amplitudes of a quantum state. Moreover, it would require about $87$ Gb memory on a classical computer if one stores all these amplitudes exactly. Concretely, a pure python program using the Ryser’s algorithm is able to compute the permanent of a matrix of size $13$ in $0.09$s. Taking into considerations the huge number of matrix permanents to be computed to reconstruct the quantum state, MPS will be $10^{6}$ times faster than the Ryser’s algorithm. We note that here we only compare our method to a brute-force simulation that uses Ryser’s algorithm, but more efficient alternatives may exist.

We also apply MPS to study BosonSampling in presence of losses, which is called lossy BosonSampling Aaronson2016 . The hardness of MPS used to represent a density operator for an open quantum system can be evaluated by the operator space entanglement entropy (OSEE), which we denote as $S^{o}$ in comparison to the unitary case in Eq.(28). Mathematically, OSEE is defined as Prosen2008a

[TABLE]

In Fig. 3(a), we plot the distribution of $S^{o}$ as a funciton of bond index $i$, while in Fig. 3(b) we plot the growth of $S^{o}$ as a function of the evolution depth $d$, for $\bar{n}=0.99,0.6,0.2$ and for bosonic particles. We can see that as the average occupation $\bar{n}$ goes down, $S^{o}$ grows much slower, which is in correspondence with the theoretical results in Ref. Aaronson2016 .

Here we point out that $S^{o}$ is essentially a different concept than $S$ and does not represent the entanglement of a density operator Prosen2008a . It just characterizes the hardness of an density operator written as an MPS. Thus, we cannot study the difference between ideal BosonSampling and lossy BosonSampling by directly comparing the entanglement entropy $S$ and the OSEE $S^{o}$. Moreover, the entropy $S$ and OSEE $S^{o}$ only indicate the difficulty of a problem when simulated with MPS algorithms, which are not necessarily equivalent to the complexity of the underlying problem.

To this end, we note that MPS is an approximate representation of the quantum state, in that each time a two-body gate operation is performed, an SVD is done and the singular values under the threshold $t$ are thrown away. Therefore it is possible that such small errors during the evolution will accumulate after many steps of evolution. The number of two-body gate operations in the circuit is approximately $M^{2}/2$. Thus in the worst case the error $\epsilon$ for the output probability calculated by MPS would be $\epsilon\approx M^{2}t/2$ (number of two-body operations times the truncation threshold for each two-body operation). For the case $M=20$ and $N=10$, the value of $\epsilon$ would be around $2\times 10^{-6}$ at most. To test the overall quality of the MPS method, we generate $10^{3}$ random optical quantum circuits. For each circuit, we calculate $K$ output probabilities according to $K$ randomly chosen basis, where we have used $K=10^{4}$ in our simulation. Each output probability computed with MPS is denoted as ${\rm perm}^{{\rm MPS}}_{j,k}$ in the bosonic case and ${\rm det}^{{\rm MPS}}_{j,k}$ in the fermionic case, where $j$ labels the $j$-th random optical circuit and $k$ labels the $k$-th random basis. In Fig. 4, We compare ${\rm perm}^{{\rm MPS}}_{j,k}$ with the samples from a brute force application of the Ryser’s algorithm, denoted as ${\rm perm}^{{\rm Ryser}}_{j,k}$ in the bosonic case, and compare ${\rm det}^{{\rm MPS}}_{j,k}$ with ${\rm det}_{j,k}$ computed with the python built-in $\det$ function in the fermionic case. Each point on Fig. 4 corresponds to

[TABLE]

We can see that the numerical difference between our results with the standard algorithm reaches smaller than $10^{-9}$, indicating a remarkable accuracy of our approach. We could also increase this precision by a lower $t$, which would also mean that more resources would be required.

IV Conclusion

In summary, we have used a U(1) symmetric MPS to study the dynamical evolution of photons in the optical quantum circuit. We monitor the entanglement entropy growth during the dynamical process, and show that the entropy of a bosonic circuit grows much faster than that of a fermionic circuit. We have also studied lossy BosonSampling, showing that as the loss increases, the operator space entanglement entropy grows slower, which means that classically simulating an optical circuit with larger loss using MPS could be much easier than an ideal one. MPS provides a generic platform to study universal multi-particle quantum walks. More investigations into interesting quantum phenomena in other types of multi-particle quantum walks could be carried out in future works. Other BosonSampling protocols that could significantly reduce the physical resource requirements, and validation methods based on high-order correlation measurements, such as the work in Ref. Walschaers2016 ; Huang2017 , could also be further investigated.

Acknowledgements.

We thank Dario Poletti and Zu-En Su for helpful discussions. C.G. acknowledges support from National Natural Science Foundation of China under Grants No. 11504430 and No. 11805279. H.-L. Huang acknowledges support from the Open Research Fund from State Key Laboratory of High Performance Computing of China (HPCL) (No. 201901-01)