Qubit lattice algorithms

KEYWORDS:

• Unitary
• qubits
• entanglement
• quantum computing

1. Introduction

Qubit lattice[Q1] algorithms (QLA) (1–7) can be viewed as a cousin of the more familiar lattice Boltzmann scheme (LB) (8) to solve classical physics problems computationally in a more efficient manner than standard computational (CFD) codes. In particular LB algorithms for nonlinear fluid systems achieved important status because of the simplicity of the code and their inherent parallelization on classical supercomputers. In LB, the fluid equations are modeled by a truncated kinetic description in which the permitted phase space velocity is minimized subject to recovering the required fluid equations and retaining appropriate symmetries to leading order by Chapman-Enskog expansions. An immediate advantage of the LB approach is that the notorious computational headaches of resolving the nonlinear convective derivative for turbulent flows in the momentum equation, $\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v}$, is now replaced by a simple linear kinetic advection (streaming) term ${\xi }_{\mathbf{i}}\cdot \mathrm{\nabla }{f}_i(\mathbf{x}\mathbf{,}{\xi }_{\mathbf{i}}\mathbf{,}\mathbf{t}\mathbf{)}$. ${\xi }_{\mathbf{i}}$ define the chosen discrete phase space lattice velocities. The fluid nonlinearities are recovered by incorporating quadratic algebraic nonlinearities into the LB collision operator. The initial collision operator was of BGK form $(f_i-f_i^{\mathrm{eq}})/\tau$. The relaxation distribution function $f_i^{\mathrm{eq}}$ is a function of the macroscopic variables, and τ is a relaxation time. The very early LB algorithms were restricted to low Mach number and low Reynolds number flows. However a literal explosion of research in LB has moved the field into multispeed-lattices, multiple relaxation times, entropic LB, supersonic flows using adaptive space-time reference frame based on the actual local fluid velocity and temperature (9–11)

The QLA employs a non-commuting sequence of unitary collide-stream operators acting on a chosen qubit basis set as a mesoscopic model that recovers the equations of interest in the continuum limit. In the beginning, QLA's were introduced for the solution of the Schrodinger equation (1–3). For quantum entanglement one requires at least 2 qubits/lattice site. The four on-site basis kets are (1) $|0⟩\otimes |0⟩=|1000⟩,|0⟩\otimes |1⟩=|0100⟩,|1⟩\otimes |0⟩=|0010⟩,|1⟩\otimes |1⟩=|0001⟩.$(1) The unitary collision operator required to recover the kinetic energy operator in the Schrodinger equation is just the square-root-of-swap gate (2) $U_1=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& (1-i)/2& (1+i)/2& 0\\ 0& (1+i)/2& (1-i)/2& 0\\ 0& 0& 0& 1\end{array}\right].$(2) It is convenient to work with the $2\times 2$ sub-block, C, entangling the on-site qubits $\{q_1,q_2\}$. $C^2$ is just the standard swap gate on these 2 qubits, (3) $C^2\left[\begin{array}{c}{q}_1\\ q_2\end{array}\right]={\left[\begin{array}{cc}(1-i)/2& (1+i)/2\\ (1+i)/2& (1-i)/2\end{array}\right]}^2\left[\begin{array}{c}{q}_1\\ q_2\end{array}\right]=\left[\begin{array}{c}{q}_2\\ q_1\end{array}\right]$(3) while $C^4=\mathcal{I}$, the identity operator.

The unitary streaming operator $S_1$ simply shifts the amplitude of qubit $q_1$ at lattice site x to site $x+\delta$, while the transpose $S_1^T$ shifts $q_1$ to lattice site $x-\delta$. Similarly for the shift operator $S_2$ acting on the qubit $q_2$. The operators U and S do not commute.

If one considers the following sequence of non-commuting unitary interleaved collision-streaming operators for a time advancement of $\mathrm{\Delta }t$ for the qubits $q_1,q_2$: (4) $\left[\begin{array}{c}{q}_1(t+\mathrm{\Delta t})\\ q_2(t+\mathrm{\Delta t})\end{array}\right]=S_2^T.C.S_2.C.S_2^T.C.S_2.C\cdot S_1^T.C.S_1.C.S_1^T.C.S_1.C\left[\begin{array}{c}{q}_1(t)\\ q_2(t)\end{array}\right]$(4) one recovers the free Schrodinger equation for the x-dependent wave function ψ determined from the zero moment of our mesoscopic qubit representation $\psi (x,t)=q_1(x,t)+q_2(x,t)$. The error is $O({\delta }^2)$, provided we have diffusion ordering, with $\mathrm{\Delta }t=O({\delta }^2)$. If one had performed just two collide-stream operators on each qubit then the scheme would become just $O(\delta )$. By using the alternate-direction implicit (ADI) scheme for the y-direction and z-direction independently, we recover the full 3D free Schrodinger equation.

