Efficiently simulating the work distribution of multiple identical bosons with boson sampling

Abstract

Boson sampling has been theoretically proposed and experimentally demonstrated to show quantum computational advantages. However, it still lacks the deep understanding of the practical applications of boson sampling. Here we propose that boson sampling can be used to efficiently simulate the work distribution of multiple identical bosons. We link the work distribution to boson sampling and numerically calculate the transition amplitude matrix between the single-boson eigenstates in a one-dimensional quantum piston system, and then map the matrix to a linear optical network of boson sampling. The work distribution can be efficiently simulated by the output probabilities of boson sampling using the method of the grouped probability estimation. The scheme requires at most a polynomial number of the samples and the optical elements. Our work opens up a new path towards the calculation of complex quantum work distribution using only photons and linear optics.

Graphical abstract

Keywords

quantum simulation / quantum work distribution / boson sampling / linear optics

1 Introduction

Work in nonequilibrium systems is a fundamental research topic and has stimulated many research interests in statistical physics [1-3]. Quantum work distribution is a key quantity in the thermodynamic analysis of any quantum system, and determines many important thermodynamic properties, such as the free energy difference [1] and nonequilibrium work relation [2]. Quantum work distribution in a thermally isolated system can be effectively determined by the beginning-time and end-time energy measurements [4]. Previously, there are many explorations on the work distribution in the nonequilibrium quantum system, both theoretically [5-9] and experimentally [10-12]. However, these results mainly focused on the single-particle systems.

In recent years, multi-particle work distribution in nonequilibrium processes has received more and more attentions [13-16]. The calculation of work distribution for an identical multi-particle system involves the transition probability between multi-particle eigenstates, which may be formidable difficulty due to the interference influence of these particles [17-19]. The transition probability between the eigenstates of multiple identical bosons (fermions) associates with a permanent (determinant) of the corresponding transition amplitude matrix between the single-boson (single-fermion) eigenstate. The determinant can be efficiently calculated on a classical computer, while evaluating the permanent is a so-called “#P-complete” hard [20,21]. This implies that the calculation of work distribution with multiple bosons may be classically difficult.

Quantum boson sampling, a remarkably quantum computational supremacy candidate [22], was proposed by Aaronson and Arkhipov in 2011 [23]. Boson sampling emerges as a powerful paradigm to efficiently solve the output probability distribution of photons in a linear optical network, and provides several practical applications in graph theory [24-26], decision and function problems [27,28], quantum chemistry [29-33], random number generation [34], and image encryption [35]. The early proof-of-principle demonstration of boson sampling has experimentally confirmed that the sampling result is related to the matrix permanent [36-41]. Recently, boson samplings with tens [42,43] or even more than a hundred of photons [44] have been experimentally realized, which fully shows a quantum computational advantage over the classical computer. Besides, the scalable implementations of boson sampling utilizing photonic modes of trapped ions have also been proposed [45,46]. Note that when the boson number is N\le 30 and mode number is \mathcal{O}(N^2), certainly a classical computer could simulate the output probability of boson sampling. But when N achieves approximately 50 or larger, and mode number is approximately \mathcal{O}(N^2), it would completely beyond the capability of the classical computers [23].

Though the output probability of boson sampling is associated with the permanent of a matrix, it is completely different from the problem of predicting or estimating the permanent via boson sampling. In fact, it is infeasible to directly estimate the individual output probability using boson sampling, as the detected probability is exponentially small and one has to collect exponentially many samples to achieve a reasonable accuracy [23]. This is one of the main obstacles for limiting the practical applications of boson sampling. Fortunately, it is found that if these output probabilities are grouped and their sums are estimated, the polynomial samples rather than exponential ones are required for solving the certain problems [29-32,47]. In this way, boson sampling could be a potentially effective method to solve some practical problems.

In this paper, we investigate how to use a boson sampling system to simulate the work distribution of multiple non-interaction identical bosons in a one-dimensional quantum piston. We firstly present a general theory for the work distribution of multiple bosons and establish a connection between the work distribution and boson sampling. Specifically, we numerically calculate the transition amplitude matrix of the single-boson eigenstates and program it into an optical network of boson sampling. The work distribution can be efficiently simulated by the output probabilities of boson sampling using the method of the grouped probability estimation (GPE) instead of the individual probability estimation (IPE). We also analyze the effect of system parameters on the work distribution and finally present a feasibility analysis of the scheme in terms of the resource cost. Our proposed way opens up a new possibility for studying the work distribution problem of quantum thermodynamics via the linear quantum optical network of boson sampling.

2 Results

2.1 Work distribution with multiple bosons

We consider that N identical bosons in a quantum system are driven by a varied external work parameter \lambda from initial time t=0 to final time t=\tau. The work parameter \lambda could be the position of the quantum piston or the spring coefficient of the harmonic oscillator. Suppose that at initial time t=0, the parameter is \lambda (0) ={\lambda }_0 and the system is prepared in a thermal state with a heat bath at an inverse temperature \beta =1/(k_{\mathrm{B}}T). Here k_{\mathrm{B}} is the Boltzmann constant and T is the temperature of the system. After preparing a thermal state, we perform a projective measurement over the energy eigenstates of the system. The probability of bosons initially occupying in the energy eigenstate |i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩ is expressed by

