Qubit Lattice Algorithms based on the Schrodinger-Dirac representation of Maxwell Equations and their Extensions

Qubit lattice algorithms (QLA) were first being developed in the late 1990’s to solve the Schrodinger equation [1-3] using unitary collision and streaming operators acting on some qubit basis. QLA recovers the Schrodinger equation in the continuum limit to second order in the spatial lattice grid spacing. Because the lattice node qubits are entangled by the unitary collision operator (much like in the formation of Bell states), QLA is encodable onto a quantum computer with an expected exponential speed-up over a classical algorithm run on a supercomputer. Moreover, since QLA is extremely parallelizeable on a classical supercomputer, it provides an alternate algorithm for solving difficult problems in computational classical physics.

We then applied these QLA ideas to the study of the nonlinear Schrodinger equation (NLS) [4], by incorporating the cubic nonlinearity in the wave function, $|\psi|^{2}\psi$, as an external potential operator following the unitary collide-stream operator sequence on the qubits. While the inclusion of such nonlinear terms poses no problem for a hybrid classical-quantum computer, it remains a very important and difficult research topic for their implementation on a quantum computer. The accuracy of the QLA for NLS was tested for soliton-soliton collisions in long-term integration and compared to exact analytic solutions, and while the QLA is second order, it seemed to behave like a symplectic integrator. The QLA was then extended to the totally integrable vector Manakov solitons [5] to handle inelastic soliton scattering. The Manakov solitons are solutions to a coupled set of NLS equations.

Following these successful benchmarking simulations, we moved into QLA for two (2D) and three (3D) dimensional NLS equations - where now there are no exact solutions to these nonlinear equations. In the field of condensed matter, these higher dimension NLS equations are known as the Gross-Pitaevskii equations and give the mean field representation of the ground state wave function $\psi$ of a zero-temperature Bose-Einstein condensate (BEC). For scalar quantum turbulence in 3D we [6] observed a triple energy cascade on a $5732^{3}$ grid, with the low-k (”classical”) regime exhibiting a Kolmogorov $k^{-5/3}$ cascade in the compressible kinetic energy while the incompressible kinetic energy exhibited a long k-range of $k^{-3}$ spectrum. Similar results were found for both 2D and 3D scalar quantum [7-9], while results for spinor BECs can be found in [10-12]. A somewhat related, but significantly different, approach is that of the quantum lattice Boltzmann method [13-14].

Here we will discuss a QLA for the solution of Maxwell equations in a tensor dielectric medium [15-18], and present some simulations results of the scattering of a 1D electromagnetic pulse off 2D localized dielectric objects. This can be viewed as a precursor to examining the scattering of electromagnetic pulses off plasma blobs in the exterior region of a tokamak.

There has been much interest in rewriting the Maxwell equations in operator form and exploit its similarity to the Schrodinger-Dirac equation from the early 1930’s (e.g., see the references in [19]). For homogeneous media the qubit representation of the electric and magnetic fields, E, H, leads to a Dirac equation in a fully unitary representation. However, when the media becomes inhomogeneous, a Dyson map [20] is required to yield a unitary Schrodinger-Dirac equation for the evolution of the electromagnetic qubit field representation. In particular, one can use the fields $(n_{x}E_{x},n_{y}E_{y},n_{z}E_{z},B_{x},B_{y},B_{z})$, where $n_{i}$ is the refractive index in the $i^{th}$-direction.

A QLA is developed for this representation of the Maxwell equations in Sec. 3. This particular algorithm is a generalization of that used for the NLS equations. The initial value problem is then solved for the case of an electromagnetic pulse propagating in the $x$-direction and scattering from different 2D localized dielectric objects with refractive index $n(x,y)$ in Sec. 4. In particular, we have examined both polarizations of the pulse and $\nabla\cdot\mathbf{B}=0$. In Sec. 5 we consider the case in which the medium is dissipative. This brings in the field of open quantum systems and interactions with an environment. For illustration we consider a simplified cold electron-ion dissipative fluid model in an electromagnetic field. Kraus operators are determined by a multidimensional analog of the quantum amplitude damping channel. Some estimates on the quantum gates required are given as well as a quantum circuit diagram illustrating the implementation of Kraus operators on the open Schrodinger equation. We summarize our results in Sec. 6.

