Dynamical phase transitions in sampling complexity

TL;DR

This paper investigates how the complexity of approximately simulating quantum systems via sampling changes over time, revealing a dynamical phase transition in sampling difficulty related to physical phases.

Contribution

It establishes bounds on sampling complexity over time for bosonic systems and identifies a phase transition in classical simulability linked to physical phases.

Findings

Sampling complexity has an upper bound scaling with system size.

A lower bound indicates when sampling becomes computationally hard.

Systems in the Anderson-localized phase are always easy to sample from.

Abstract

We make the case for studying the complexity of approximately simulating (sampling) quantum systems for reasons beyond that of quantum computational supremacy, such as diagnosing phase transitions. We consider the sampling complexity as a function of time $t$ due to evolution generated by spatially local quadratic bosonic Hamiltonians. We obtain an upper bound on the scaling of $t$ with the number of bosons $n$ for which approximate sampling is classically efficient. We also obtain a lower bound on the scaling of $t$ with $n$ for which any instance of the boson sampling problem reduces to this problem and hence implies that the problem is hard, assuming the conjectures of Aaronson and Arkhipov [Proc. 43rd Annu. ACM Symp. Theory Comput. STOC '11]. This establishes a dynamical phase transition in sampling complexity. Further, we show that systems in the Anderson-localized phase are always easy to sample from at arbitrarily long times. We view these results in the light of classifying phases of physical systems based on parameters in the Hamiltonian. In doing so, we combine ideas from mathematical physics and computational complexity to gain insight into the behavior of condensed matter, atomic, molecular and optical systems.

Equations

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}[s]=\frac{1}{r!s!}|\mathrm{Per}(A){|}^2,$$

$$|[a_i(t),a_j^{†}(0)]|=|R_{ij}(t)|\le \min \left(1,\exp \left(\frac{vt-{\ell }_{ij}}{\xi }\right)\right),$$

$$|R_{ij}|\le \exp \left(\frac{-{\ell }_{ij}}{\xi }\right).$$

$$\underset{{\mathcal{D}}_{DP}}{\mathrm{Pr}}[s]=\sum _{\sigma }\frac{1}{s!}{\mathcal{P}}_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{\mathcal{P}}_{{\mathsf{in}}_2,{\mathsf{out}}_{\sigma (2)}}\dots {\mathcal{P}}_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}},$$

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}[s]=\frac{1}{r_1!r_2!\dots r_m!s_1!s_2!\dots s_m!}|⟨\mathrm{vac}|b_1^{s_1}{b}_2^{s_2}\dots b_m^{s_m}{a}_1^{{†}^{r_1}}{a}_2^{{†}^{r_2}}\dots a_m^{{†}^{r_m}}|\mathrm{vac}⟩{|}^2$$

$$=\frac{1}{r!s!}|⟨\mathrm{vac}|({\hat{U}}^{†}{a}_1^{s_1}\hat{U})({\hat{U}}^{†}{a}_2^{s_2}\hat{U})\dots ({\hat{U}}^{†}{a}_m^{s_m}\hat{U})a_1^{{†}^{r_1}}{a}_2^{{†}^{r_2}}\dots a_m^{{†}^{r_m}}|\mathrm{vac}⟩{|}^2$$

$$=\frac{1}{r!s!}|⟨\mathrm{vac}|{\left(\sum _{k_1=1}^m{R}_{1k_1}^{†}{a}_{k_1}\right)}^{s_1}\dots {\left(\sum _{k_m=1}^m{R}_{m{k}_m}^{†}{a}_{k_m}\right)}^{s_m}{a}_1^{{†}^{r_1}}{a}_2^{{†}^{r_2}}\dots a_m^{{†}^{r_m}}|\mathrm{vac}⟩{|}^2,$$

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}[s]=\frac{1}{r!s!}{\left|\sum _{\sigma }\prod _i{R}_{{\mathsf{out}}_{\sigma (i)},{\mathsf{in}}_i}^{†}\right|}^2$$

$$=\frac{1}{r!s!}{\left|\sum _{\sigma }\prod _i{A}_{\sigma (i),i}^{†}\right|}^2,$$

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}[s]=\frac{1}{r!s!}{\left|\mathrm{Per}(A^{†})\right|}^2=\frac{1}{r!s!}|\mathrm{Per}(A){|}^2,$$

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}(s)=|\varphi {|}^2,\mathrm{with}$$

$${\varphi }^{\ast }=\frac{1}{\sqrt{r!s!}}\sum _{\sigma }{R}_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{R}_{{\mathsf{in}}_2,{\mathsf{out}}_{\sigma (2)}}\dots R_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}}$$

$$=\frac{1}{\sqrt{s!}}\sum _{\sigma }{A}_{1,\sigma (1)}{A}_{2,\sigma (2)}\dots A_{n,\sigma (n)},$$

$$\underset{{\mathcal{D}}_U}{\mathrm{Pr}}(s)=\frac{1}{s!}\sum _{\sigma }|R_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{|}^2|R_{{\mathsf{in}}_2,{\mathsf{out}}_{\sigma (2)}}{|}^2\dots |R_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}}{|}^2+$$

$$\frac{1}{s!}\sum _{\sigma \ne \tau }{R}_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{R}_{{\mathsf{in}}_2,{\mathsf{out}}_{\sigma (2)}}\dots R_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}}(R_{{\mathsf{in}}_1,{\mathsf{out}}_{\tau (1)}}{R}_{{\mathsf{in}}_2,{\mathsf{out}}_{\tau (2)}}\dots R_{{\mathsf{in}}_n,{\mathsf{out}}_{\tau (n)}}{)}^{\ast }.$$