(1)P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩)={\displaystyle \frac{1}{Z^{{\lambda }_0}}}\text{exp}(-\beta n_i^{{\lambda }_0}{E}_i^{{\lambda }_0}).

Here Z^{{\lambda }_0} is a partition function, which is given by Z^{{\lambda }_0}= \sum _i\text{exp}(-\beta n_i^{{\lambda }_0}{E}_i^{{\lambda }_0}). E_i^{{\lambda }_0} is the ith eigenenergy at initial time t=0, and its corresponding eigenstate is | i^{{\lambda }_0}⟩. n_i^{{\lambda }_0} represents the initial occupation number of bosons in the ith eigenstate |i^{{\lambda }_0}⟩.

Then the system is detached from the heat bath, and the work parameter of the system is changed from \lambda (0)={\lambda }_0 at initial time t=0 to \lambda (\tau )={\lambda }_{\tau } at final time t=\tau to perform the work on the system. Finally, we apply the second measurement to project the system into its energy eigenstates again. During the process, the system is evolved under the unitary dynamics and the work distribution of the system can be written as [48,49]

(2)\begin{array}{rl}\rho (W)=& \sum _i\sum _{f}P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩)\\ & \times P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩\to |f^{{\lambda }_{\tau }}:n_f^{{\lambda }_{\tau }}⟩)\\ & \times \delta (W-n_f^{{\lambda }_{\tau }}{E}_f^{{\lambda }_{\tau }}+n_i^{{\lambda }_0}{E}_i^{{\lambda }_0}).\end{array}

Here E_f^{{\lambda }_{\tau }} is the fth eigenenergy at final time t=\tau, and its corresponding eigenstate is | f^{{\lambda }_{\tau }}⟩. n_f^{{\lambda }_{\tau }} represents the final occupation number of bosons in the fth eigenstate |f^{{\lambda }_{\tau }}⟩. P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩\to |f^{{\lambda }_{\tau }}:n_f^{{\lambda }_{\tau }}⟩) is a transition probability from the initial multi-boson eigenstate |i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩ to the final multi-boson eigenstate |f^{{\lambda }_{\tau }}:n_f^{{\lambda }_{\tau }}⟩. The work distribution of multiple particles can be further rewritten as

(3)\rho (W)={\displaystyle \sum _{\mathit{F}}P(\mathit{I})\cdot P(\mathit{F}|\mathit{I})\cdot \delta (W-n_f^{{\lambda }_{\tau }}{E}_f^{{\lambda }_{\tau }}+n_i^{{\lambda }_0}{E}_i^{{\lambda }_0}).}

Here the probability of initial thermal states is described by a vector P(\mathit{I}) and P(\mathit{I})=(P( | 1^{{\lambda }_0}: n_1^{{\lambda }_0}⟩),P(|2^{{\lambda }_0}:n_2^{{\lambda }_0}⟩), \cdots ). \mathit{I}=(|1^{{\lambda }_0}:n_1^{{\lambda }_0}⟩,|2^{{\lambda }_0}:n_2^{{\lambda }_0}⟩,\cdots ) is the initial eigenstate vector and its corresponding initial boson number distribution is given by {\mathit{n}}_{\mathit{I}}^{{\mathit{\lambda }}_{\mathbf{0}}}=(n_1^{{\lambda }_0},n_2^{{\lambda }_0},\cdots ). \mathit{F}=(|1^{{\lambda }_{\tau }}:n_1^{{\lambda }_{\tau }}⟩,|2^{{\lambda }_{\tau }}:n_2^{{\lambda }_{\tau }}⟩,\cdots ) is the final eigenstate vector and its corresponding final boson number distribution is given by {\mathit{n}}_{\mathit{F}}^{{\mathit{\lambda }}_{\mathit{\tau }}}=(n_1^{{\lambda }_{\tau }},n_2^{{\lambda }_{\tau }},\cdots ). P(\mathit{F}|\mathit{I}) denotes the transition probability from the initial multi-boson eigenstate vector \mathit{I} to the final multi-boson eigenstate vector \mathit{F}. From Eq. (3), one can see clearly that the work distribution is mainly determined by two factors. One is the initial thermal distribution probability in Eq. (1), which can be easily calculated on the classical computer. The other is the transition probability P(\mathit{F}|\mathit{I}) between the multi-boson eigenstates, which is a classically difficult problem to calculate. Due to the interference of multiple identical bosons, the transition probability between multi-boson eigenstates can be constructed from the permanent of the transition amplitude matrix between single-boson eigenstates and it is expressed by [13,17,18]

(4)P(\mathit{F}| \mathit{I})={\displaystyle \prod _{i=1}{\displaystyle \frac{1}{n_i^{{\lambda }_0}!}}{\displaystyle \prod _{f=1}{\displaystyle \frac{1}{n_f^{{\lambda }_{\tau }}!}}|\text{Per}({\mathrm{\Lambda }}^{({\mathit{n}}_{\mathit{F}}^{{\mathit{\lambda }}_{\mathit{\tau }}},{\mathit{n}}_{\mathit{I}}^{{\mathit{\lambda }}_{\mathbf{0}}})}){|}^2.}}