Finally, in this introduction, we quickly review the entanglement of qubits - and in particular for the 2-qubit Bell state [21]

The most general 1-qubit states are $\{a_{0}|0\rangle+a_{1}|1\rangle\}$ , and $\{b_{0}|0\rangle+b_{1}|1\rangle\}$ with normalization $|a_{0}|^{2}+|a_{1}|^{2}=1=|b_{0}|^{2}+|b_{1}|^{2}$. The tensor product of these two 1-qubits yields a space of the form

However the Bell state, Eq. (1) is not part of this tensor product space: to remove the $|01\rangle$ state from Eq. (2) either $a_{0}=0$ or $b_{1}=0$. This in turn would remove either the $|00\rangle$ or the $|11\rangle$ states, respectively. Now consider the unitary collision operator

acting on the subspace basis $\{|00\rangle,|11\rangle\}$. The choice of $\theta=\pi/4$ yields the Bell state $B_{+}$ - a maximally entangled state. It is the quantum entanglement of states that will give rise to exponential speed-up of a quantum algorithm. The QLA is a sequence of interleaved unitary collision-streaming operators that entangle the qubits and then spread that entanglement throughout the lattice.

Consider the subset of Maxwell equations

and treat $\nabla\cdot\mathbf{B}=0$ and $\nabla\cdot\mathbf{D}=0$ as initial constraints that remain satisfied in the continuum limit for all times. (This, of course, follows immediately from taking the divergence of Eq. (4) ).

However, if the medium is spatially inhomogeneous, then $\left[\mathbf{W}^{-1},\mathbf{M}\right]\neq 0$ and the evolution equation for the $\mathbf{u}$-field is not unitary.

QLA consists of a sequence of unitary collision and streaming operators on a 2D spatial lattice which will recover the continuum Maxwell equations, Eq. (14) to second order in the spatial grid size, $\delta$. In particular we need to have 6 qubits/lattice site to represent the fields components in Eq. (13). QLA permits us to handle the x- and y-dependence separately. Let us first consider the x-dependence and recover the $\partial q_{i}/\partial x$ - terms. From Eq. (14) we see coupling between $q_{1}\leftrightarrow q_{5}$, $q_{2}\leftrightarrow q_{4}$, Hence we introduce the local entangling collision operator

The collision angles $\theta_{1}$ and $\theta_{2}$ need to be chosen to recover the refractive index factors before the corresponding spatial derivatives,

The first of the unitary streaming operators will stream qubits $q_{1},q_{4}$ one lattice unit in either direction while leaving the other four qubits fixed : $S_{14}^{\pm}$. The other unitary streaming operator will act on qubits $q_{2},q_{5}$ : $S_{25}^{\pm}$. The final unitary collide-stream sequence, $\mathbf{U_{X}}$ in the x-direction that leads to a second order scheme in $\delta$ can be shown to be

It should be noted that if only applies the first 4 collide-stream sequence in Eq. (17) then the algorithm would only be first order accurate.

Similarly, to recover the $\partial q_{i}/\partial y$-terms one would collisionally entangle qubits $q_{0}\leftrightarrow q_{5}$, $q_{2}\leftrightarrow q_{3}$ with

with corresponding collision angles $\theta_{0}$ and $\theta_{2}$. $\theta_{2}$ is given in Eq. (16), and

The streaming operator $S^{y}_{03}$ will act on qubits $q_{0},q_{3}$ only and similarly for the operator $S^{y}_{25}$. The final unitary collide-stream second-order accurate or the y-direction for Maxwell equations is

We still need to recover the spatial derivatives on the refractive index components in Eq. (14). To obtain the $\partial n_{z}/\partial x$ and $\partial n_{y}/\partial x$ terms we introduce the (non-unitary) sparse potential matrix

with collision angles

while the corresponding (non-unitary) sparse potential matrix to recover the $\partial/\partial y$-derivatives in the refractive index components is

with collision angles

Thus the final discrete QLA, that models the 2D Maxwell equations, Eq. (14), to $O(\delta^{2})$, advances the lattice qubit-vector $\mathbf{Q}(t)$ to $\mathbf{Q}(t+\Delta t)$ is