$$\underset{{\mathcal{D}}_{DP}}{\mathrm{Pr}}(s)=\frac{1}{s!}\sum _{\sigma }{\mathcal{P}}_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{\mathcal{P}}_{{\mathsf{in}}_2,{\mathsf{out}}_{\sigma (2)}}\dots {\mathcal{P}}_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}},$$

$$\epsilon$$

$$\le \sum _s\frac{1}{2}\sum _{\sigma \ne \tau }|R_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}\dots R_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}}||R_{{\mathsf{out}}_{\tau (1)},{\mathsf{in}}_1}\dots R_{{\mathsf{out}}_{\tau (n)},{\mathsf{in}}_n}|$$

$$=\sum _s\frac{1}{2}\sum _{\sigma ,\rho }|R_{{\mathsf{in}}_1,{\mathsf{out}}_{\sigma (1)}}{R}_{{\mathsf{out}}_{\sigma (1)},{\mathsf{in}}_{\rho (1)}}|\dots |R_{{\mathsf{in}}_n,{\mathsf{out}}_{\sigma (n)}}{R}_{{\mathsf{out}}_{\sigma (n)},{\mathsf{in}}_{\rho (n)}}|,$$

$$\sum _{\sigma }\sum _{\tau }\prod _i|R_{{\mathsf{out}}_{\tau (i)},{\mathsf{in}}_i}|=\sum _{\sigma }\sum _{\tau }\prod _i|R_{{\mathsf{out}}_i,{\mathsf{in}}_{{\tau }^{-1}(i)}}|\underset{i\to \sigma (i)}{\overset{\text{rearrange}}{\to }}\sum _{\sigma }\sum _{\tau }\prod _i|R_{{\mathsf{out}}_{\sigma (i)},{\mathsf{in}}_{{\tau }^{-1}(\sigma (i))}}|$$

$$\epsilon \le \frac{1}{2}\sum _j\sum _{\rho \ne \mathsf{Id}}\prod _i|R_{{\mathsf{in}}_i,j_i}||R_{j_i,{\mathsf{in}}_{\rho (i)}}|,$$

$$\sum _{\begin{array}{c}{j}_i,\rho \\ \rho (i)\ne i\end{array}}|R_{{\mathsf{in}}_i,j_i}||R_{j_i,{\mathsf{in}}_{\rho (i)}}|=C_i\le c\forall i\mathrm{and}$$

$$\sum _{\begin{array}{c}{j}_i,\rho \\ \rho (i)=i\end{array}}|R_{{\mathsf{in}}_i,j_i}||R_{j_i,{\mathsf{in}}_{\rho (i)}}|=\sum _{j_i}|R_{{\mathsf{in}}_i,j_i}{|}^2=D_i=1\forall i.$$

$$\epsilon \le \frac{1}{2}\sum \prod _{i\in {\mathcal{I}}_C}{C}_i\prod _{k\in {\mathcal{I}}_D}{D}_k,$$

$$\epsilon \le \frac{1}{2}\sum _{l=|{\mathcal{I}}_C|=2}^n\left(\genfrac{}{}{0ex}{}{n}{l}\right)c^l=(c+1{)}^n-nc-1.$$

$$\epsilon$$

$$=\left(\genfrac{}{}{0ex}{}{n}{2}\right)(1+h{)}^{n-2}{c}^2\mathrm{for}\mathrm{some}h\in [0,c](\mathrm{by}\mathrm{Taylo}{\mathrm{r}}^{\prime }\mathrm{s}\mathrm{theorem})$$

$$\epsilon$$

$$\le O\left(\exp \left[2\frac{vt-L}{\xi }+2(d-1)\log L\right]\right).$$

$$C_i=\sum _{j_i}\sum _{\begin{array}{c}\rho (i)\ne i\end{array}}\left|R_{{\mathsf{in}}_i,j_i}\right|\left|R_{j_i,{\mathsf{in}}_{\rho (i)}}\right|.$$

$$C_i=\sum _{\begin{array}{c}j\\ \ell \le L\end{array}}\sum _{\begin{array}{c}k\in \mathsf{in}\\ k\ne i\end{array}}\left|R_{i,j}\right|\left|R_{j,k}\right|+\sum _{\begin{array}{c}j\\ \ell >L\end{array}}\sum _{\begin{array}{c}k\in \mathsf{in}\\ k\ne i\end{array}}\left|R_{i,j}\right|\left|R_{j,k}\right|.$$

$$\sum _{\begin{array}{c}j\\ \ell \le L\end{array}}\left|R_{i,j}\right|\sum _{\begin{array}{c}k\in \mathsf{in}\\ k\ne i\end{array}}\left|R_{j,k}\right|$$

$$\le \sqrt{\left(\sum _{\begin{array}{c}j\\ \ell \le L\end{array}}{1}^2\right)\left(\sum _{\begin{array}{c}j\\ \ell \le L\end{array}}{\left|R_{i,j}\right|}^2\right)}{\mathrm{e}}^{vt/\xi }\sum _{k:\Vert {\vec{x}}_k\Vert \ge 1}{\mathrm{e}}^{(-2L\Vert {\vec{x}}_k\Vert +L)/\xi }$$

$$\le a{L}^{d/2}{e}^{(vt+L)/\xi }\sum _{k:||{\vec{x}}_k||\ge 1}{e}^{-2L||{\vec{x}}_k||/\xi }$$

$$\le ab{e}^{(vt-L)/\xi }{L}^{d/2}.$$