Here matrix function Per(\mathrm{\Lambda }) represents the permanent of a matrix \mathrm{\Lambda }. \mathrm{\Lambda } = (⟨f^{{\lambda }_{\tau }}|\hat{U}|i^{{\lambda }_0}⟩ ) is a transition amplitude matrix between single-boson eigenstates. \hat{H}(t) is the Hamiltonian of the system, and \hat{U} denotes an evolutionary unitary operator of the system satisfying the time-dependent Schrödinger equation \mathrm{i}\hslash {\mathrm{\partial }}_t\hat{U}(t)= \hat{H}(t)\hat{U}(t ). {\mathrm{\Lambda }}^{({\mathit{n}}_{\mathit{F}}^{{\mathit{\lambda }}_{\mathit{\tau }}},{\mathit{n}}_{\mathit{I}}^{{\mathit{\lambda }}_{\mathbf{0}}})} denotes a sub-matrix of \mathrm{\Lambda } by taking n_f^{{\lambda }_{\tau }} (f=1,2,\cdots) copies of the fth column and n_i^{{\lambda }_0} (i=1,2,\cdots) copies of the ith row of \mathrm{\Lambda }. Since the total number of bosons is conserved, i.e., \sum _{i=1}{n}_i^{{\lambda }_0}=\sum _{f=1}{n}_f^{{\lambda }_{\tau }}=N, {\mathrm{\Lambda }}^{({\mathit{n}}_{\mathit{F}}^{{\mathit{\lambda }}_{\mathit{\tau }}},{\mathit{n}}_{\mathit{I}}^{{\mathit{\lambda }}_{\mathbf{0}}})} occupies a dimension of N\times N.

As shown in Fig.1, N identical bosons in a one-dimensional quantum piston are an interesting example to understand the work distribution. At t=0, the bosons are prepared into a thermal state in a stretchable box of length {\lambda }_0 [see Fig.1(a)], and their population in the eigenstates is schematically shown in Fig.1(c). Then, the box is stretched to the length {\lambda }_{\tau } ({\lambda }_{\tau }> {\lambda }_0) at a constant speed v [see Fig.1(b)] and at this time the population of bosons in the eigenstates is schematically presented in Fig.1(d). In Appendix A, we provide an explicit expression for constructing the transition amplitude matrix \mathrm{\Lambda } between the single-boson eigenstates [6,50]. Note that calculating the transition probability between the multi-boson eigenstates in Eq. (4) is a classically difficult problem, because the computational complexity of the permanent of a general complex matrix is #P-hard [20,21]. To solve the work distribution problem, in the following text, we would use the boson sampling to simulate the work distribution with multiple identical bosons.

2.2 Mapping work system into boson sampling

Boson sampling is considered as there are N indistinguishable bosons which are scattered into a linear unitary network with M optical modes. We denote the input photon state in the Fock basis as |\mathit{T}⟩=|t_1, t_2,\cdots ,t_M⟩. Each t_i denotes boson occupation-number in the ith optical mode and |\mathit{T}⟩ describes N= \sum _{i=1}^M{t}_i bosons distribution in each mode. These photons are sent through a linear optical network that is characterized by a unitary transformation A. According to the linear mapping relation a_i^{†} \to \sum _{j=1}^M{A}_{ij} b_j^{†} between the input mode creation operator a^{†} and the output mode creation operator b^{†} of the network, the probability of getting an output state | \mathit{S}⟩=|s_1,s_2,\cdots ,s_M⟩ in the Fock basis is mathematically described by [36-39]

(5)\begin{array}{rl}P(\mathit{S}|\mathit{T})=& |⟨\mathit{S}|A|\mathit{T}⟩{|}^2\\ =& \prod _{j=1}^M\frac{1}{s_j!}\prod _{i=1}^M\frac{1}{t_i!}|\text{Per}(A^{(\mathit{S}\mathbf{,}\mathit{T})}){|}^2.\end{array}

Remarkably, the probability P(\mathit{S}| \mathit{T}) for each of input states and output states is proportional to a permanent of sub-matrix of A. Combining Eq. (2), Eq. (4), and Eq. (5), we present a corresponding relationship between the work distribution and boson sampling in Tab.1. As shown in Tab.1, one can see that the space of the work distribution with N bosons and M-dimensional transition amplitude matrix is isomorphic to the space of boson sampling with N bosons and M-dimensional optical network. Therefore, the work distribution can be obtained by sampling from a great quantity of matrix permanents, equivalent to the boson sampling problem.

As shown in Fig.2, we design a schematic setup for simulating the work distribution of multiple identical bosons via boson sampling. We first carefully prepare the input Fock state to simulate the initial thermal state of the multiparticle bosonic system. The input Fock state |\mathit{T}{⟩}_n in the nth optical mode represents the nth energy eigenstate of the initial thermal state with the probability P_n that can be calculated through Eq. (1). The similar method to prepare the thermal state has been experimentally demonstrated [12]. We then construct the optical network based on the transition amplitude matrix \mathrm{\Lambda } between the single-boson eigenstates. The matrix elements ⟨f^{{\lambda }_{\tau }} | \hat{U}|i^{{\lambda }_0}⟩ described in Eq. (A5) can be calculated numerically, which depend on the parameters {\lambda }_0, {\lambda }_{\tau }, and v. The transition amplitude matrix \mathrm{\Lambda } is a unitary matrix in principle because the dimension of the matrix \mathrm{\Lambda } could be infinite and the matrix \mathrm{\Lambda } satisfies the normalization,