It is very interesting to realize that if the collide-stream operators commuted, then Equation (4(4) $\left[\begin{array}{c}{q}_1(t+\mathrm{\Delta t})\\ q_2(t+\mathrm{\Delta t})\end{array}\right]=S_2^T.C.S_2.C.S_2^T.C.S_2.C\cdot S_1^T.C.S_1.C.S_1^T.C.S_1.C\left[\begin{array}{c}{q}_1(t)\\ q_2(t)\end{array}\right]$(4) ) simply degenerates into the identity evolution: (5) $\left[\begin{array}{c}{q}_1(t+\mathrm{\Delta t})\\ q_2(t+\mathrm{\Delta t})\end{array}\right]=\left[\begin{array}{c}{q}_1(t)\\ q_2(t)\end{array}\right]$(5) The non-commutation of the collide-stream operators is essential in our QLA in order to recover the desired equations of interest.

The effect of the potential energy term in the Schrodinger equation is readily modeled in QLA by the introduction of an external collision operator. In this way one can consider the scattering of a wave packet from an external potential well, or the 1D Nonlinear Schrodinger equation (NLS), or the 1D Korteweg-de Vries (KdV) equation. These last 2 problems are examples of nonlinear physics – and this poses a stumbling block for immediate application onto a quantum computer which is based on linear operator theory alone. However, both 1D NLS and 1D KdV are exactly soluble nonlinear equations – with solitons as solutions. In particular exact soliton-soliton collision solution provided excellent benchmarking for our QLA on classical computers, on which there are no problems handling the nonlinear terms.

We have extended QLA to both 2D and 3D NLS – now non-integrable equations that can give a mean-field theory for the evolution of a ground-state Bose-Einstein condensate (BEC) wave function. Very high precision QLA simulations were performed on quantum turbulence in these systems and their further generalization to spinor BECs. We achieved these results because QLA is ideally parallelized on classical supercomputers, with no degradation in parallel performance as the number of cores are increased. In particular, on the IGM BlueGene Mira supercomputer at Argonne, using OpenMP on a grid of ${5120}^3$ we attained strong scaling of 94.1% as the number of cores were increased from a base of 65, 536 to 524, 288 and achieving 1.17 PetaFlops (runs performed in 2016).

In some sense, QLA can be viewed as an example of a quantum random walk process.

2. QLA for plasma physics

Recently we have turned our attention to QLA for electromagnetic wave propagation in plasmas. Our approach has been to consider Maxwell equations and progressively include plasma effects into the constitutive equations. Eventually this will move into including the plasma evolution equations for mass, momentum and energy along with the evolution equations for the Maxwell electromagnetic fields $\mathbf{E},\mathbf{H}$. With constitutive equations (in a coordinate representation in which the Hermitian tensor dielectric $\mathit{ϵ}$ is diagonal) for an inhomogeneous non-magnetic medium (6) $\mathbf{D}=\mathit{ϵ}\cdot \mathbf{E},\mathbf{B}={\mu }_0\mathbf{H}.$(6) the Maxwell equations can be written in a unitary representation on discerning the required Dyson map (12) on a 6-qubit representation of the electromagnetic field, (7) $\mathbf{U}={(n_x{E}_x,n_y{E}_y,n_z{E}_z,{\mu }_0^{1/2}\mathbf{H})}^T=\mathbf{Q}$(7) where the refractive indices $\{n_x,n_y,n_z\}=\{{ϵ}_x^{1/2},{ϵ}_y^{1/2},{ϵ}_z^{1/2}\}$.

For 2D x-y spatially dependent fields and media, the unitary Maxwell equations take the form (8) $\begin{array}{rl}\frac{\mathrm{\partial }{q}_0}{\mathrm{\partial }t}& =\frac{1}{n_x}\frac{\mathrm{\partial }{q}_5}{\mathrm{\partial }y},\frac{\mathrm{\partial }{q}_1}{\mathrm{\partial }t}=-\frac{1}{n_y}\frac{\mathrm{\partial }{q}_5}{\mathrm{\partial }x},\frac{\mathrm{\partial }{q}_2}{\mathrm{\partial }t}=\frac{1}{n_z}[\frac{\mathrm{\partial }{q}_4}{\mathrm{\partial }x}-\frac{\mathrm{\partial }{q}_3}{\mathrm{\partial }y}]\\ \frac{\mathrm{\partial }{q}_3}{\mathrm{\partial }t}& =-\frac{\mathrm{\partial }(q_2/n_z)}{\mathrm{\partial }y},\frac{\mathrm{\partial }{q}_4}{\mathrm{\partial }t}=\frac{\mathrm{\partial }(q_2/n_z)}{\mathrm{\partial }x},\frac{\mathrm{\partial }{q}_5}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }(q_1/n_y)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }(q_0/n_x)}{\mathrm{\partial }y}\end{array}$(8) where $\mathbf{Q}=\{q_0,q_1,q_2,q_3,q_4,q_5{\}}^T$. There are considerable number of internal reflections/transmissions into the vacuum region since there is a sharp (but continuous) spatial gradient in the refractive index near the boundary of the vacuum-dielectric region. It should be noted that QLA is an initial value algorithm. There are no internal boundary conditions imposed at the vacuum-dielectric interface. There are no Fresnel-like boundary conditions imposed, but the reflection and transmission wave fronts are self-consistently generated in QLA.