$$\sum _{\begin{array}{c}j\\ \ell >L\end{array}}\left|R_{i,j}\right|\sum _{\begin{array}{c}k\in \mathsf{in}\\ k\ne i\end{array}}\left|R_{j,k}\right|$$

$$\le \sum _{\begin{array}{c}j\\ \ell >L\end{array}}\left|R_{i,j}\right|\left(1+{\mathrm{e}}^{vt/\xi }\sum _{k:||{\vec{x}}_k||\ge 1}{\mathrm{e}}^{(L-2L\Vert {\vec{x}}_k\Vert )/\xi }\right)$$

Full text

In the quest towards building scalable and fault-tolerant quantum computers, demonstration of a quantum speedup over the best possible classical computers is an important milestone and is termed quantum computational supremacy Preskill (2012). There are several candidates for tasks where such a speedup could be demonstrated Terhal and DiVincenzo (2002, 2004); Bremner et al. (2010); Aaronson and Arkhipov (2011); Jozsa and den Nest (2013); Fefferman and Umans (2016); Bremner et al. (2016); Farhi and Harrow (2016); Aaronson and Brod (2016); Fefferman et al. (2017); Gao et al. (2017); Harrow and Montanaro (2017), where the problem is to simulate a quantum system in the sense of approximate sampling. However, there has also been some debate about the required system size before one can claim quantum computational supremacy, due to improved simulation techniques and algorithms Neville et al. (2017); Clifford and Clifford (2018); Dalzell et al. (2018). In this Letter, we consider the impact of the field of quantum computational supremacy on other areas of physics and show that studying the complexity of simulating physical systems is useful for understanding phase transitions.

Here we consider the classical complexity of approximate sampling, referred to as “sampling complexity”. This is the task of producing samples from a distribution close to the probability distribution occurring in a quantum system upon measurement in a standard basis. This task is a good notion of what it means to simulate physics on a classical computer since it captures how well a computer can mimic an experiment in which one can measure the output at several sites. When we consider sampling complexity as a function of system parameters, the system can be classified as easy in some regimes and hard in some others. Since the designations “easy” and “hard” are exhaustive and there is no smooth way to go from one regime to another, we posit that this transition from easy to hard happens abruptly, a phenomenon very reminiscent of phase transitions. Just like order parameters are zero on one side of a phase transition and nonzero on another, sampling complexity is different on either side of the transition, and can be used to draw phase boundaries as a function of time and other system parameters. These boundaries can be different from those drawn by more conventional order parameters, signifying new physics in otherwise well-studied systems. Indeed, phase transitions in average-case complexity have been studied both in the classical Kirkpatrick and Selman (1994) and quantum regime Laumann et al. (2010).

In this Letter, we show that transitions in sampling complexity Seshadreesan et al. (2015) are indicative of physical transitions. We consider a system of $n$ bosons hopping from one site to another on a lattice of $m$ sites and study the sampling complexity as a function of time for an initial product state. We show that it goes from easy to hard as the scaling of evolution time $t$ with the number of bosons $n$ increases, thereby exhibiting a dynamical phase transition Heyl et al. (2013); Heyl (2017). We find that the timescale at which the complexity changes is the timescale when interference effects start becoming relevant. We conjecture that in general, this is linked to the Ehrenfest timescale at which quantum effects in a system become considerable Rozenbaum et al. (2017). We also show that systems in the Anderson-localized phase are always easy to simulate. We use the Lieb-Robinson bound Lieb and Robinson (1972) as an ingredient in our proof of easiness.

Setup.— The model consists of free bosons hopping on a lattice in $d$ dimensions (denoted $d$-D) with sites labeled by indices $i,j$. Our results in this paper can be applied to linear optics as well, with the bosonic sites being replaced by photonic modes. The Hamiltonian is given by $H=\sum_{i,j=1}^{m}J_{ij}(t)a_{i}^{\dagger}a_{j}$, where $a_{i}^{\dagger}$ is the creation operator of a boson at the $i$’th site. $J(t)$, which can be time-dependent in general, is an $m\times m$ Hermitian matrix that encodes the connectivity of the lattice. One way to show hardness in the boson sampling proposal by Aaronson and Arkhipov (AA) Aaronson and Arkhipov (2011) is to generate any linear optical unitary $U$ acting on the bosonic sites, in particular, Haar-random unitaries. Any $U$ can be generated through a free boson Hamiltonian by taking $H=i\log U$ and evolving it for unit time. However, this Hamiltonian can require arbitrarily long-range hops on the lattice in general. Since such long-ranged Hamiltonians may not be realistic, we consider the sampling complexity of a bosonic Hamiltonian with nearest-neighbor hops. $J_{ij}$ is nonzero only if $i=j$ or $i$ and $j$ label adjacent sites. We further restrict $J_{ij}$ to satisfy $|J_{ij}|\leq 1$ in order to set an energy scale. Our proofs remain valid even if we allow the diagonal terms $J_{ii}$ to be unbounded.

One can efficiently solve the equations of motion $i\dot{a_{i}}^{\dagger}(t)=[a_{i}^{\dagger}(t),H(t)]$ on a classical computer to obtain $a_{i}^{\dagger}(t)=\sum_{k}a_{k}^{\dagger}(0)R_{ki}(t)$ for some transformation matrix $R$. From here onward, we shall take $R(t)$ to be the input to the problem, since it can be determined from the input Hamiltonian $H$ and time $t$ in time $\mathrm{poly}(m,\log t)$.