(6){\displaystyle \sum _i|⟨f^{{\lambda }_{\tau }}|\hat{U}|i^{{\lambda }_0}⟩{|}^2={\displaystyle \sum _i⟨f^{{\lambda }_{\tau }}|\hat{U}|i^{{\lambda }_0}⟩⟨i^{{\lambda }_0}|{\hat{U}}^{†}|f^{{\lambda }_{\tau }}⟩=1,}}

(7){\displaystyle \sum _f|⟨f^{{\lambda }_{\tau }}|\hat{U}|i^{{\lambda }_0}⟩{|}^2={\displaystyle \sum _f⟨f^{{\lambda }_{\tau }}|\hat{U}|i^{{\lambda }_0}⟩⟨i^{{\lambda }_0}|{\hat{U}}^{†}|f^{{\lambda }_{\tau }}⟩=1.}}

However, in reality we have to restrict the dimension of the matrix \mathrm{\Lambda } to be finite to encode the matrix into a finite dimensional unitary optical network of boson sampling, which causes the matrix \mathrm{\Lambda } to become near-unitary and introduces an encoding error. We can truncate the size of the matrix \mathrm{\Lambda } to make it as unitary as possible to reduce the encoding error. We evaluate the error by calculating the unitary fidelity of the matrix \mathrm{\Lambda }, which is defined as [51]

(8)\mathcal{F}={\displaystyle \frac{1}{d}}|\mathrm{t}\mathrm{r}\sqrt{I_d^{1/2}\sigma I_d^{1/2}}|.

Here d is the dimension of \mathrm{\Lambda } and I_d is a d-order identity matrix. \sigma = \mathrm{\Lambda }{\mathrm{\Lambda }}^{\prime } and {\mathrm{\Lambda }}^{\prime } is the Hermitian conjugate of \mathrm{\Lambda }. The encoding error \mathcal{E} is described by \mathcal{E}=1-\mathcal{F}. We numerically calculate the matrix elements and truncate the size of matrix \mathrm{\Lambda } when the encoding error is 0.5% with the fixed parameters {\lambda }_0=1, {\lambda }_{\tau }= 2 and a varied speed v. A relationship between the matrix dimension and the expansion speed v is plotted in Fig.3(a). We also present a relationship between the matrix dimension and the final length of the box in Fig.3(b). As shown in Fig.3, one can see that the dimension of the matrix \mathrm{\Lambda } increases almost linearly with the acceleration of the piston speed v or the final length {\lambda }_{\tau } of the box. Besides, we note that if the truncated \mathrm{\Lambda } matrix as a submatrix is embedded into a unitary matrix [52] or increasing the dimension of the truncated \mathrm{\Lambda } matrix, the encoding error can be further reduced. In general, the dimension of the optical network will become bigger as the encoding error decreases.

The next step is to encode the transition amplitude matrix \mathrm{\Lambda } into a linear optical network. There are mainly two configurations of optical network to realize arbitrary unitary matrix, one is the triangle-shaped network [53] and the other is the square-shaped network [54]. As shown in Fig.2(a), here we use the square-shaped network to realize the transition amplitude matrix, because the symmetry design of the network is more robust against the photon loss and has minimal optical depth and better stability [54]. As shown in Fig.2(b), the crossing between two optical modes a and a+1 in the interferometer consists of two 50:50 beam splitters and two phase shifters, which can be expressed mathematically by a matrix T_{a,a+1}( \theta ,\phi ) [54]. T_{a,a+1}(\theta ,\phi ) is called an elimination matrix and it is obtained by replacing the entries of an identity matrix with the same size as \mathrm{\Lambda } at the ath and (a+1 )th rows and the ath and (a+1)th columns with

(9)\left(\begin{array}{cc}{\mathrm{e}}^{\mathrm{i}\phi }\cos \theta & -\sin \theta \\ {\mathrm{e}}^{\mathrm{i}\phi }\sin \theta & \cos \theta \end{array}\right),

and the rest of the other entries remain unchanged. Based on Gaussian elimination method, the matrix \mathrm{\Lambda } can be diagonalized into a diagonal matrix D by multiplying a series of T_{a,a+1} and its inverse matrix T_{a,a+1}^{-1}. The matrix \mathrm{\Lambda } is realized physically in an optical network by choosing suitable values of parameters \theta and \phi of T_{a,a+1} and the phase values of diagonal matrix D at each output port. The resource overhead for encoding a d\times d unitary matrix into the optical network requires d(d- 1) 50:50 beam splitters and d^2 phase shifters. In boson sampling system, the phase shifters to realize the diagonal matrix D can be removed as only final photon number is sampled, which will not affect the result but can reduce the resource cost. From Fig.3, one can see that the total resource cost presents at a polynomial hierarchy with the expansion speed v of the piston or the final length {\lambda }_{\tau } of the box in terms of the number of required optical elements.

