Trotterized adiabatic quantum simulation and its application to a simple all-optical system

As first proposed for the adiabatic quantum information processing by Wu et al (2002 Phys. Rev. Lett. 89 057904), the Trotterization technique is a very useful tool for universal quantum computing, and in particular, the adiabatic quantum simulation of quantum systems. Given a boson Hamiltonian involving arbitrary bilinear interactions, we propose a static version of this technique to perform an optical simulation that would enable the identification of the ground state of the Hamiltonian. By this method, the dynamical process of the adiabatic evolution is mapped to a static linear optical array which is robust to the errors caused by dynamical fluctuations. We examine the cost of the physical implementation of the Trotterization, i.e. the number of discrete steps required for a given accuracy. Two conclusions are drawn. One is that the number of required steps grows much more slowly than the system size if the number of non-zero matrix elements of Hamiltonian is not too large. The second is that small fluctuations of the parameters of optical elements do not affect the first conclusion. This implies that the method is robust against the certain type of errors as we considered. Last but not least, we present an example of implementation of the simulation on a photonic chip as well as an optimized scheme. By such examples, we show a reduction of the costs compared to its classical counterpart and the potential for further improvement, which promotes a more general application.


Introduction
The reason for simulating a quantum system using another quantum system is to obtain information about an uncontrollable system from a controllable one which is similar to the former. It has attracted a lot of attention ever since proposed by Richard P Feynman [1], and developed by Lloyd [2]. Recent studies [3][4][5][6][7][8][9][10][11][12][13][14][15] show that quantum simulation can provide alternative approaches to finding solutions by encoding them to the ground state of a Hamiltonian. Some of the simulation strategies have been proven to be capable of dealing with classically intractable problems, for example NP-complete problems [4,12].
One major obstacle to realizing the quantum simulation of a particular system is the difficulties in preparing the ground state of a Hamiltonian. In a number of quantum systems, it is relatively easy to find the ground state of some particular Hamiltonian, but very difficult to find the one required to solve a specific problem about which we are concerned. A great deal of effort has been expended developing the strategies and technologies for ground state preparation, both experimentally and theoretically [16][17][18][19][20]. Among those preparation strategies, adiabatic evolution is well-known for its applicability to many different types of systems. In principle, if one prepares the ground state of some Hamiltonian, then one can then obtain the ground state of a target Hamiltonian by starting with the ground state that one can prepare and slowly evolving the system from the prepared Hamiltonian to the desired one. Such a scheme is guaranteed by adiabatic theorem and now termed adiabatic quantum computing (AQC) [21]. AQC has been verified in several experiments [22][23][24][25][26] and is considered a promising candidate for universal quantum computing [27]. In the implementation of AQC, the crucial step is to adiabatically connect the problem Hamiltonian (whose ground state encodes the solution) with the initial, prepared Hamiltonian. Fortunately, the Trotterizaion technique provides a way to achieve such connection. With this technique, one can decompose the total evolution into short-time operations during which the system Hamiltonian is approximately time-independent for each step. The dynamical control of the system can be implemented by a sequence of such operations. This dramatically lowers the difficulty of realizing AQC by reducing the control requirements. In general, the whole Trotterized-AQC (TAQC) protocol can be described as follows [5]. (i) Prepare the ground state |y ñ 0 of Hamiltonian H 0 . (ii) Find the problem Hamiltonian H p whose ground state encodes the solution. (iii) Set the total Hamiltonian ( ) and ( ) = g t t Twhere t is the time and T is the period for the entire evolution. Then decompose the evolution operator into a sequence of steps using the Trotter-Suzuki formula, which is the key ingredient and given by is the evolution operator from 0 to T, k is a large integer so that t = T k is a small time segment, and  is time ordering operator. (iv) Finally, obtain the solution by measuring the state |y ñ f which is the simulation of | ( )| y y ñ = ñ U T ad 0 using U(T) implemented according to equation (1). For operators A and B and a sufficiently small δ, the Trotter-Suzuki formula implies ( ) It was introduced for the simulation of complex time-independent Hamiltonians in [2]. The application of the formula to an adiabatic strategy involving a time-dependent Hamiltonian in a TAQC protocol as described above, was first proposed in [5] and experimentally implemented in [25].
Here, we propose an optical implementation of a Trotterized adiabatic simulation-one type of TAQC algorithm using linear optical elements. Linear optics is a promising system for many types of quantum information processing. A logical qubit can be encoded in the polarization, frequency, spatial modes or other degrees of freedom of a photon which can be preserved for a relatively long time and is controllable [28][29][30][31][32][33][34][35][36][37][38]. Just as important for our purposes, the operations of the system are static so that the dynamics are discretized. We consider the matrix diagonalization problem which is classically classified as NP-hard. In order to obtain the solution, we propose a method for reaching the ground state of a boson Hamiltonian with arbitrary bilinear interactions. Finding the ground state of a Hamiltonian is QMA-hard [39] and, in special cases, reduces to the diagonalization problem in certain Hilbert subspaces. In our case, this reduction provides an important example of optical simulation. We analyze the dependence of the implementation cost, given by the Trotter Number (parameter k in the decomposition (1)), on the system size. We also study the effects of fluctuations of the parameters of the simulation by using a randomized trotter formula (RTF) [12]. The definition of RTF is and g a is a random number. When g a is deleted, the decomposition (2) reduces to the standard one (1). In our case, the fluctuation of τ corresponds to imperfections of experimental optical elements. Unlike some of the recent investigations on evolution errors [40,41] which affect the structures of Hamiltonians, τ a represents an inaccurate evolution time in each piece rather than the perturbation on the structure or parameters of the Hamiltonians, meaning that the basic functioning of the optical elements is preserved. We show numerically that such error will add little extra cost to the simulation for a given accuracy. Also, we would like to point out that the scaling of the ideal settings will likely be decimated by other types of implementation errors such as those studied by [40,41]. This paper is organized as follows. In section 2, we introduce our simulation proposal of AQC via a static optical circuit. We present our analysis of the dependence of the Trotter number on the system size in section 3. We then provide three examples and show the numerical results of them in section 4. Section 5 illustrates an implementation using photonic chip. Section 6 demonstrates a potential optimized scheme for the Trotterization method enabled by our proposal. Section 7 concludes.