The $m$ sites in the problem are numbered from 1 to $m$, and together with $n$ ancilla sites, are arranged in a lattice of side length $(m+n)^{1/d}$ in $d$ dimensions. The initial state has $n$ bosons equally spaced in the lattice as shown in Fig. 1. We take $m=c_{1}n^{\beta}$, where $\beta$ controls the sparsity of occupied sites in the lattice and can be set to 5 as required for the hardness of boson sampling Aaronson and Arkhipov (2011). The minimum spacing between any two bosons in the initial state is $2L=\left(\frac{m+n}{n}\right)^{1/d}>c_{1}^{1/d}n^{\frac{\beta-1}{d}}$. The quantity $L$ is an important length scale in the problem. The ancillas in the lattice, marked in cyan, are not counted as part of the $m$ sites and are present in order to accelerate the time required to construct an arbitrary unitary, which is useful for the hardness result. Their presence does not change the scaling of quantities like $L$ with $n$.

The input states are described by vectors of the form $r=(r_{1},\ldots,r_{m})$, specifying the number of bosons on each site, so that $r_{1}+\ldots r_{m}=n$. Measurement in the boson number basis defines a distribution $\mathcal{D}_{U}$, which we aim to sample from. The probability of finding an output state $s=(s_{1},s_{2},\ldots,s_{m})$ is given by

[TABLE]

where $r!:=r_{1}!\ldots r_{m}!$ (with $s!$ defined similarly), $A_{n\times n}$ is a matrix formed by taking $s_{i}$ copies of the $i$’th column and $r_{j}$ copies of the $j$’th row of $R$ in any order, and $\mathrm{Per}(A)$ denotes the permanent of $A$ (see Refs. Aaronson (2011); SM_ for details).

For the particular choice of initial states described in Fig. 1, the task is to sample from a distribution that is close to $\mathcal{D}_{U}$ in variation distance when given a description of the unitary $R(t)$. We now formalize the notion of efficient sampling.

2. $f(n)=\Omega(g(n))\Leftrightarrow g(n)=O(f(n))$.

3. $f(n)=\Theta(g(n))$ if $f(n)=\Omega(g(n))$ and $f(n)=O(g(n))$.).*

We call the sampling problem easy if there exists an efficient sampler for the problem in the stated regime. Conversely, the problem is hard if there cannot be an efficient sampler. Since a negative statement such as the inexistence of an algorithm is difficult to prove, our practical definition of hardness of a sampling problem is if it is at least as hard as boson sampling. This enables us to use the results of AA Aaronson and Arkhipov (2011) to claim the hardness of sampling in some regimes. In doing so, our hardness results ultimately rely on the truth of AA’s conjectures. One of these conjectures concerns the hardness of additively approximating $|\mathrm{Per}(G)|^{2}$ for Gaussian-random $G$ with high probability (for which AA give reasonable evidence). The second, more widely believed conjecture is that the polynomial hierarchy, an infinite tower of complexity classes, does not “collapse”, i.e. is truly infinite.

We restrict our attention to two special cases where we can show the existence of an efficient sampling algorithm: i) when the system evolves for a time smaller than the system timescale $L/v$ (where $v$ is the Lieb-Robinson velocity of information spreading in the lattice, defined more precisely in Eq. (2)), and ii) when there is Anderson localization in the system Anderson (1958). These two cases correspond to a promise on the input unitary $R$. We now state our main results.

The intuition behind this Theorem is that when the time is smaller than the Lieb-Robinson timescale of particle interference $L/v$, the dynamics is approximately classical (in the sense that the particles are distinguishable).

An important technical achievement of this paper is to prove rigorously that this intuition is correct. This is done in Lemma 3 by showing that the approximation works.

For this result, we show that we can apply any unitary $R$ after the stated time. Therefore, if an efficient sampler exists for this problem, then we have an efficient boson sampler, which cannot exist by our assumption.

The result for the case of Anderson localization comes out as a corollary from Theorem 1.A:

The easiness of sampling for Anderson-localized systems is analogous to results showing efficient simulation of localized systems according to various definitions Friesdorf et al. (2015); Huang (2015); Pollmann et al. (2016); Yu et al. (2017); Chapman and Miyake (2017).

Easiness at short times.— In this section, we prove Theorem 1.A and Corollary 2. First, let us examine the promise we have on the unitary $R$ in both cases. We use the Lieb-Robinson bound Lieb and Robinson (1972) on the speed of information propagation in a system. Applying the bound to our Hamiltonian, we get

[TABLE]

where $\ell_{ij}$ is the distance between two sites $i$ and $j$, $v$ is the upper bound to the velocity of information propagation called the Lieb-Robinson velocity and $\xi$ is a length scale). Note that the results from Ref. Eisert and Gross (2009) do not apply since we have free bosons here and we work in the single-particle subspace. The Lieb-Robinson velocity is at most $4(1+2de)$ Hastings (2010) when $|J_{ij}|\leq 1$ and $\xi=1$.

When the Hamiltonian is Anderson-localized, the unitary $R$ satisfies the following promise at all times Fröhlich et al. (1985):

[TABLE]

Here, $\xi$ is the maximum localization length among all eigenvectors. Eq. (3) can be viewed as a consequence of a Lieb-Robinson bound with zero velocity Hamza et al. (2012). On account of the zero-velocity Lieb-Robinson bound, all results for the time-dependent case can be ported to the Anderson localized case, setting $v=0$.

We give an algorithm that efficiently samples from the output distribution for short times $t<\frac{0.9L}{v}=O(n^{(\beta-1)/d})$, given the promise in Eq. (2). The algorithm outputs a sample from a distribution $\mathcal{D}_{DP}$, obtained by assuming that the bosons are distinguishable particles, ignoring the effects of interference. The algorithm is described in more detail in Ref. SM_ .