provided we have diffusion ordering in the space-time lattice: i.e., $\Delta t=\delta^{2}$. It is this ordering that requires us to have the unitary collision angles to be $O(\delta)$, Eqs. (16) and (19), and the external potential angles $O(\delta^{2})$, Eqs. (22). We note that computationally QLA is more accurate if we employ the external potentials twice: once half-way through the collide-stream sequence and then at the end.

We now present detailed QLA simulations to the initial value problem of the scattering of a 1D electromagnetic pulse from a localized dielectric object. In particular we consider a 1D Gaussian pulse propagating in the $x$-direction towards a localized dielectric object of refractive index $n(x,y)$. The initial pulse has non-zero field components $E_{z},B_{y}$,, Fig. 1, and scatters from either a localized cylindrical dielectric, Fig. 2a, or a conic dielectric object, Fig 2b. These simulations were performed for $\delta=0.1$.

It should be noted that QLA is an initial value algorithm : the refractive index profiles are smooth (e.g., hyperbolic tangents for the dielectric cylinder with boundary layer thickness $\approx 10$ lattice units) and so no internal boundary conditions are imposed at any time in the simulation.

So far we have treated Maxwell equations as a closed system based on the energy conservation dictated from the Hermiticity and positive definiteness of the constitutive matrix $\mathbf{W}$, Eq. (7), since we have restricted ourselves to perfect materials. However, when we wish to consider actual materials there is dissipation. This immediately defeats any attempt to pursue a unitary representation in the original Hilbert space. The obvious question is : can we embed our dissipative system into a higher dimension closed Hilbert space, and thus recover unitary evolution in this new space and build an appropriate QLA that can be encoded onto quantum computers? To accomplish this, we resort to open quantum system theory [21] to describe classical dissipation as the observable result of interaction between our system of interest and its environment.

For a closed quantum system, the time evolution of a pure state $|\psi(t)\rangle$ is given by the unitary evolution from the Schrodinger equation : $|\psi(t)\rangle=U(t)|\psi(0)\rangle$ with $U=exp[-itH_{0}]$ unitary for the Hermitian Hamiltonian $H_{0}$. The evolution of the density matrix, $\rho=|\psi\rangle\langle\psi|$, is governed by the corresponding von Neumann equation : $\rho(t)=U(t)\rho(0)U^{\dagger}(t)$. The density matrix formulation is required when dealing with composite systems. Kraus realized that the density matrix retains its needed properties if one generalized its evolution operator to

where the only restriction on the set of so-called Kraus matrices $K_{k}$ is that the sum of $K_{k}^{\dagger}K_{k}$ is the identity matrix. The evolution of the density matrix, Eq. (32), is no longer unitary for $k\geq 2$.

The Kraus representation [21] is most useful when dealing with quantum noisy operations due to interaction with an environment. For those problems in which this noisy operation translates into a dissipative process, the Hamiltonian for the system in the Schrodinger representation has both a Hermitian part, $H_{0}$, and an anti-Hermitian part, $i\,H_{1}$, that models the dissipation. A simple but non-trivial example is the $1D$ Maxwell equations (without sources) for a homogeneous scalar medium with electrical losses,

with complex permittivity $\epsilon=\epsilon_{R}+i\epsilon_{I}.$ $\epsilon^{*}=\epsilon_{R}-i\epsilon_{I}.$ $\hat{p}_{x}=-i\partial_{x}$ is the momentum operator. Introducing the Dyson map $\bm{\rho}=diag(|\epsilon|/\sqrt{\epsilon_{R}},\sqrt{\mu_{0}})$ into Eq. (33) and after some algebraic manipulations the evolution equation can be written as

where the state vector $\mathbf{Q}=\bm{\rho}\mathbf{u}$, where $\mathbf{u}$ is defined in Eq. (6). $\delta=\epsilon_{I}/\epsilon_{R}$ is the loss angle, $v_{\delta}$ is the phase velocity $v_{\delta}=1/\sqrt{\epsilon_{R}\mu_{0}(1+\delta^{2})}$ and $\sigma_{x},\sigma_{y},\sigma_{z}$ are the Pauli matrices.