Adiabatic evolution simulated via optical elements
We consider a general model with .
and J mn is the coupling coefficient between the mth mode and nth mode. In the one-photon subspace, the Hamiltonian (3) can represent a matrix which has no additional constraints other than being Hermitian. So the process of finding its ground state is equivalent to diagonalizing a general Hermitian matrix. In our proposal, the bosonic modes are mapped to the spatial modes of photons. Hence, † b i corresponds to a photon propagating along an optical path labeled by i, and b i corresponds to the absence of the photon from the path. To implement a TAQC, one must design a physical realization of the adiabatic evolution. We now discuss the details of such a realization. First, applying the decomposition (1) to the evolution of a Hamiltonian (3), we have We can utilize the Hermiticity of where J Re mn ( J Im mn ) is the real (imaginary) part of J mn . Given the commutators and Trotter-Suzuki formula, every multiplier of expression (4) can be separated into three exponential operators, and each one can be further decomposed as  Next, we demonstrate how to simulate these operators using an array of optical devices. However, we note that it is possible to implement the same set of elements using a photonic chip [35,42,43] which should enable a speed-up over the classical counterpart. Here we introduce the basic methods and present a prototype for such a chip. An example of the chip design is provided in the section 5. We primarily use two common linear elements, phase shifters (PSs) and beam splitters (BSs), shown by figures 1(a) and (b). The mathematical descriptions of PS and BS are (See for example [44].) † c and † d are two different spatial modes. f is the phase shifted by a PS and θ defines the reflection (transmission) rate of a BS through q and ( ) q U bs match the form of equations (6) and (7). The factors of the forms ( . Superscripts s and l denote the optical modes. Factor ( ) Im mn m n m n can be implemented by one BS, , where m and n denote optical modes. The factor ( ) i R e mn m n m n can be implemented by a combination of four PSs and a BS, and the connection between the Lie group SU(2) and boson operators. An illustration of the above combination is given by figure 1(c).The whole implementation of the simulation is described by figure 2. For ease of illustration, figure 2 only shows nearest-neighbor interactions. However, it is in principle possible to implement any type of bilinear interaction.

Parameter dependence analysis
Our objective is the simulation of large quantum systems which are difficult to simulate using classical computers. However, when system size grows, more resources may be required to obtain the same level of simulation accuracy. Therefore, it is important to examine the variation of the resources with the system size. We next investigate this resource dependence in terms of the number of required segments (k) when the number of bosonic modes (N) increases. Although the error in the Trotterized time evolution for time-dependent Hamiltonians has been analyzed in previous work [45][46][47], we here provide a different perspective to look into the problem, which is directly related to the optical system studied here.
We note that, as shown by the decomposition (1) and the expression (4), the accuracy of simulation increases when Trotter number k grows. Also, the number of optical elements required to perform the simulation is proportional to k (see figure 2). Now consider the difference between the ideal adiabatic evolution and the Trotterized one as measured by | | | y y D = -á ñ 1 ad f 2 . The function U d (T), which is the discrete form of U(T) obtained using a finite-difference Schödinger equation, is sin , sin cos . Obviously, ( ) Also, consider the commutator of H 0 and H p which, to a large extent, describes the error when applying equations (6) and (7). By Taylor expansion, we can find the difference of . Then, the dominant factor of Δ can be calculated by adding up the leading terms in above expressions. More specifically, from U(T) to U d (T), we have the leading error given by (4)), the leading error is given by which is, in principle, sufficient to describe thek N relation, the leading error is given by 3 . In the second to last approximation, higher order terms are neglected. After some simplification, we pick out the leading terms and obtain where E pg (E 0g ) is the ground state energy of H p (H 0 ). This expression comes from the fact that |y ñ ad (|y ñ 0 ) is the ground state of H p (H 0 ). Because H 0 is diagonal, E 0g is independent of N. In general, E pg is a function of N determined by the structure of H p . The real part of the overlap | y y á ñ ad 0 is bounded by one. So we can rewrite Δ in the following form: where A, B, C are constants independent of N. We can conclude from equation (14) that for a given Δ, the dependence of the Trotter number k on system size N is determined by the ground state energy of H p to the leading order in the approximation. Thus determining the form of H p will enable the determination of the dependence of E g on N and therefore the relation between k and N for a given Δ. Notice that we only focus on the simulation accuracy of an AQC procedure based on the Trotter formula, rather than the effectiveness of the AQC when preparing the ground state. The latter is determined by the evolution time T as well as the band structure of the total Hamiltonian, while the former rests on the Trotterization precision affected by the decomposition of the exponential operator. So, in contrast to standard AQC, equation (14) involves the square of the ground state energy. Next, we consider the effect of the fluctuation of the optical elements on this relation. The analysis process is basically the same as before, except that we replace τ by τ a (given by RTF (2)). The fluctuation is modeled by zero-mean random number g a where a is an integer subscript. Such a replacement generates more terms than equations (9)- (11). The resultant expression of the additional terms to leading order are given by The consideration of RTF only changes the duration of a single part of evolution, so the parameters of the Hamiltonian in one part do not change. The error scaling of equation (15) can be analyzed by exploring the magnitude of the summations. Based on the following: (1) the H 0 (t) and H p (t) here are linear in t and (2) g a is a small amplitude fluctuation, each summation is estimated using the averages of the random numbers over finite samples. More specifically, one has