2.3 The effect of the system parameters on the work distribution

The temperature and the speed play an important role in the piston system during the work process. On one hand, from Eq. (1), one can see that the temperature of the system affects the initial distribution of the bosons. The work distribution for two bosons with different temperature is presented in Fig.4. As shown in Fig.4, as the temperature increases, the probability that the bosons populate higher energy levels will increase, which makes the higher energy levels need to be considered. As a result, the dimension of the transition amplitude matrix will become larger, making the calculation of the work distribution more complicated.

On the other hand, the moving speed of the piston is related to the transition probability P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩\to|f^{{\lambda }_{\tau }}:n_f^{{\lambda }_{\tau }}⟩). In the low speed limit, one can get the result P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩\to |f^{{\lambda }_{\tau }}:n_f^{{\lambda }_{\tau }}⟩)\to {\delta }_{i,f}(v\to 0) based on the quantum adiabatic theorem. Fig.5(a) shows approximately the initial energy distribution: the highest peak represents the energy from the beginning to the end of the ground state, the second highest peak corresponds to the energy of the first excited state, etc. As shown in Fig.5, with the increase of the piston speed, the work distribution for two bosons becomes more complex. This is caused by the higher energy level transition of the bosons when the speed becomes faster.

2.4 An example

To understand the work distribution simulated by boson sampling well, as an interesting example, we calculate the work distribution of three bosons in 1D quantum piston in detail. We consider three bosons are trapped in a box with the initial length {\lambda }_0=1, the stretching speed v=0.4, the final length {\lambda }_{\tau }=2, and the temperature \beta =0.1. We numerically calculate the matrix elements of \mathrm{\Lambda } based on Eq. (A5) and obtain a 5\times 5 dimensional near-unitary matrix {\mathrm{\Lambda }}_5 with a unitary fidelity \mathcal{F}=0.9992 (see Appendix B). We decompose the matrix {\mathrm{\Lambda }}_5 into the product of a diagonal matrix D and a series of elimination matrices T_{a,a+1} based on Gaussian elimination method [54]. The result of the decomposition is expressed as

(10)\begin{array}{r}{\mathrm{\Lambda }}_5=D{T}_{3,4}^{(5)}{T}_{4,5}^{(4)}{T}_{1,2}^{(5)}{T}_{2,3}^{(4)}{T}_{3,4}^{(3)}{T}_{4,5}^{(2)}{T}_{1,2}^{(3)}{T}_{2,3}^{(2)}{T}_{3,4}^{(1)}{T}_{1,2}^{(1)}.\end{array}

As shown in Fig.6, the matrix {\mathrm{\Lambda }}_5 is programmed into an optical network of the boson sampler. The values of phase shifter angles \theta and \phi in the network are calculated and presented in Tab.2. The diagonal matrix D can be ignored without affecting the final output probability of photons.

Finally, we calculate the cumulative work distribution by the definition as follows:

(11)\chi (W)={\displaystyle {\int }^W\rho (W^{\prime })\mathrm{d}{W}^{\prime }.}

The result of cumulative work distribution based on Eq. (11) is plotted in the curve of Fig.7. In real experiments, the noise and error are inevitable. Therefore, we evaluate the effect of the noise on the cumulative work distribution by adding a random noise \mathcal{N}\in (-0.01,0.01) to the angles of the beam splitters. Explicitly, we randomly choose 100 groups of noise terms and calculate the cumulative work distribution under the noise effects (see error bars in Fig.7). We find that the ratio between the error bar and the cumulative work distribution curve is 1% to 2%.

3 Feasibility analysis

Before we discuss the feasibility of the scheme, at first we clarify the relation between boson sampling and the matrix permanent. Boson sampling is a classically difficult problem as the sampled probability involves the permanent of a matrix. However, this does not mean that boson sampling can directly simulate the matrix permanent. Scaling up the system size of boson sampling, the probability of individual output event will become exponentially small, which causes that one has to obtain an exponential number of observations of the event to maintain the accuracy. This method of individual probability estimation is obviously infeasible. Remarkably, if one groups and sums the individual output probabilities so that the sum probability is polynomially small rather than exponentially small, then the number of needed samples would become polynomial size [29,47]. The GPE can deal with the work distribution problem well.

We next give the explicit method of GPE and evaluate the feasibility of the method by analyzing the required total sample numbers and the reasonable accuracy. Based on the delta function in Eq. (3), the output outcomes can be grouped and sum them when W=m_f{E}_f^{{\lambda }_{\tau }}- n_i{E}_i^{{\lambda }_0}, resulting in

(12)\begin{array}{rl}G(W)=& \sum _{\mathit{m}}P(\mathit{m}|\mathit{n})\delta (W-m_f{E}_f^{{\lambda }_{\tau }}+n_i{E}_i^{{\lambda }_0})\\ =& \sum _{\mathit{m}\in \mathcal{G}(W)}P(\mathit{m}|\mathit{n}).\end{array}