The Schrodinger-Dirac equations are the backbone of the work presented here on Maxwell equations both in lossless inhomogeneous and lossy dielectric media. In both cases, straightforward application of unitary algorithms fail : in the first case, somewhat surprisingly one finds that even though a Dyson map points to the required electromagnetic field variables in a tensor dielectric, its implementation has till now defied a fully unitary representation. Our current QLA approach requires some external non-unitary operators that recover the terms involving the spatial derivatives on the refractive indices of the medium. These sparse matrices can be modeled by the sum of a linear combination of unitaries (LCU), which can then be encoded onto a quantum computer [23, 24]. In the second case, handling dissipative systems immediately forces us to consider an open quantum system interacting with its environment. Typically this forces us into a density matrix formulation and a clever introduction of what is known as Kraus operators [21, 25]. The beauty of the Kraus representation is that even though the system of interest is interacting with the environment, the Kraus operators do not need detailed information on the environment.

We presented detailed 2D scattering of a 1D electromagnetic pulse off localized dielectric objects. QLA is an initial value scheme. No internal boundary conditions are imposed at the vacuum-dielectrix interface. For dielectrics will large spatial gradients in the refractive index, QLA simulations show strong internal reflection/transmission within the dielectric object. These lead to quite complex time evolution of wavefronts from the dielectric objects. On the other hand, for weak spatial gradients in the refractive index there are negligible reflections from the vacuum-dielectric interface. This is reminiscent of WKB-like effects in the ray tracing approximation.

In considering the dissipative counterpart, one must now include both the system and its environment in order to get a closed system with unitary representation. The Kraus operators are the most general scheme that will retain the properties of the density matrix in time. The probability of obtaining the desired non-unitary evolution of the open system after the measurement operator $\hat{P}_{0}=|0\rangle\langle 0|\otimes I_{d\times d}$ is

The form of the unitary operator $\hat{\mathcal{U}}_{diss}$, Eq. (42) implies that it can be decomposed into r two-level unitary rotations $\hat{R}_{y}(\theta_{i})$ with $\cos(\theta_{i}/2)=e^{-\gamma_{i}\delta t}.$ Then, the quantum circuit implementation of $\hat{\mathcal{U}}_{diss}$ requires $O(rlog_{2}^{2}d)$ CNOT - and $\hat{R}_{y}(\theta_{i})$ quantum gates, so that there is a improvement in the circuit depth of $O(r/d^{2})$. Our multidimensional amplitude damping channel approach is directly related to the Sz. Nagy dilation by a rotation. The Sz. Nagy dilation [26] is the minimal unitary dilation containing the original dissipative (non-unitary) system.

This research was partially supported by Department of Energy grants DE-SC0021647, DE-FG02-91ER-54109, DE-SC0021651, DE-SC0021857, and DE-SC0021653. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 - EUROfusion). Views and opinions expressed, however, are those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. E. K. is supported by the Basic Research Program, NTUA, PEVE. K.H is supported by the National Program for Controlled Thermonuclear Fusion, Hellenic Republic. 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.

[1] B. M. Boghosian & W. Taylor IV, ”Quantum lattice gas models for the many-body Schrodinger equation”, Internat. J. Modern Phys. C8, 705-716 (1997).

[2] [1] B. M. Boghosian & W. Taylor IV, ”Simulating quantum mechanics on a quantum computer”, Physica D 120, 30-42 (1998).

[3] J. Yepez & B.M. Boghosian, ”An efficient and accurate quantum lattice-gas model for the many-body Schroedinger wave equation,” Computer Physics Communications, 146, (2002) 280-294.

[4] G. Vahala, J. Yepez & L. Vahala, ”Quantum lattice gas representation of some classical solitons”, Phys. Lett A 310, 187-196 (2003)

[5] G. Vahala, L. Vahala & J. Yepez, ”Inelastic vector soliton collisions: a lattice-based quantum representation”, Pbhil. Trans: Math, Phys. and Eng. Sciences, Roy. Soc. 362, 1677-1690 (2004)

[6] J. Yepez, G. Vahala, L. Vahala & M. Soe, ”Superfluid turbulence from quantum Kelvin waves to classical Kolmogorov cascades”, Phys. Rev. Lett. 103, 084501 (2009)