Analysis.— We prove the correctness of the algorithm by showing that the variation distance between the distributions is upper bounded by an inverse exponential in $n$. When the bosons are distinguishable, their dynamics is given by a Markov process, described by the matrix $\mathcal{P}_{kl}=|R(t)|_{kl}^{2}$. The probability of getting an outcome $s$ is given by

[TABLE]

where the sum is over all permutations $\sigma$ mapping the input bosons to the output ones. In the above equations, $\mathsf{in}_{i}$ is the site index of the $i$’th boson in nondecreasing order in the input and $\mathsf{out}$ is defined similarly at the output (see Ref. SM_ for an example). We now state a result on how close $\mathcal{D}_{DP}$ is to the true distribution $\mathcal{D}_{U}$.

This Lemma makes intuitive sense: because of the Lieb-Robinson bound Eq. (2) and the fact that the initially occupied sites are separated by a minimum distance $\Theta(L)$ from each other, it takes a time $t=\Theta(L/v)$ for the bosons to start interfering considerably. Therefore, the classical and quantum distributions agree exponentially closely in $L$ when $t\leq 0.9L/v$. For a proof of Lemma 3, see Ref. SM_ . Assuming this Lemma, we now show Theorem 1.A.

Hardness at longer times.— If we allow the system to evolve for a longer amount of time, we can use the time-dependent control to effect any arbitrary unitary and implement any boson sampling instance in the system. We can perform phase gates on a site $k$ by setting $J_{kk}$ to be nonzero for a particular time, with the hopping terms and other diagonal terms set to zero. We can apply a nontrivial two-site gate between adjacent sites, for example the balanced beamsplitter unitary on the sites 1 and 2, $U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\ 1&1\end{pmatrix}$, by setting $H=-i(a_{1}^{\dagger}a_{2}-a_{1}a_{2}^{\dagger})$ for time $t=\frac{\pi}{4}$. One can also apply arbitrary unitaries using arbitrary on-site control $J_{ii}(t)$ and fixed, time-independent nearest-neighbor hopping $J_{ij}(t)=1$.

Using the constructions in Refs. Reck et al. (1994); Clements et al. (2016), we can effect any arbitrary $m\times m$ unitary on the $m$ sites with $O(m^{2})$ depth from a nontrivial beamsplitter and arbitrary single-site phase gates. A construction of AA that employs ancillas to obtain the desired final state rather than applying the full unitary on all the sites can be used to achieve a depth $O(nm^{1/d})$. Each of the $n$ columns of the unitary are implemented in time $O(m^{1/d})$, which corresponds to the timescale set by the Lieb-Robinson velocity and the distance between the furthest two sites in the system (see note 333AA’s construction applied each column of the unitary in $O(\log m)$-depth, whereas we can only apply it in $O(m^{1/d})$ depth because of the spatial locality of the Hamiltonian and the Lieb-Robinson bounds that follow.).

Outlook.— We have defined the sampling problem for local Hamiltonian dynamics and given upper and lower bounds for the scaling of time $t(n)$ with the number of bosons $n$ for which the problem is efficiently simulable or hard to classically simulate, respectively. Our results are captured in Fig. 2 that illustrates the complexity phase diagram of the system. For time-independent systems, we observe that boson sampling is classically easy for all times if the system is Anderson-localized. In a future work Deshpande et al. , we show that there is a class of static, local Hamiltonians that generate a hard-to-sample output distribution at some time that is not formally infinity. This means that sampling complexity distinguishes Anderson-localized and delocalized systems, which makes it similar to an order parameter that distinguishes different phases.

The case with nonzero Lieb-Robinson velocity [Eq. (2)] shows two regimes of the scaling of $t$ with $n$ where sampling is provably easy/hard. We have shown that sampling is easy when $t\leq t_{\mathrm{easy}}=\Theta(n^{\frac{\beta-1}{d}})$ and hard when $t\geq t_{\mathrm{hard}}=\Theta(n^{1+\frac{\beta}{d}})$. Since our definitions of easiness and hardness are exhaustive, we argue that there must exist a constant $c$ such that sampling is efficient for $t<\Theta(n^{c})$ and hard otherwise, illustrating a dynamical phase transition. Our proof implies that $c\in[\frac{\beta-1}{d},\frac{\beta+d}{d}]$ and we show in future work Deshpande et al. that the transition is sharp (see note 444The transition here is sharp if there is a specific time $t_{\mathrm{samp}}=\Theta(n^{c})$ such that sampling is easy for all $t<t_{\mathrm{samp}}$ and hard otherwise. The transition is coarse if sampling is easy for all $t=\Theta(n^{c})$ and hard for all $t$ such that $\lim_{n\rightarrow\infty}\frac{t}{n^{c}}\rightarrow\infty$.). The transition is between two regimes, one for short times in which the system’s dynamics is essentially indistinguishable from classical dynamics; and the other in which quantum mechanical effects dominate to such an extent as to forbid an efficient classical simulation.

This result may be viewed as a generalization of a similar result in Ref. Brod (2015), where it was shown that exact boson sampling for depth-4 circuits is hard. The results there are not directly comparable to ours, since Ref. Brod (2015) assumes $\beta=1$, whereas our results need $\beta>1$ (easiness) and $\beta\geq 2\ \mathrm{or}\ 5$ (hardness). The reason we get easiness even after polynomial time is that we deal with approximate sampling, a less stringent notion of simulation. In addition, adapting the hardness proof of Ref. Brod (2015) into our setup, we obtain that the system is hard to exactly sample from at timescale $\Omega(8L/v)=\Omega(n^{\frac{\beta-1}{d}})$ Deshpande et al. , partially answering what happens in the subregion with the question mark in Fig. 2.

