Algorithmic approach to simulate Hamiltonian dynamics and an NMR simulation of Quantum State Transfer
We propose an iterative algorithm to simulate the dynamics generated by any -qubit Hamiltonian. The simulation entails decomposing the unitary time evolution operator (unitary) into a product of different time-step unitaries. The algorithm product-decomposes in a chosen operator basis by identifying a certain symmetry of that is intimately related to the number of gates in the decomposition. We illustrate the algorithm by first obtaining a polynomial decomposition in the Pauli basis of the -qubit Quantum State Transfer unitary by Di Franco et. al. (Phys. Rev. Lett. 101, 230502 (2008)) that transports quantum information from one end of a spin chain to the other; and then implement it in Nuclear Magnetic Resonance to demonstrate that the decomposition is experimentally viable and well-scaled. We furthur experimentally test the resilience of the state transfer to static errors in the coupling parameters of the simulated Hamiltonian. This is done by decomposing and simulating the corresponding imperfect unitaries.
pacs:03.67.Ac, 76.60.-k, 03.67.Lx
Introduction: Feynman Feynman (1982) has stated that it should be possible to manipulate the Hamiltonian of one quantum system to simulate the dynamics of another, exponentially faster than a classical computer. This serves as one of the main motivations for building a quantum information processor (QIP) Lloyd (1996). More precisely, simulation is the desire to mimic the unitary dynamics of a -qubit Hamiltonian ; and if the physical resources required for the simulation scale polynomially with the problem size , then it is said to be efficient Lloyd (1996).
The Hamiltonian of any QIP is the sum of an intrinsic part and a time dependent part Khaneja et al. (2001) that can be experimentally controlled such that the following holds to arbitrary precision Kitaev (2002):
where is the simulation time and is the Dyson time-ordering operator. A QIP allows universal control if it can simulate any unitary dynamics Vlasov (2001). However finding the requisite efficient set of controls is in general a challenge Ryan et al. (2008).
In practice, the simulation is constructed by decomposing into a product of unitary evolution operators : Kitaev (2002) for a sequence of time-steps with . A polynomial product decomposition (PPD) of is a decomposition such that scales polynomially with the problem size . Clearly, a PPD is necessary for a simulation to be efficient. Efficiency also entails that each (or gates) be implemented so that the amount of physical resources (spatial cost) and the implementation time (temporal cost) are small – i.e. the total cost incurred scales polynomially with . Moreover, inevitable decoherence processes in the QIP Viola et al. (1999) could impose further efficiency constraints.
For certain unitaries a PPD may not exist in principle Kitaev (2002). In this case, or when the PPD of is not known, one resorts to approximate methods. This entails resolving into a finer time-steps , and then by assuming one expands onto a product decomposition (PD). Consequently, both the total cost (the number of time steps) and the precision become dependent: if one increases (decreases) the precision, then the total temporal cost also increases (decreases). For instance, the use Lloyd (1996) of Suzuki-Trotter expansion Suzuki (1990) limits the precision to Vlasov (2001). This has motivated numerical optimization methods to achieve the desired balance between the precision and time Khaneja et al. (2005); Fortunato et al. (2002), but their inherent computationally intensive nature restricts them to small , and they generally cannot be easily extended to arbitrary .
In this paper, we not only propose an algorithm to exactly product-decompose any , but perhaps more importantly the algorithm also allows one to search for an efficient PPD. The algorithm first studies the support of in a basis B by representing it as a vector , since , where and . A PD, , is then obtained by systematically using the unitaries (gates) generated by the basis operators to iteratively rotate . The PD is derived for arbitrary in an inductive manner, by extrapolating the symmetry in the PD for small . In principle, our algorithm allows any basis at the start: if a PPD is not obtained for one choice, then a different basis maybe tried. Note however, that any random search for such bases is bound to be inefficient since the optimal search algorithm of Grover Grover (1996) is non-polynomial and hence inefficient. Since a PPD is realized when the PD consists of a polynomial number of gates, this may be interpreted to mean that has a large symmetry in the sense that it is the composition of a “few” rotations: the Quantum Fourier Transform unitary Simon (1997); Weinstein et al. (2001) being a famous example.
We note that our algorithm also enables one to understand how errors and noise affect the simulation of the unitary . The errors manifest themselves by reducing the symmetry of hence increasing the size of the PD. This can be used to develop appropriate techniques to control the errors.
Other algorithms for simulation have been suggested previously. Ref. Vartiainen et al. (2004) for example uses Given’s rotations Givens (1958); Householder (1958) to product-decompose in terms of CNOT and controlled phase gates. Here, we have explored the usefulness of the Pauli operator basis Sorensen et al. (1984) for the decomposition. The advantage is that gates from this basis can be implemented in certain spin-based QIPs in a time optimal manner Khaneja et al. (2001, 2002). Moreover, other optimal control techniques Khaneja et al. (2005); Fortunato et al. (2002) can also be modularized in our algorithm to further reduce the total cost of simulation. Note, however, that the Pauli basis will not lead to a PPD for all , and for a given QIP, there may exist several efficient bases. Recently Samal et al. (2011) employed our algorithm to experimentally simulate the Quantum No-Hiding theorem Pati and Braunstein (2000).
We explain the essential ingredients in the algorithm, including the role of symmetries, by explicitly product-decomposing (and simulating) the unitary that causes quantum state transfer (QST) in a linear spin- chain Christandl et al. (2004). QST allows the chain to act equivalent to a “wire” in a spin-based QIP architecture Bose (2003). Simulation is motivated because the spin chains, as required for QST, are normally hard to manufacture Bose (2003). We first obtain a PPD in the Pauli basis of a remarkable QST protocol by Di Franco et al Di Franco et al. (2008) that inherently doesn’t require the chain to be initialized Di Franco et al. (2008); Cappellaro et al. (2007). We then experimentally simulate it in a Nuclear Magnetic Resonance (NMR) QIP Mádi et al. (1997); Zhang et al. (2006, 2007).
We also use the PPD to investigate the protocol’s robustness against certain kind of errors Strauch and Williams (2008); Zwick et al. (2011). This is done by introducing specific errors in the protocol which is then simulated by our algorithm. We explicitly demonstrate how the error results in a PD that scales faster than the PPD of the error-free protocol.
The QST protocol Di Franco et al. (2008): Given a -qubit chain and couplings between the qubits and , the Ising Hamiltonian of the chain is given as,
where denote spin- Pauli matrices, eg. . Let be the initial state of the chain, then QST is achieved by evolving for :
where . A state locally equivalent to is then recovered at the end of the chain by measuring the first qubit in the Z-basis.
The algorithm: Consider a -qubit as an example. It can be represented as a column vector in the -dimensional vector space with basis ; and where the subscript denotes that it supports the hyperplane where the vector lies. The coefficients can be deduced via the orthogonality of the basis: eg. . More importantly, Clifford algebra ensures that the basis forms a group up to a phase, which we denote as . We define the norm (squared length) of in space as,
We recursively use the fact that one can express in an effective two-component form by choosing two orthogonal subspaces of , for instance, and . is necessarily chosen to be a proper subgroup of , and is a coset. This allows us to decompose the norm as,
The algorithm obtains the PD in terms of the unitaries (gates) generated by the basis elements of B. To do so, it rotates to , and since is closed under the product operation, the elements of have no further role to play. The desired rotations are obtained by first product-decomposing (for example) into , where the gates and are chosen with appropriate angles and such that has no support in – we call such gates orthogonal rotations. To obtain the angles (equivalently, the time of evolution), note that the product in the vector representation is simply the following linear transformation of ,
where is a permutation matrix such that . We first choose that maximizes , the increase in norm of in :
where is a reflection matrix about . The optimal is then a function of Ajoy et al. (2011):
since . Similarly when . The orthogonal rotations rotate to . The PD is completed by obtaining such that .
The above can be generalized to a -qubit unitary . Now basis B contains operators, and to within a phase, forms a group under product operation. In summary, the PD of the is achieved by iteratively rotating the vector via orthogonal rotations to smaller and smaller (subgroup) subspaces until the PD is completed.
The rotations are essentially deduced by diagonalizing in two-component form, analogous to determining the eigenvectors in “coupled two-level problems” in quantum mechanics. More importantly, the order of the rotations are chosen to minimize the number of terms (gates) in the PD. This is done through dynamic programming Bellman (2003): at each step, the orthogonal rotation that causes the maximum transfer of norm to the desired subspace is chosen first (or, the that has the maximum “off-diagonal” contribution to the two-component matrix is extinguished first).
In order to minimize the time to obtain a PD for a of arbitrary size, we will product decompose for small sizes, and identify the structural symmetry in the form of the PD. We will then use the pattern identified to inductively extended the PD to any size. We now demonstrate the above by explicitly product-decomposing the QST unitary.
The PPD of QST unitary: First consider the -qubit QST , which can be represented as the vector in the space spanned by :
Divide into two orthogonal subspaces: a subgroup , and its coset space . The algorithm generates a single orthogonal rotation in the space that rotates to . Now a further iteration gives the following PD:
Similarly, the PD for the -qubit ,
The above PD’s are seen to exhibit the following pattern: the QST is mirror symmetric about the center of the chain Karbach and Stolze (2005). For any equal division of into and , there exists a single orthogonal rotation in that causes maximum transfer of the norm such that the resulting unitary has no support in the coset space. This property is invariant under further iterations. Thus the decomposition of scales linearly with . Inductively, the PPD of an -qubit is given as,
NMR simulation of : We used liquid state NMR to simulate the QST protocol Zhang et al. (2006, 2007); Cappellaro et al. (2007) in a -qubit system , where H, C, and F are the qubits. Our motivation is to study experimental viability of the PD and its scaling. Experiments were performed at room temperature in a 11.7 T magnetic field with the resonance frequencies MHz (H), MHz (C), and MHz (F). The couplings between the qubits were Hz, Hz and Hz.
The gates in the PPD (Eq. 9) are implemented by using NMR pulse sequences Sorensen et al. (1984) as the controls. To do so, note that the rotations generated by any basis operators, such as , can be realized by using Vlasov (2001): . Thus orthogonal rotations by elements of B are implemented using hard pulses alone, by Hahn echo refocusing Braunschweiler and Ernst (1983); Tseng et al. (1999), or by a combination of both. Piecing together these sequences leads to the implementation of (Fig. 1) Ajoy et al. (2011). In Eq. 9 the multiqubit operators act only on contiguous qubit positions. Hence they can be constructed only using nearest-neighbor couplings. The non nearest-neighbor couplings are refocused during the simulation time and H, C and F (in that order) form the requisite -qubit chain.
Let be the initial state of the first qubit, which is to be transferred along the spin chain. Although the protocol allows the state of the second qubit to be arbitrary, we fix its initial state to for ease of measurement. Thus the initial state of the chain is , and is created from the psuedopure state (created by spatial averaging Cory et al. (1998)) using a X pulse on the first qubit. The pulse sequence corresponding to the PPD of (Fig. 1) is then applied. Since projective measurement is not possible in an bulk-ensemble QIP like NMR Ryan et al. (2008), the readout of the final state of the spin chain is done by applying a CNOT gate (see Price et al. (1999)) to force the system to . The state of the third qubit is measured with the receiver coil in the Y direction. The results (Fig. 1i) shows excellent agreement with Ajoy et al. (2011), which indeed confirms the QST.
To see the role of symmetry in the study of the robustness of the PPD (thus of the protocol) to coupling errors, we consider the simplest case of static disorder in the engineered couplings Strauch and Williams (2008); Burrell et al. (2009), where the deviations of the couplings from their ideal values are independent of qubit position, i.e . A representative example of the pulse sequence for this simulation is shown in Fig. 1 Ajoy et al. (2011). The experimental results (Fig. 1ii) show that the state recovered is corrupted by an additional phase, and the degree of corruption increases with the error of couplings. Hence, the transfer fidelity – characterized by how closely the final state reproduces the initial state – decreases with the error of couplings.
Fig. 1 demonstrates that the PD of the protocol with coupling error destroys the structural symmetry Bacon and Dam (2010) of the error-free PD Ajoy et al. (2011). Consequently, the size of the PD scales faster with the system size, hence incurring higher simulation cost. In conclusion, our work strongly suggests that the the structural symmetry of in a given basis is intimately related to the scaling of the corresponding PD. Any algorithm that unravels the symmetry in can therefore most economically simulate it.
Authors thank A.D. Patel and P. Cappellaro for discussions, and K. P. Yogendran for editing the manuscript.
- Feynman (1982) R. P. Feynman, Int. J. Theor. Phys 21, 467 (1982).
- Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
- Khaneja et al. (2001) N. Khaneja, R. Brockett, and S. J. Glaser, Phys. Rev. A 63, 032308 (2001).
- Kitaev (2002) A. Kitaev, Russian Math. Survey 52, 1191 (2002).
- Vlasov (2001) A. Y. Vlasov, Phys. Rev. A 63, 054302 (2001).
- Ryan et al. (2008) C. A. Ryan, C. Negrevergne, M. Laforest, E. Knill, and R. Laflamme, Phys. Rev. A 78, 012328 (2008).
- Viola et al. (1999) L. Viola, E. Knill, and S. Lloyd, Phys. Rev. Lett. 82, 2417 (1999).
- Suzuki (1990) M. Suzuki, Phys. Lett. A 146, 319 (1990).
- Khaneja et al. (2005) N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Heruggen, and S. J. Glaser, J. Magn. Res 172, 296 (2005).
- Fortunato et al. (2002) E. Fortunato, M. Pravia, N. Boulant, G. Teklemariam, T. Havel, and D. Cory, J. Chem. Phys. 116, 7599 (2002).
- Grover (1996) L. Grover, in Proceedings of 28th Annual ACM Symposium on Theory of Computing (STOC) (1996), pp. 212–219.
- Simon (1997) D. R. Simon, SIAM J. Comput. 26, 1474 (1997).
- Weinstein et al. (2001) Y. S. Weinstein, M. A. Pravia, E. M. Fortunato, S. Lloyd, and D. G. Cory, Phys. Rev. Lett. 86, 1889 (2001).
- Vartiainen et al. (2004) J. J. Vartiainen, M. Möttönen, and M. M. Salomaa, Phys. Rev. Lett. 92, 177902 (2004).
- Givens (1958) W. Givens, J. SIAM 6, 26 (1958).
- Householder (1958) A. S. Householder, J. ACM 5, 339 (1958).
- Sorensen et al. (1984) O. Sorensen, G. Eich, M. Levitt, G. Bodenhausen, and R. Ernst, Progress in Nuclear Magnetic Resonance Spectroscopy 16, 163 (1984).
- Khaneja et al. (2002) N. Khaneja, S. J. Glaser, and R. Brockett, Phys. Rev. A 65, 032301 (2002).
- Samal et al. (2011) J. R. Samal, A. K. Pati, and A. Kumar, Phys. Rev. Lett. 106, 080401 (2011).
- Pati and Braunstein (2000) A. Pati and S. Braunstein, Nature (London) 404, 164 (2000).
- Christandl et al. (2004) M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
- Bose (2003) S. Bose, Phys. Rev. Lett. 91, 207901 (2003).
- Di Franco et al. (2008) C. Di Franco, M. Paternostro, and M. S. Kim, Phys. Rev. Lett. 101, 230502 (2008).
- Cappellaro et al. (2007) P. Cappellaro, C. Ramanathan, and D. G. Cory, Phys. Rev. Lett. 99, 250506 (2007).
- Mádi et al. (1997) Z. Mádi, B. Brutscher, T. Schulte-Herbrüggen, R. Brüschweiler, and R. Ernst, Chemical Physics Letters 268, 300 (1997).
- Zhang et al. (2006) J. Zhang, X. Peng, and D. Suter, Phys. Rev. A 73, 062325 (2006).
- Zhang et al. (2007) J. Zhang, N. Rajendran, X. Peng, and D. Suter, Phys. Rev. A 76, 012317 (2007).
- Strauch and Williams (2008) F. W. Strauch and C. J. Williams, Phys. Rev. B 78, 094516 (2008).
- Zwick et al. (2011) A. Zwick, G. A. Álvarez, J. Stolze, and O. Osenda, Phys. Rev. A 84, 022311 (2011).
- Ajoy et al. (2011) A. Ajoy, R. K. Rao, A. Kumar, and P. Rungta, Supplementary online material (2011).
- Bellman (2003) R. Bellman, Dynamic Programming (Dover Publications, 2003).
- Karbach and Stolze (2005) P. Karbach and J. Stolze, Phys. Rev. A 72, 030301 (2005).
- Braunschweiler and Ernst (1983) L. Braunschweiler and R. Ernst, Journal of Magnetic Resonance (1969) 53, 521 (1983).
- Tseng et al. (1999) C. H. Tseng, S. Somaroo, Y. Sharf, E. Knill, R. Laflamme, T. F. Havel, and D. G. Cory, Phys. Rev. A 61, 012302 (1999).
- Cory et al. (1998) D. Cory, M. Price, and T. Havel, Physica D 120, 82 (1998).
- Price et al. (1999) M. D. Price, S. S. Somaroo, A. E. Dunlop, T. F. Havel, and D. G. Cory, Phys. Rev. A 60, 2777 (1999).
- Burrell et al. (2009) C. K. Burrell, J. Eisert, and T. J. Osborne, Phys. Rev. A 80, 052319 (2009).
- Bacon and Dam (2010) D. Bacon and W. V. Dam, Communications of the ACM 53, 84 (2010).