Numerical results
In this section, we numerically evaluate this dependence for some particular cases to show that our approximations are justified. Unlike the investigations in [45][46][47] which focus on the norm of the Hamiltonian matrices, we characterize the influence by ground state energy as we deduced in the last section. In the specific examples, the structures of the problem Hamiltonians are categorized by the density (or equivalently by its sparsity). This provides a way to estimate the relative computational complexity of a matrix calculation and is often used in computational physics. We first let H 0 be a diagonal matrix whose entries are sorted, equallyspaced and range from 0.5 to - The distribution of non-zero off-diagonal elements of H p is described by the density of a matrix, which is defined as the number of non-zero matrix elements divided by the total number of matrix elements.
We consider three off-diagonal examples of H p . The first one only involves the nearest-and next-nearestneighbor interaction, i.e. H p is a pentadiagonal matrix with density ( ) -N N 5 6 2 . The second off-diagonal part forms a sparse matrix with fixed density 0.5. The locations of non-zero entries are random. The third one is the case with full non-zero-off-diagonal entries which means the total density is 1. The values of the off-diagonal entries in all three types of H p randomly vary from 0 to 1. The simulation results are shown in  (14). We simulate the change of E pg with N in each of the three cases. For the first type of H p (bottom, blue), E pg only varies when N increases, so k is nearly constant. For the other two, E pg 2 is found to be linear for both cases and the slope of the fitting lines are 9.20×10 −4 for the third and 5.09×10 −4 for the second, which is approximately half of the former. The linearity of E pg 2 also indicates thatẼ N pg which means that E pg 2 is dominant, especially when N is large. Therefore, from figure 3, we can see that the slope of the third (top, red) is approximately twice the slope of the second (middle, green). The solid lines and the dashed lines are nearly coincident which supports our analysis of the fluctuations. In general, such a linear dependence is compatible with the results in [47].