Recent studies Rozenbaum et al. (2017); Banerjee and Altman (2017); He and Lu (2017); Syzranov et al. (2017) have also studied transitions based on time-dependent features of the out-of-time-ordered correlator. This raises the question of the connection between sampling complexity and scrambling time Hosur et al. (2016); *Roberts2017 in quantum many-body systems Huang et al. (2016); *Fan2017; *Chen2016a; *Swingle2017 and fast scramblers like the Sachdev-Ye-Kitaev Sachdev and Ye (1993); Kitaev (2015); Maldacena and Stanford (2016); Maldacena et al. (2016); Banerjee and Altman (2017) models, and black holes Sekino and Susskind (2008); Lashkari et al. (2013); Swingle (2016); Brown et al. (2016a); *Brown2016a, where one can explore the connection to recent conjectures on complexity in the dual CFT Aaronson and Chen (2016); Brown et al. (2016a); *Brown2016a; Brown and Susskind (2018). Further, it would be interesting to study sampling complexity in various other physical settings Marzolino and Prosen (2016), like systems with topological or many-body-localized phases and systems in their ground state to explore the connection with Hamiltonian complexity.

We therefore propose another motivation for considering sampling complexity beyond the field of quantum computational supremacy: complexity provides a natural way to classify phases of matter that is complementary to traditional approaches based on symmetries and topology. This is akin to how the study of entanglement in many-body physics has helped us understand phases of matter Amico et al. (2008); Zeng et al. (2015) and characterize thermalization and localization Bardarson et al. (2012); Yang et al. (2017).

Coming to experimental implementations, platforms such as ultracold atoms in optical lattices Morsch and Oberthaler (2006); Bloch et al. (2012) and superconducting circuits Buluta and Nori (2009); Devoret and Schoelkopf (2013) are ideal for experimentally studying the transition by comparing the distribution sampled by the algorithm and the experimental distribution. The ingredients required, which have been realized in several groups Bakr et al. (2009, 2010); Sherson et al. (2010); Mazza et al. (2015); Miranda et al. (2015); Yamamoto et al. (2017), are: 1. Preparation of the initial state Weitenberg et al. (2011); Wang et al. (2015) of the type shown in Fig. 1 (see note 555All of our arguments go through even if the initial state is not regular as in Fig. 1 but still has the property that the initially occupied sites are separated by a minimum distance of $\Theta(L)$ from each other.). 2. Evolution under a Hamiltonian with either arbitrary time-dependent nearest-neighbor hopping strength or fixed nearest-neighbor hopping strength $J_{ij}(t)=1$ together with arbitrary time-dependent on-site potential $J_{ii}(t)$ Weitenberg et al. (2011); Fukuhara et al. (2013), and 3. Single-site resolved measurement of occupation number of the sites Bakr et al. (2010); Sherson et al. (2010). Cold atoms in quantum gas microscopes can be controlled at the single site level, enabling all three ingredients above. To maintain integrability, we can turn off the Hubbard interaction for the bosonic atoms by tuning to a Feshbach resonance Paul et al. (2016); Xu et al. (2003); Herbig (2003); Dürr et al. (2004).

New architectures like optical tweezers Kaufman et al. (2014) are also promising since they allow for deterministic creation of desired initial states and feature tunable interactions in time Endres et al. (2016); Bernien et al. (2017). Similarly, superconducting circuits have been proposed for quantum simulation of quantum walks Flurin et al. (2017) and the Bose-Hubbard model Hacohen-Gourgy et al. (2015); Deng et al. (2016); Roushan et al. (2017). The gmon qubit architecture, which was used in Ref. Roushan et al. (2017), naturally allows for time-dependent variation of coupling strengths Chen et al. (2014), and raises the prospect of an experimental study of the transition.

S1 Supplemental material

In this supplemental material, we give expressions for the output probabilities in the distribution $\mathcal{D}_{U}$ in a boson sampling experiment. We then explicitly present the algorithm and derive the expression for $\mathcal{D}_{DP}$. We then derive an upper bound to the variation distance between them, proving Lemma 3 of the main text.

Expression for output probabilities.— In this section, we describe the standard boson sampling set-up and derive an expression for the output probabilities of a boson sampling experiment that define the distribution $\mathcal{D}_{U}$. First, let us represent the input and output states pictorially and develop some notation.

In Fig. S1, the top line denotes the input state $|r\rangle$ and the bottom line the output $|s\rangle$. Each filled circle denotes a boson occupying the corresponding lattice site, which is labeled below the circles. The numbers marked in orange above each boson label the bosons from left to right (more generally, this is in a nondecreasing order of the site index). We will call this the boson index.

A given configuration (basis state) is completely specified by specifying the boson number in each site, such as $r=(1,0,0,0,1,0,0,0,1)$ and $s=(0,1,0,0,1,1,0,0,0)$ in the above. It can also be specified by listing the site index for every boson index, i.e. the occupied sites. Thus the input state can be represented as $\mathsf{in}=(1,5,9)$, the output state as $\mathsf{out}=(2,5,6)$.

All $n!$ permutations of the boson indices represent valid paths that the bosons can take to the output state, and correspond to the $n!$ terms in the permanent of the matrix. In cases where there are two or more bosons in a particular site at the input or output, there are $\frac{n!}{r!s!}$ paths (and terms in the amplitude). Here, $r!:=r_{1}!r_{2}!\ldots r_{m}!$ and similarly $s!$. By taking repeated rows and columns of $R$, this has the effect of still giving $n!$ terms in total, which we identify with the $n!$ permutations in the boson indices. The expression for the probability of an outcome $s$ is (here, $b_{i}:=a_{i}(t)$):