Here P(\mathit{m}|\mathit{n}) is the transition probability from the initial photon distribution \mathit{n}=(n_1 ,\dots ,n_M) to the final distribution \mathit{m}=(m_1,\dots ,m_M). The sets \mathcal{G}(W)= \{\mathit{m}\in{\mathbb{Z}}_{\ge 0}^M|m_f{E}_f^{{\lambda }_{\tau }}-n_i{E}_i^{{\lambda }_0}=W\} for each of grouped outcomes W\in \{ 0,\dots ,W_{\text{max}}\}. To analyse the feasibility of the GPE method, we define two integer energy vectors {\mathit{E}}_{\mathit{i}} \in {\mathbb{Z}}_{\ge 0}^M and {\mathit{E}}_{\mathit{f}} \in {\mathbb{Z}}_{\ge 0}^M, in which each of elements E_i and E_f (i,f=1,2,\dots ,M) are both at most a polynomial large number, i.e., E_i\le \mathcal{O}(\text{poly}(M)) and E_f\le \mathcal{O}(\text{poly}(M) ). In fact, each initial energy E_i^{{\lambda }_0} and each final energy E_f^{{\lambda }_{\tau }} in Eq. (3) can be expressed as a floating-point number and they can be transformed as the integer numbers E_i and E_f by multiplying by a sufficiently large number. Therefore, the number of different groups is W_{\text{max}}+1\le \mathcal{O}(\text{poly}(M)), which causes the probability | G(W)| of each group greater than or equal to \mathcal{O}(1/\text{poly} (M)) rather than the exponentially small.

Based on the result of grouped probability in Eq. (12) and the work distribution expression in Eq. (2), the estimated work distribution can be rewritten as

(13){\rho }_{\text{est}}(W) ={\displaystyle \sum _{i}P(|i^{{\lambda }_0}:n_i^{{\lambda }_0}⟩)G(W).}

This allows us to simulate the work distribution by grouping the output probabilities and collecting at most the polynomial samples of boson sampling.

We next introduce a reasonable accuracy ϵ to evaluate the performance of the estimated work distribution and the ideal one, i.e.,

(14)|{\rho }_{\text{est}}(W)-{\rho }_{\text{ide}}(W)|\le ϵ.

Based on the central limit theorem and Chebyshev’s inequality, the total sample number N_{\text{tot}} to achieve a reasonable accuracy scales as Var({\rho }_{\text{est}}(W ))/{ϵ}^2. It is clear that the variance of {\rho }_{\text{est}}(W) (denote as Var({\rho }_{\text{est}}(W ))) is bounded by 1 because the value of the work distribution {\rho }_{\text{est}}(W ) is between 0 and 1. In other words, N_{\text{tot}}= \mathcal{O}(1/{ϵ}^2) is an upper bound on the total number of samples required to simulate the work distribution, and this upper bound also implies that the reasonable target accuracy becomes \mathcal{O}(1/\text{poly} (M)), instead of the exponentially small for individual output probabilities. It suggests that polynomial large samples would be still reasonable.

4 Conclusion

In conclusion, we have presented a connection between the work distribution and boson sampling. Intuitively, calculating the work distribution may be difficult because the calculation of the transition probability between the multi-boson eigenstates is a classically hard problem. We found that boson sampling can be used to efficiently simulate the work distribution by sampling the output probability of photons and using the method of the GPE. We analyzed the computational cost with this set-up, and the results showed that at most a polynomial number of observation samples and optical elements are required to achieve a reasonable accuracy. The connection provides a new possible solution for studying the work distribution that is too difficult to calculate on a classical computer. For other systems, such as multiple bosons in a harmonic oscillator, the calculation of the work distribution also is not an easy problem because calculating the transition probability between the multi-boson eigenstates is necessary [13]. The scheme we developed is also suitable to simulate the work distribution of multiple identical bosons in both a contraction piston system and a harmonic oscillator system.

References

[1] C. Jarzynski. Nonequilibrium equality for free energy differences. Phys. Rev. Lett., 1997, 78(14): 2690

[2] G. E. Crooks. Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences. Phys. Rev. E, 1999, 60(3): 2721

[3] T. Hatano, S. I. Sasa. Steady-state thermodynamics of Langevin systems. Phys. Rev. Lett., 2001, 86(16): 3463

[4] D. Kafri, S. Deffner. Holevo’s bound from a general quantum fluctuation theorem. Phys. Rev. A, 2012, 86(4): 044302

[5] S. Deffner, E. Lutz. Nonequilibrium work distribution of a quantum harmonic oscillator. Phys. Rev. E, 2008, 77(2): 021128

[6] H. T. Quan, C. Jarzynski. Validity of nonequilibrium work relations for the rapidly expanding quantum piston. Phys. Rev. E, 2012, 85(3): 031102

[7] C. Jarzynski, H. T. Quan, S. Rahav. Quantum-classical correspondence principle for work distributions. Phys. Rev. X, 2015, 5(3): 031038

[8] L. Zhu, Z. Gong, B. Wu, H. T. Quan. Quantum-classical correspondence principle for work distributions in a chaotic system. Phys. Rev. E, 2016, 93(6): 062108

[9] M. Łobejko, J. Luczka, P. Talkner. Work distributions for random sudden quantum quenches. Phys. Rev. E, 2017, 95(5): 052137