An implementation of the Trotterization circuit using photonic chip
We utilize the following example to show that our proposal can be implemented using a quantum chip composed of fewer gates than the corresponding classical circuit. If the slowly-varying control functions are ( ) =f t t T 1 and ( ) = g t t T, as we consider in the main text, the total circuit can be built using two blocks. The first one, denoted by A, is the first piece of the Trotterization circuit (figure 2), which simulates the operator † i s s s s . The second one, denoted by B, simulates the difference of the two adjacent pieces, given by s s s s l l l l m n mn m n . According to our proposal, block A can be implemented by using only PSs and block B can be implemented by a combination of PSs and BSs (similar to a piece in the Trotterization circuit). Therefore, the ath piece in the Trotterization circuit can be realized by arranging block A and B based on the scheme shown by figure 4.
In such a manner, the Trotterization circuit we proposed can be transformed into a simpler one containing only block A and B, together with a set of control strategies to route the photons. To see that such controls are available, see the literature describing recent work on photonic chips [42,43]. An example is presented in figure 5. The photon state propagates in optical fibers or waveguides. The photon-number-resolving detector (PNRD) is used to obtain the photon information without destroying their states. Then, a control program can  adjust the phase of the photon via a tunable phase shifter based on the photon information. Therefore, the photon can be routed to block A or B by a Mach-Zehnder interferometer as required by the algorithm. A long optical distance might be used to make the photon state analysis in the router more convenience.
In the implementation shown in figure  if we restrict our consideration to two-mode interactions. Among the elements, the router for each mode requires one PS, two PNRDs and three BSs. Therefore, the total number of elements is 6N. Hence, number of gates in the circuit in figure 5 is , which implies the time complexity of the circuit is O(N 2 ). This is smaller than the corresponding classical algorithm which has a lower bound of O(N 3 ).

An implementation of adiabatic simulation using high-order Trotter-Suzuki formulas
Recently, several strategies for improving the performance of the Trotterization method have been proposed. In the proposals, high-order Trotter formulas as well as linear combinations of are introduced to decompose the adiabatic evolution of the Hamiltonian, which can boost the quality of the simulation or computation [48][49][50].
Here, we would like to point out that these can be directly implemented by our scheme. We illustrate our setup below by considering a simple example.
The Taylor series of the exponential operators shows that the second order error of the  [48]. In our case, the scaling of the error in each Trotter step can be suppressed as discussed previously. Therefore, the total number of Trotter steps for a required accuracy will be smaller. We numerically evaluate the cases in section 4 with following replacement e e e e . 1 7

t T H t T H t T H t T H t T H t T H
The number of Trotter steps when applying the replacement for different types of H p are all one order smaller than the results in section 4. Also, such a strategy is easy to be implemented in an optical system. As shown in section 5, one Trotter piece can be decomposed into a series of blocks. Realizing the evolution indicated by replacement (17) merely requires an additional series of the blocks arranged in a different order. Hence, the gate number only increases by a constant factor. A schematic of the setup is shown by figure 6.
There are other forms of high-order Trotter-Suzuki formulas as well. Some applications have been discussed in detail [49,50]. Here, we would like to emphasize that our optical implementation of the adiabatic simulation is compatible with any strategy composed of linear operations for error suppression or speed-up. One interesting benefit of the linear strategy is that the additional implementation costs can be limited to polynomial scaling. Fortunately, optical systems naturally capture the basic requirements for linear information processing, providing a good platform for implementing such strategies.

Conclusions
AQC provides an important potential avenue for quantum computing. However, some issues still need to be addressed. For example, there is no evidence that AQC can solve optimization problems (even in the ideal case) more efficiently than classical algorithms. Therefore, there is still the question of how one can improve Figure 6. A schematic demonstration of implementing replacement (17). Every Trotter step is completed by two similar blocks with different orders. All steps are connected by two beam splitters. Like figure 5, only the propagation of one mode is marked out. Block A and Block B are the same as those in figure 4. It is worth mentioning that above is only for a functional display. The scheme can also be integrated into a quantum chip as shown in figure 5. optimization algorithms running on AQC. Another very important challenge is to find suitable platforms for AQC and/or adiabatic simulation. This has been our main concern in the present work.
We proposed a scheme to adiabatically reach the ground state of a boson Hamiltonian with arbitrary bilinear interactions using Trotterization. The whole process is implemented by a linear optical design which is robust against errors caused by fluctuations in the accuracy of the individual elements. To the best of our knowledge, Trotterization is the only way a dynamical quantum process can be completely mapped to a static circuit. We also analyzed the dependence of implementation cost on the system size when the simulation accuracy is approximately fixed. Corresponding analytical and numerical results show that the cost of the simulation, represented by Trotter number k, and system size N is determined by the structure of problem Hamiltonian. When the structure is relatively simple, such as the case when H p is relatively sparse, the cost will grow more slowly than system size. Moreover, we found that imperfect experimental conditions, modeled by the fluctuation of evolution time, do not significantly affect the relation between k and N which means that the simulation is robust against such type of errors. This was verified numerically.
Following our theoretical and numerical analysis, we presented an example of the implementation of the proposal on a photonic chip. This shows that our scheme can be performed by a set of gates of order N 2 , which is a clear reduction of the costs compared to the algorithm's classical counterpart. Last but not least, we proposed to implement an optimized scheme for our simulation proposal based on high-order Trotter formula. By example, we showed that this kind of scheme is effective and can be directly realized by our design which indicates a further improvement of the all-optical adiabatic simulation.