[7] G. Vahala, J. Yepez, L. Vahala, M. Soe, B. Zhang & S. Ziegeler, ”Poincare recurrence and spectral cascades in three-dimensional quantum turbulence”, Phys. Rev. E 84, 046713 (2011)

[8] B. Zhang, G. Vahala, L. Vahala, &M. Soe, ”Unitary quantum lattice algorithm for two dimensional quantum turbulence”, Phys. Rev. E 84, 046701 (2011)

[9] G. Vahala, B. Zhang, J. Yepez, L. Vahala & M. Soe, ”Unitary qubit lattice gas representation of 2D and 3D quantum turbulence”, in Advanced Fluid Dynamics, ed. H. W. Oh, (InTech , 2012), Chpt 11., p. 239 - 272.

[10] L. Vahala, M. Soe, G. Vahala & J. Yepez, ”Unitary qubit lattice algorithms for spin-1 Bose Einstein condensates”, Rad. Effects Def. Solids 174, 46-55 (2019).

[11] G. Vahala, L. Vahala & M.Soe, ”Qubit unitary lattice algorithms for spin-2 Bose-Einstein Condensates,I - Theory and Pade Initial Conditions”, Rad. Effects Def. Solids 175, 102-112 (2020).

[12] G. Vahala, M. Soe & L. Vahala, ”Qubit unitary lattice algorithms for spin-2 Bose-Einstein Condensates,II - vortex reconnection simulations and non-Abelian vortices”, Rad. Effects Def. Solids 175, 113-119 (2020).

[13] S. Palpacelli & S. Succi ”The Quantum Lattice Boltzmann Equation: Recent Developments”, Comm. Computational Phys. 4, 980-1007 (2008).

[14] S. Succi, F. Fillion-Gourdeau, & S. Palpacelli, ”Quantum lattice Boltzmann is a quantum walk”, EPJ Quantum. Technol 2, 12 (2015).

[15] G. Vahala, L. Vahala, M. Soe, & A. K. Ram, ”Unitary quantum lattice simulations for Maxwell equations in vacuum and in dielectric media”, J. Plasma Phys. 86, 905860518 (2020).

[16] G. Vahala, J. Hawthorne, L. Vahala, A. K. Ram & M. Soe, ”Quantum lattice representation for the curl equations of Maxwell equations”, Rad. Effects and Defects in Solids 177, 85 (2022).

[17] G. Vahala, M. Soe, L.Vahala, A. Ram, E. Koukoutsis, & K. Hizanidis, ”Qubit Lattice Algorithm Simulations of Maxwell’s Equations for Scattering from Anisotropic Dielectric Objects”, e-print arXiv:2301.13601 (2023).

[18] A. Oganesov, G. Vahala, L. Vahala & M. Soe, ”Effect of Fourier transform on the streaming in quantum lattice gas algorithms”, Rad. Effects and Defects in Solids. 173, 169 (2018).

[19] S. A. Khan & R. Jagannathan, ”A new matrix representation of the Maxwell equations based on the Riemann-Silberstein-Weber vector for a linear inhomogeneous medium”, arXiv:2205.09907v2.

[20] E. Koukoutsis, K. Hizanidis, A. K. Ram & G. Vahala, ”Dyson maps and unitary evolution for Maxwell equations in tensor dielectric media”, Phys. Rev. A 107, 042215 (2023).

[21] M. A. Nielsen & I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2010 (10th Ed.)

[22] M. Suzuki, ”Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivatives with applications to many-body probloems, Commun. Math. Phys. 51,183 (1976)

[23] A. M. Childs & N. Wiebe, ”Hamiltonian Simulation Using Linear Combinations of Unitary Operations”, Quantum Inf. and Comp. 12 , 901 (2012).

[24] A. M. Childs, R. Kothari & R. D. Somma, ”Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision”, SIAM J. on Comp. 46, 1920 (2017).

[25] W. F. Stinespring, ”Positive Functions on $C^{*}$-Algebras”, Proc. Amer. Math. Soc. 6, 211 (1955).

[26] V. Paulsen, Completely Bounded Maps and Operator Algebras, Cambridge Univ. Press, 2003.