[10] T. B. Batalhão, A. M. Souza, L. Mazzola, R. Auccaise, R. S. Sarthour, I. S. Oliveira, J. Goold, G. De Chiara, M. Paternostro, R. M. Serra. Experimental reconstruction of work distribution and study of fluctuation relations in a closed quantum system. Phys. Rev. Lett., 2014, 113(14): 140601

[11] S. An, J. N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z. Q. Yin, H. T. Quan, K. Kim. Experimental test of the quantum Jarzynski equality with a trapped-ion system. Nat. Phys., 2015, 11(2): 193

[12] D. J. Sivananda, N. Roy, P. C. Mahato, S. S. Banerjee. Exploring the non-equilibrium fluctuation relation for quantum mechanical tunneling of electrons across a modulating barrier. Phys. Rev. Res., 2020, 2(4): 043237

[13] Z. Gong, S. Deffner, H. T. Quan. Interference of identical particles and the quantum work distribution. Phys. Rev. E, 2014, 90(6): 062121

[14] Q. Wang, H. T. Quan. Understanding quantum work in a quantum many-body system. Phys. Rev. E, 2017, 95(3): 032113

[15] B. Wang, J. Zhang, H. T. Quan. Work distributions of one-dimensional fermions and bosons with dual contact interactions. Phys. Rev. E, 2018, 97(5): 052136

[16] J.GooldF. PlastinaA.GambassiA.Silva, The role of quantum work statistics in many-body physics, in: Thermodynamics in the Quantum Regime, Springer, 2018, pp 317–336

[17] M. C. Tichy, M. Tiersch, F. Mintert, A. Buchleitner. Many-particle interference beyond many-boson and many-fermion statistics. New J. Phys., 2012, 14(9): 093015

[18] M. C. Tichy. Interference of identical particles from entanglement to boson-sampling. J. Phys. At. Mol. Opt. Phys., 2014, 47(10): 103001

[19] J. D. Urbina, J. Kuipers, S. Matsumoto, Q. Hummel, K. Richter. Multiparticle correlations in mesoscopic scattering: Boson sampling, birthday paradox, and Hong‒Ou‒Mandel profiles. Phys. Rev. Lett., 2016, 116(10): 100401

[20] L. G. Valiant. The complexity of computing the permanent. Theor. Comput. Sci., 1979, 8(2): 189

[21] S.Aaronson, A linear-optical proof that the permanent is #P-hard, Proc. R. Soc. A 467(2136), 3393 (2011)

[22] A. P. Lund, M. J. Bremner, T. C. Ralph. Quantum sampling problems, boson-sampling and quantum supremacy. npj Quantum Inf., 2017, 3: 15

[23] S.AaronsonA. Arkhipov, The computational complexity of linear optics, in: Proceedings of the forty-third annual ACM symposium on Theory of computing, 2011, pp 333–342

[24] X. Gu, M. Erhard, A. Zeilinger, M. Krenn. Quantum experiments and graphs (II): Quantum interference, computation, and state generation. Proc. Natl. Acad. Sci. USA, 2019, 116(10): 4147

[25] J. M. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, J. Hundal, T. Isacsson, R. B. Israel, J. Izaac, S. Jahangiri, R. Janik, N. Killoran, S. P. Kumar, J. Lavoie, A. E. Lita, D. H. Mahler, M. Menotti, B. Morrison, S. W. Nam, L. Neuhaus, H. Y. Qi, N. Quesada, A. Repingon, K. K. Sabapathy, M. Schuld, N. Su, J. Swinarton, A. Száva, K. Tan, P. Tan, V. D. Vaidya, Z. Vernon, Z. Zabaneh, Y. Zhang. Quantum circuits with many photons on a programmable nanophotonic chip. Nature, 2021, 591(7848): 54

[26] J. Bao, Z. Fu, T. Pramanik, J. Mao, Y. Chi, Y. Cao, C. Zhai, Y. Mao, T. Dai, X. Chen, X. Jia, L. Zhao, Y. Zheng, B. Tang, Z. Li, J. Luo, W. Wang, Y. Yang, Y. Peng, D. Liu, D. Dai, Q. He, A. L. Muthali, L. K. Oxenlowe, C. Vigliar, S. Paesani, H. Hou, R. Santagati, J. W. Silverstone, A. Laing, M. G. Thompson, J. L. O’Brien, Y. Ding, Q. Gong, J. Wang. Very-large scale integrated quantum graph photonics. Nat. Photonics, 2023, 17(7): 573

[27] G. M. Nikolopoulos, T. Brougham. Decision and function problems based on boson sampling. Phys. Rev. A, 2016, 94(1): 012315

[28] G. M. Nikolopoulos. Cryptographic one-way function based on boson sampling. Quantum Inform. Process., 2019, 18(8): 259

[29] J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, A. Aspuru-Guzik. Boson sampling for molecular vibronic spectra. Nat. Photonics, 2015, 9(9): 615

[30] J. Huh, M. H. Yung. Vibronic boson sampling: Generalized Gaussian boson sampling for molecular vibronic spectra at finite temperature. Sci. Rep., 2017, 7(1): 7462

[31] Y. Shen, Y. Lu, K. Zhang, J. Zhang, S. Zhang, J. Huh, K. Kim. Quantum optical emulation of molecular vibronic spectroscopy using a trapped-ion device. Chem. Sci. (Camb.), 2018, 9(4): 836