[TABLE]

where $R_{ij}$ describes the action of $\hat{U}$ on the annihilation operators at a site: $b_{i}=a_{i}(t)=\sum_{k}R^{\dagger}_{ik}(t)a_{k}(0)$. Now define the matrix $A^{\dagger}$ to be the one obtained by taking $s_{i}$ copies of the $i$’th row and $r_{j}$ copies of the $j$’th column of $R^{\dagger}$. For concreteness, this can be done by first considering the rows and repeating a row $i$ of $R^{\dagger}$ whenever $s_{i}>1$, or not picking it if $s_{i}=0$, to convert it into an $n\times m$ matrix. We can then do the same with columns to convert it into an $n\times n$ matrix. However, the order ultimately does not matter since the quantity that emerges, the permanent, is symmetric under exchange of rows or columns. We have

[TABLE]

where the sum is over all permutations $\sigma$. This finally gives us

[TABLE]

where $\mathrm{Per}(A)$ is the permanent of $A$.

Algorithm.— The sampling algorithm is given below. It is easy to see that it implements one step of a Markov process of $n$ distinguishable bosons walking on a lattice.

Note that $\mathcal{P}$ from line 1 is a doubly stochastic matrix describing the classical Markov process. To see that the runtime is polynomial in $n$, note that the loop is over $n$ boson indices. Line 1 takes time $O(m\log m)=\tilde{O}(n^{\beta})$, giving a total runtime of $\tilde{O}(nm)=\tilde{O}(n^{1+\beta})$. The notation $\tilde{O}$ suppresses factors of $\log n$.

Bound on variation distance.— Here we derive a bound on the variation distance $\|\mathcal{D}_{U}-\mathcal{D}_{DP}\|=\frac{1}{2}\sum_{s}|\Pr_{\mathcal{D}_{U}}(s)-\Pr_{\mathcal{D}_{DP}}(s)|$. Rewriting the actual probability in terms of the amplitudes, we have

[TABLE]

where $A$ is the $n\times n$ matrix formed by taking the appropriate number of copies of each row and column of $R_{m\times m}$. We have set $r!=1$ since our input state has bosons in distinct sites. Continuing,

[TABLE]

The probability distribution $\mathcal{D}_{DP}$ that the algorithm samples from is given by the first line of Eq. (S10):

[TABLE]

where the sum is over all the $n!$ ways of assigning the $n$ input states to the $n$ output states. As before, the $s!$ is to account for overcounting when two distinct permutations in the boson index refer to the same site index in the output state.

The expression for the probability is proportional to the permanent of the matrix with the positive entries $\mathcal{P}_{\mathsf{in}_{i},\mathsf{out}_{j}}$, and can hence be efficiently approximated Jerrum et al. (2004). Note that the algorithm does not explicitly calculate this probability but only samples from the distribution. We can now prove Lemma 3 of the main text.

Proof of Lemma 3.

The variation distance is given by

[TABLE]

where $\rho=\tau^{-1}\circ\sigma\neq\mathsf{Id}$, the identity permutation. The last equality comes from rearranging the terms in $|R_{\mathsf{out}_{\tau(1)},\mathsf{in}_{1}}\ldots R_{\mathsf{out}_{\tau(n)},\mathsf{in}_{n}}|$ so that the terms involving $R_{\mathsf{in}_{i},\mathsf{out}_{\sigma(i)}}$ and $R_{\mathsf{out}_{\sigma(i)},j}$ (for some $j$) are together:

[TABLE]

Summing over all outcomes $s$ (or configurations $\mathsf{out}$), Eq. (S14) is equivalent to

[TABLE]

where the sum $j$ is over ordered tuples $(j_{1},\ldots j_{n})$, representing the intermediate lattice sites that the bosons in positions $(\mathsf{in}_{1},\ldots\mathsf{in}_{n})$ jump to, before jumping back to positions $(\mathsf{in}_{\rho(1)},\ldots\mathsf{in}_{\rho(n)})$.

We can proceed to break the sum in Eq. (S16) based on the number of fixed points of the permutation $\rho$, that is, the number of indices $i$ such that $\rho(i)=i$. We bound these quantities separately as follows:

[TABLE]

The variation distance is therefore bounded above:

[TABLE]

where the sum is over subsets $\mathcal{I}_{C}$ of the indices representing the input state, $\mathsf{in}$. $\mathcal{I}_{D}$ is the complement of $\mathcal{I}_{C}$ and $|\mathcal{I}_{D}|$ is the number of fixed points. Suppose we find an upper bound $c$ for $C_{i}$ in Eq. (S17), we then have

[TABLE]

In Lemma 4, we show that $c=\eta L^{d-1}\mathrm{e}^{(vt-L)/\xi}$ for some constant $\eta$. Continuing from Eq. (S19),

[TABLE]

Now, plugging in the value of $c$ and assuming that $vt\leq 0.9L$ and $\beta>1$, we get

[TABLE]

In the first line, we use the inequality $\log(1+x)\leq x$ and in the second, $\exp\left[(n-2)\times\eta L^{d-1}\mathrm{e}^{(vt-L)/\xi}\right]=O(1)$ since $vt<L$ and $|vt-L|=\Omega(n^{\beta-1})$. This completes the proof of Lemma 3. ∎

Lemma 4**.**

For all constant dimensions $d$, $c=\eta L^{d-1}\mathrm{e}^{(vt-L)/\xi}$.

Proof.

Recall that

[TABLE]