2.1. QLA for dispersive media

The results presented here assumed the instantaneous (so-called ‘optical’) response of the medium to the given electromagnetic field $\mathbf{u}=\{\mathbf{E},\mathbf{H}{\}}^T$. Dispersion can be included in the constitutive equations by introducing the susceptibility kernel $\hat{G}$ (13) (9) $\mathbf{d}(\mathbf{r},t)=\hat{W}(\mathbf{r})\mathbf{u}(\mathbf{r},t)+{\int }_0^t\hat{G}(\mathbf{r},t-\tau )\mathbf{u}(\mathbf{r},\tau ).$(9) Here $\mathbf{d}=\{\mathbf{D},\mathbf{B}{\}}^T$, and $\hat{W}=\mathrm{diag}(\mathit{ϵ}(\mathbf{r}),{\mu }_0)$. It is convenient to move to the high frequency limit in which the diagonal optical response reduces to that in a vacuum ${\hat{W}}_{\mathrm{vac}}=\mathrm{diag}({\mathit{ϵ}}_0,{\mu }_0)$. The polarization of the medium is now incorporated into the integral kernel $\hat{G}$. The QLA will be unitary, but would now need to introduce extra qubits at each lattice site for the dynamics of the polarization and the handling of integrals of the form (10) ${\int }_0^t{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{e^{-\mathrm{i\omega }(t-\tau )}}{{\omega }_e^2-{\omega }^2}\mathbf{E}(\mathbf{r},\tau )\mathrm{d}\omega \mathrm{d}\tau$(10) where ${\omega }_e$ is a characteristic frequency of the medium.

2.2. QLA for dissipative media

The inclusion of dissipation into the plasma equations requires a fundamental change in the QLA representation. One method is to employ the techniques of quantum information science (QIS) when dealing with quantum system that interact with the environment. Such quantum systems are termed ‘open’ systems since there is an energy exchange between it and the environment. However by considering the extended system of the dissipative quantum system and its environment one again recovers a closed system with energy conservation. This will then permit a unitary representation for the extended system.

The non-unitary time evolution of the open system takes the form (13) (11) $i\frac{\mathrm{\partial }\mathit{\psi }}{\mathrm{\partial }t}=[{\hat{D}}_0-i{\hat{D}}_{\mathrm{diss}}]\mathit{\psi }$(11) where ${\hat{D}}_0$ and ${\hat{D}}_{\mathrm{diss}}$ are Hermitian. $i{\hat{D}}_{\mathrm{diss}}$ is anti-Hermitian and represents the dissipation in the system. Thus (12) $\mathit{\psi }(\mathrm{\Delta }t)=\hat{\mathcal{U}}(\mathrm{\Delta t})\mathit{\psi }(0)=\exp -i\mathrm{\Delta t}[{\hat{D}}_0-i{\hat{D}}_{\mathrm{diss}}]\mathit{\psi }(0).$(12) A first crude approximation is to apply the zeroth order-Suzuki-Trotter decomposition to Equation (12(12) $\mathit{\psi }(\mathrm{\Delta }t)=\hat{\mathcal{U}}(\mathrm{\Delta t})\mathit{\psi }(0)=\exp -i\mathrm{\Delta t}[{\hat{D}}_0-i{\hat{D}}_{\mathrm{diss}}]\mathit{\psi }(0).$(12) ) (13) $\exp -i\mathrm{\Delta t}[{\hat{D}}_0-i{\hat{D}}_{\mathrm{diss}}]=e^{-i\mathrm{\Delta t}{\hat{D}}_0}{e}^{-\mathrm{\Delta }t{\hat{D}}_{\mathrm{diss}}}+\mathit{O}\mathit{(}\mathit{\Delta }{\mathit{t}}^{\mathit{2}}\mathit{)}\mathit{,}$(13) which neglects the effect of the non-commutation of ${\hat{D}}_0$ and ${\hat{D}}_{\mathrm{diss}}$. This approximation, extremely common in QIS, permits the treatment of the non-unitary operator $\exp [-\mathrm{\Delta t}{\hat{D}}_{\mathrm{diss}}]$ separately from the unitary operator $\exp [-i\mathrm{\Delta t}{\hat{D}}_0]$.

In QIS, one technique to handle the non-unitary operators is the introduction of Kraus operators ${\hat{K}}_{\mu }$ for the evolution of the open system density matrix ${\rho }_S$ (14) ${\rho }_S(t)=\sum _{\mu }{\hat{K}}_{\mu }{\rho }_S(0){\hat{K}}_{\mu }^{†}.$(14) The individual Kraus operators are arbitrary, but subject to the sum constraint be unitary: (15) $\sum _{\mu }{\hat{K}}_{\mu }^{†}{\hat{K}}_{\mu }=\mathcal{I}.$(15) By defining the Kraus operator ${\hat{K}}_0$ (16) ${\hat{K}}_0=\exp \{-\mathrm{\Delta t}{\hat{D}}_{\mathrm{diss}}\}$(16) and introducing a second Kraus operator ${\hat{K}}_1$ such that ${\hat{K}}_1^{†}{\hat{K}}_1=\mathcal{I}-{\hat{K}}_0^{†}{\hat{K}}_0$ one can extend to a higher dimensional Hilbert space (which in QIS is called a dilation) in which the dilation operator ${\hat{\mathcal{U}}}_{\mathrm{diss}}$ is unitary (17) ${\hat{\mathcal{U}}}_{\mathrm{diss}}=\left[\begin{array}{cc}{\hat{K}}_0& -{\hat{K}}_1^{†}\\ {\hat{K}}_1& \hat{\mathcal{X}}{\hat{K}}_0\hat{\mathcal{X}}\end{array}\right].$(17) The operator $\hat{\mathcal{X}}$, is an appropriate extension of the Pauli $\hat{X}$ operator to the appropriate dimensionality.

2.3. QLA for Plasma-Maxwell system

To develop a QLA to simulate directly the nonlinear two fluid plasma equations for mass, momentum and energy is not at all straightforward. Extending the work of Meng and Yang (14), one may utilize the Madelung transformation on the 3D NLS BEC equation (known in condensed matter at the Gross-Pitaevskii equation) for the scalar wavefunction ψ: (18) $\psi =\sqrt{\rho }{e}^{\mathrm{i\varphi }}.$(18) One can immediately identify $\rho =\sqrt{{\psi }^{\ast }\psi }$ as the fluid density, and $\mathbf{v}=-\mathrm{\nabla }\varphi$ as the fluid velocity. The scalar BEC equation for ψ can be rewritten in the form of fluid conservation equation of mass, momentum and energy, not unlike the compressible Navier-Stokes system. The major difference between the quantum NLS equation and the classical Navier-Stokes system is the addition of an extra quantum pressure term in the momentum equation and the fact that for scalar NLS quantum vortices are singular, with $\rho \to 0$ at the vortex core. A classical vortex, on the other hand is non-singular. To remove these quantum irregularities, Meng & Yang achieve extra degrees of freedom by moving to a quaternion representation of the qubit fields whose zeroth moment yields the wave function, and thereby remove the quantum pressure term and singular vortices.

It seems very worthwhile to extend this approach to spinor BECs, with the Hamiltonian augmented to incorporate the electromagnetic field. We already have very well benchmarked and parallelized QLA's for spinor BECs and separate QLA's for Maxwell equations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was partially supported by Department of Energy grants DE-SC0021647, DE-FG0291ER-54109, DE-SC0021651, DE-SC0021857 and DE-SC0021653. This work has been carried out partially within the framework of the EUROfusion Consortium. E.K has received funding from the Euratom research and training program WPEDU under grant agreement no. 101052200 as well as from the National Program for Controlled Thermonuclear Fusion, Hellenic Republic. K.H is supported by the National Program for Controlled Thermonuclear Fusion, Hellenic Republic. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0020430.

Notes on contributors

George Vahala

George Vahala is a professor of Physics at WIlliam & Mary, Lattice algorithms, and applications to plasma physics.

Efstratios Koukoutsis

Efstartios Koukoutsis is a PhD student in National Technical University of Athens, with research in developing quantum algorithms for dispersive dissipative plasma models.

Min Soe

Min Soe is a professor at Rogers State University - parallelization of algorithms and 3D graphics.

Kyriakos Hizanidis

Kyriakos Hizanidis is a professor of Electrical and Computer Engineering at National Technical University of Athens - Clifford algebra representation of Maxwell equations, quantum computing algorithms.

Linda Vahala

Linda Vahala is an associate professor of Electrical & Computer Engineering, Old Dominion University - quantum and qubit lattice algoritihms, symmetry properties.

Abhay K. Ram

Abhay K. Ram Plasma Science and Fusion Center, MIT - wave propagation in plasmas, unitary algorithms.