[32] C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. M. Girvin, R. J. Schoelkopf. Efficient multiphoton sampling of molecular vibronic spectra on a superconducting bosonic processor. Phys. Rev. X, 2020, 10(2): 021060

[33] L. Banchi, M. Fingerhuth, T. Babej, C. Ing, J. M. Arrazola. Molecular docking with Gaussian boson sampling. Sci. Adv., 2020, 6(23): eaax1950

[34] J.ShiT. ZhaoY.WangC.YuY.Lu R.ShiS. ZhangJ.Wu, An unbiased quantum random number generator based on boson sampling, arXiv: 2206.02292 (2022)

[35] J. Shi, T. Zhao, Y. Wang, Y. Feng, J. Wu. Chaotic image encryption based on boson sampling. Adv. Quantum Technol., 2023, 6(2): 2200104

[36] M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T. C. Ralph, A. G. White. Photonic boson sampling in a tunable circuit. Science, 2013, 339(6121): 794

[37] J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Kolthammer, X. M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, I. A. Walmsley. Boson sampling on a photonic chip. Science, 2013, 339(6121): 798

[38] M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, P. Walther. Experimental boson sampling. Nat. Photonics, 2013, 7(7): 540

[39] A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvao, N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, F. Sciarrino. Integrated multimode interferometers with arbitrary designs for photonic boson sampling. Nat. Photonics, 2013, 7(7): 545

[40] H. Wang, Y. He, Y. H. Li, Z. E. Su, B. Li, H. L. Huang, X. Ding, M. C. Chen, C. Liu, J. Qin, J. P. Li, Y. M. He, C. Schneider, M. Kamp, C. Z. Peng, S. Höfling, C. Y. Lu, J. W. Pan. High-efficiency multiphoton boson sampling. Nat. Photonics, 2017, 11(6): 361

[41] H. Wang, W. Li, X. Jiang, Y. M. He, Y. H. Li, X. Ding, M. C. Chen, J. Qin, C. Z. Peng, C. Schneider, M. Kamp, W. J. Zhang, H. Li, L. X. You, Z. Wang, J. P. Dowling, S. Höfling, C. Y. Lu, J. W. Pan. Toward scalable boson sampling with photon loss. Phys. Rev. Lett., 2018, 120(23): 230502

[42] H. Wang, J. Qin, X. Ding, M. C. Chen, S. Chen, X. You, Y. M. He, X. Jiang, L. You, Z. Wang, C. Schneider, J. J. Renema, S. Höfling, C. Y. Lu, J. W. Pan. Boson sampling with 20 input photons and a 60-mode interferometer in a 10¹⁴-dimensional Hilbert space. Phys. Rev. Lett., 2019, 123(25): 250503

[43] H. S. Zhong, H. Wang, Y. H. Deng, M. C. Chen, L. C. Peng, Y. H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X. Y. Yang, W. J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. L. Liu, C. Y. Lu, J. W. Pan. Quantum computational advantage using photons. Science, 2020, 370(6523): 1460

[44] H. S. Zhong, Y. H. Deng, J. Qin, H. Wang, M. C. Chen, L. C. Peng, Y. H. Luo, D. Wu, S. Q. Gong, H. Su, Y. Hu, P. Hu, X. Y. Yang, W. J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. L. Liu, J. J. Renema, C. Y. Lu, J. W. Pan. Phase programmable Gaussian boson sampling using stimulated squeezed light. Phys. Rev. Lett., 2021, 127(18): 180502

[45] H. K. Lau, D. F. V. James. Proposal for a scalable universal bosonic simulator using individually trapped ions. Phys. Rev. A, 2012, 85(6): 062329

[46] C. Shen, Z. Zhang, L. M. Duan. Scalable implementation of boson sampling with trapped ions. Phys. Rev. Lett., 2014, 112(5): 050504

[47] C.OhY.Lim Y.WongB. FeffermanL.Jiang, Quantum-inspired classical algorithm for molecular vibronic spectra, arXiv: 2202.01861 (2022)

[48] H.Tasaki, Jarzynski relations for quantum systems and some applications, arXiv: cond-mat/0009244 (2000)

[49] P.TalknerE. LutzP.Hänggi, Fluctuation theorems: Work is not an observable, Phys. Rev. E 75, 050102(R) (2007)

[50] S. W. Doescher, M. H. Rice. Infinite square-well potential with a moving wall. Am. J. Phys., 1969, 37(12): 1246

[51] A. Gilchrist, N. K. Langford, M. A. Nielsen. Distance measures to compare real and ideal quantum processes. Phys. Rev. A, 2005, 71(6): 062310

[52] A.BjörklundB.GuptN.Quesada, A faster Hafnian formula for complex matrices and its benchmarking on the Titan supercomputer, J. Exp. Algor. 24, 11 (2019)

[53] M. Reck, A. Zeilinger, H. J. Bernstein, P. Bertani. Experimental realization of any discrete unitary operator. Phys. Rev. Lett., 1994, 73(1): 58

[54] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, I. A. Walsmley. Optimal design for universal multiport interferometers. Optica, 2016, 3(12): 1460

RIGHTS & PERMISSIONS

Higher Education Press