Since we are looking for a bound that applies for all $i$, let us, for convenience, make the following changes in notation: $\mathsf{in}_{i}\rightarrow i,j_{i}\rightarrow j,\mathsf{in}_{\rho(i)}\rightarrow k$, denoting a boson starting at position $i$, jumping to $j$ and then to $k$, where $i$ and $k\neq i$ are site indices belonging to $\mathsf{in}$. We split the sum in $C_{i}$ based on the distance between $i$ and $j$, $\ell_{ij}=:\ell$.

[TABLE]

Consider the first term:

[TABLE]

Here in the first line, we have used the Lieb-Robinson bound Eq. (2). In the second line, we use the Cauchy-Schwarz inequality and the fact that $\ell=\ell_{ij}\leq L$. In the second line, $\vec{x}_{k}$ is the position vector of site $k$ relative to site $i$, re-scaled by $2L$. Therefore the sum over $\vec{x}_{k}$ is over all vectors with integer coordinates. In the last line, we use Lemma 6, to be proven below. $a$ and $b$ are constants independent of $n$ that depend on the dimension $d$ and the length scale $\xi$.

Now, in the second term for $C_{i}$ in Eq. (S26), the intermediate site $j$ is not necessarily close to $i$. Therefore, there are terms where $j$ is close to $k\neq i$ and one has to treat these terms carefully. For these terms, we use the trivial Lieb-Robinson bound of 1 in Eq. (2) rather than $\exp\left(\frac{vt-\ell_{jk}}{\xi}\right)>1$.

[TABLE]

In the second line, we use 1 as a Lieb-Robinson bound for $|R_{j,k}|$ when $k=k^{*}$, the site belonging to $\mathsf{in}$ that is closest to $j$. All other $k$’s have distances from $j$ bounded below by $2L\|\vec{x}_{k}\|-L$, where $\vec{x}_{k}$ is now the re-scaled position vector of a site $k$ with the origin at $k^{*}$. We apply Lemma 6 in the third line and Lemma 5 in the fifth. Collecting everything, we have

[TABLE]

∎

Lemma 5** ($d$-dimensional sum).**

The sum $\sum_{\|\vec{x}\|\geq L}\mathrm{e}^{-\|\vec{x}\|/\xi}$ over points $\vec{x}$ with integer coordinates is upper bounded by $a_{d}\mathrm{e}^{-L/\xi}\left(\xi L^{d-1}\right)$ for large enough $\frac{L}{\xi}$ for some dimension-dependent constant $a_{d}$.

Proof.

We can view the sum over a lattice of vectors with integer coordinates as a Riemann sum and bound the corresponding $d$-dimensional integral. Consider the quantity

[TABLE]

where $\Gamma(d,x)=\int_{x}^{\infty}y^{d-1}\mathrm{e}^{-y}\mathrm{d}y$ is the incomplete $\Gamma$ function. We can lower bound the integral by the Riemann sum $\sum_{\Delta}V_{\Delta}\mathrm{e}^{-\|\vec{y}\|/\xi}$, where the sum is over cells $\Delta$ with volume $V_{\Delta}$ centered at lattice points $\vec{x}$. $\vec{y}$ is the point in the cell $\Delta$ with the highest norm $\|\vec{y}\|$. Further, the point with the highest norm is not too distant from the one at the center: $\|\vec{y}\|\leq\|\vec{x}\|+\frac{\sqrt{d}}{2}$. Therefore, we have

[TABLE]

We now need an upper bound on the incomplete $\Gamma$ function $\Gamma(d,x)$ for large $x$ Olver et al. :

[TABLE]

Combining Eq. (S38) and Eq. (S39), we get:

[TABLE]

∎

Using a similar method, we also get bounds on a related sum.

Lemma 6**.**

For $\vec{x}\in\mathbb{Z}^{d}$, $f_{d}:=\sum_{\|\vec{x}\|\geq 1}\mathrm{e}^{-2L\|\vec{x}\|/\xi}\leq b_{d}\exp[{-\frac{2L}{\xi}}]$ for some dimension-dependent constant $b_{d}$.

Proof.

We prove the statement by induction on the dimension $d$. For $d=1$, the statement is seen to be true since the sum evaluates exactly:

[TABLE]

For the inductive step, consider the integral $g(d)=\int_{\|\vec{y}\|\geq 1}^{\infty}\mathrm{e}^{-2L\|\vec{y}\|/\xi}\mathrm{d}^{d}\vec{y}$. This is lower bounded by the Riemann sum represented in Fig. S2. The white dots represent vectors with at least one zero coordinate and do not enter the Riemann sum according to this way of dividing the region of integration into cells. In the following, the set of points with at least one zero coordinate is denoted $\O_{cc}$. We have:

[TABLE]

where $\Delta_{\vec{x}}$ is the volume of the cell associated with the lattice vector $\vec{x}$. In Fig. S2, the volume of most cells (whose center is at distance 1.5 or beyond from the origin) is 1. The cells near the origin have some volume $\alpha_{d}<1$ that depends on the dimension. Lower bounding all volumes $\Delta_{\vec{x}}$ by $\alpha_{d}$,

[TABLE]

Now it remains to upper bound the contribution from summing over the points $\O_{cc}$. Notice that the sum over these points is upper bounded by the sum over $d$ hyperplanes of dimenson $d-1$. From the inductive hypothesis,

[TABLE]

since there are $d$ hyperplanes of dimension $d-1$. Adding Eqs. (S45) and (S46), we get

[TABLE]

proving the lemma. In the second line we have expanded the incomplete $\Gamma$ function for large $L/\xi$. ∎