Qubit Lattice Algorithms Based on the Schrödinger-Dirac Representation of Maxwell Equations and Their Extensions

Abstract

It is well known that Maxwell equations can be expressed in a unitary Schrodinger-Dirac representation for homogeneous media. However, difficulties arise when considering inhomogeneous media. A Dyson map points to a unitary field qubit basis, but the standard qubit lattice algorithm of interleaved unitary collision-stream operators must be augmented by some sparse non-unitary potential operators that recover the derivatives on the refractive indices. The effect of the steepness of these derivatives on two-dimensional scattering is examined with simulations showing quite complex wavefronts emitted due to transmissions/reflections within the dielectric objects. Maxwell equations are extended to handle dissipation using Kraus operators. Then, our theoretical algorithms are extended to these open quantum systems. A quantum circuit diagram is presented as well as estimates on the required number of quantum gates for implementation on a quantum computer.

1. Introduction

Qubit lattice algorithms (QLA) were first being developed in the late 1990s to solve the Schrodinger equation [1, 2, 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 parallelizable 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, ${\left|\psi \right|}^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 dimensions 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, 8, 9], while results for spinor BECs can be found in Refs. [10, 11, 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, 16, 17, 18] and present some simulation 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 their similarity to the Schrodinger-Dirac equation from the early 1930s (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 $\left(n_x{E}_x,n_y{E}_y,n_z{E}_z,B_x,B_y,B_z\right)$, where $n_i$ is the refractive index in the $i^{\mathit{th}}$-direction.

A QLA is developed for this representation of the Maxwell equations in Section 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\left(x,y\right)$ in Section 4. In particular, we have examined both polarizations of the pulse and $\nabla \cdot \mathbf{B}=0$. In Section 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 Section 6.

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

$$B_{+}=\frac{1}{\sqrt{2}}\left(|00⟩+|11⟩\right).$$E1

The most general 1-qubit states are $\left\{a_0|0⟩+a_1|1⟩\right\}$, and $\left\{b_0|0⟩+b_1|1⟩\right\}$ with normalization ${\left|a_0\right|}^2+{\left|a_1\right|}^2=1={\left|b_0\right|}^2+{\left|b_1\right|}^2$. The tensor product of these two 1-qubit yields a space of the form

$$a_0{b}_0\mid 00⟩+a_0{b}_1\mid 01⟩+a_1{b}_0\mid 10⟩+a_1{b}_1\mid 11⟩.$$E2

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

$$C=\left[\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta \end{array}\right]$$E3

acting on the subspace basis $\left\{|00⟩,|11⟩\right\}$. 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 the 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.

2. The Dyson map and the generation of a unitary evolution equation for Maxwell equations

Consider the subset of Maxwell equations

$$\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t},\nabla \times \mathbf{H}=\frac{\partial \mathbf{D}}{\partial t}$$E4

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)).

For lossless media, the electric and magnetic fields satisfy the constitutive relations for a tensor dielectric nonmagnetic medium

$$\mathbf{D}=\mathit{\epsilon }\cdot \mathbf{E},\mathbf{B}={\mu }_0\mathbf{H}.$$E5

For Hermitian $\mathit{\epsilon }$, one can transform to a coordinate system in which $\mathit{\epsilon }$ is diagonal. Eq. (5) can be rewritten in matrix form

$$\mathbf{d}=\mathbf{W}\mathbf{u},\text{with}\mathbf{d}\doteq {\left(\mathbf{D},\mathbf{B}\right)}^{\mathbf{T}},\mathbf{u}\doteq {\left(\mathbf{E},\mathbf{H}\right)}^{\mathbf{T}}$$E6

where W is a $6\times 6$ Hermitian block diagonal constitutive matrix

$$\mathbf{W}=\left[\begin{array}{cc}\mathit{\epsilon }{\mathbb{I}}_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& {\mu }_0{\mathbb{I}}_{3\times 3}\end{array}\right].$$E7

${\mathbb{I}}_{3\times 3}$ is the $3\times 3$ identity matrix, and the superscript $\mathbf{T}$ in Eq. (6) is the transpose. In matrix form, the Maxwell equations, Eq. (4) become

$$i\frac{\partial \mathbf{d}}{\partial t}=\mathbf{M}\mathbf{u}$$E8

where under standard boundary conditions, the curl-matrix operator $\mathbf{M}$ is Hermitian:

$$\mathbf{M}=\left[\begin{array}{cc}{0}_{3\times 3}& i\nabla \times \\ -i\nabla \times & 0_{3\times 3}\end{array}\right].$$E9

From Eq. (6), since ${\mathbf{W}}^{-\mathbf{1}}$ exists, $\mathbf{u}={\mathbf{W}}^{-1}\mathbf{d}$, so that Eq. (8) can be written

$$i\frac{\partial \mathbf{u}}{\partial t}={\mathbf{W}}^{-\mathbf{1}}\mathbf{Mu}$$E10

If the medium is homogeneous, then ${\mathbf{W}}^{-1}$ is constant and will commute with the curl-operator $\mathbf{M}$. Under these conditions, the product ${\mathbf{W}}^{-1}\mathbf{M}$ is Hermitian and Eq. (10) gives unitary evolution for $\mathbf{u}={\left(\mathbf{E},\mathbf{H}\right)}^T$.

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

2.1 Dyson map

To determine a unitary evolution of the electromagnetic fields in an inhomogeneous dielectric medium, it [20] has been shown that there exists a Dyson map $\mathit{\rho }$: $\mathbf{u}\to \mathbf{Q}$ such that in the new field variables $\mathbf{Q}$ the resulting evolution equation will be unitary. For the Maxwell equations consider

$$\mathbf{Q}=\mathit{\rho }\mathbf{u}={\mathbf{W}}^{1/2}\mathbf{u}.$$E11

For time-independent media, the evolution equation for the new fields $\mathbf{Q}$ is

$$i\mathit{\rho }\frac{\partial \mathbf{u}}{\partial t}=\mathit{\rho }{\mathbf{W}}^{-\mathbf{1}}\mathbf{M}{\mathit{\rho }}^{-1}\mathit{\rho }\mathbf{u}\Rightarrow i\frac{\partial \mathbf{Q}}{\partial t}=\mathit{\rho }{\mathbf{W}}^{-\mathbf{1}}\mathbf{M}{\mathit{\rho }}^{-1}\mathbf{Q}$$E12

and is indeed unitary. Explicitly, the new fields, Eq. (11) and the $\mathit{\rho }$ are

$$\mathbf{Q}=\left[\begin{array}{l}{q}_0\\ q_1\\ q_2\\ q_3\\ q_4\\ q_5\end{array}\right]=\left[\begin{array}{l}{n}_x{E}_x\\ n_y{E}_y\\ n_z{E}_z\\ {\mu }_0^{1/2}{H}_x\\ {\mu }_0^{1/2}{H}_y\\ {\mu }_0^{1/2}{H}_z\end{array}\right],\mathit{\rho }=\left[\begin{array}{cccccc}{n}_x& 0& 0& 0& 0& 0\\ 0& n_y& 0& 0& 0& 0\\ 0& 0& n_z& 0& 0& 0\\ 0& 0& 0& {\mu }_0^{1/2}& 0& 0\\ 0& 0& 0& 0& {\mu }_0^{1/2}& 0\\ 0& 0& 0& 0& 0& {\mu }_0^{1/2}\end{array}\right]$$E13

The refractive index $n_i=\sqrt{{\epsilon }_i}$. Typically, we will use units where ${\mu }_0=1$. In component form, Maxwell equations for fields and with constitutive matrix restricted to spatially 2D $\left(x,y\right)$ dependence, Eq. (12) reduces to

$$\begin{array}{l}\frac{\partial q_0}{\partial t}=\frac{1}{n_x}\frac{\partial q_5}{\partial y},\frac{\partial q_1}{\partial t}=-\frac{1}{n_y}\frac{\partial q_5}{\partial x},\frac{\partial q_2}{\partial t}=\frac{1}{n_z}\left[\frac{\partial q_4}{\partial x}-\frac{\partial q_3}{\partial y}\right]\\ \frac{\partial q_3}{\partial t}=-\frac{\partial \left(q_2/n_z\right)}{\partial y},\frac{\partial q_4}{\partial t}=\frac{\partial \left(q_2/n_z\right)}{\partial x},\frac{\partial q_5}{\partial t}=-\frac{\partial \left(q_1/n_y\right)}{\partial x}+\frac{\partial \left(q_0/n_x\right)}{\partial y}\end{array}$$E14

3. A qubit lattice representation for 2D tensor dielectric media

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 sites to represent the field 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

$$C_X=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& \mathit{cos}{\theta }_1& 0& 0& 0& -\mathit{sin}{\theta }_1\\ 0& 0& \mathit{cos}{\theta }_2& 0& -\mathit{sin}{\theta }_2& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& \mathit{sin}{\theta }_2& 0& \mathit{cos}{\theta }_2& 0\\ 0& \mathit{sin}{\theta }_1& 0& 0& 0& \mathit{cos}{\theta }_1\end{array}\right]$$E15

The collision angles ${\theta }_1$ and ${\theta }_2$ need to be chosen to recover the refractive index factors before the corresponding spatial derivatives,

$${\theta }_1=\frac{\delta }{4n_y},{\theta }_2=\frac{\delta }{4n_z}.$$E16

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}}_{\mathbf{X}}$ in the x-direction that leads to a second-order scheme in $\delta$ can be shown to be

$${\mathbf{U}}_{\mathbf{X}}=S_{25}^{+x}.C_X^{†}.S_{25}^{-x}.C_X.S_{14}^{-x}.C_X^{†}.S_{14}^{+x}.C_X.S_{25}^{-x}.C_X.S_{25}^{+x}.C_X^{†}.S_{14}^{+x}.C_X.S_{14}^{-x}.C_X^{†}.$$E17

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

$$C_Y=\left[\begin{array}{cccccc}\mathit{cos}{\theta }_0& 0& 0& 0& 0& \mathit{sin}{\theta }_0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& \mathit{cos}{\theta }_2& \mathit{sin}{\theta }_2& 0& 0\\ 0& 0& -\mathit{sin}{\theta }_2& \mathit{cos}{\theta }_2& 0& 0\\ 0& 0& 0& 0& 1& 0\\ -\mathit{sin}{\theta }_0& 0& 0& 0& 0& \mathit{cos}{\theta }_0\end{array}\right],$$E18

with corresponding collision angles ${\theta }_0$ and ${\theta }_2$. ${\theta }_2$ is given in Eq. (16), and

$${\theta }_0=\frac{\delta }{4n_x}.$$E19

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

$${\mathbf{U}}_{\mathbf{Y}}=S_{25}^{+y}.C_Y^{†}.S_{25}^{-y}.C_Y.S_{03}^{-y}.C_Y^{†}.S_{03}^{+y}.C_Y.S_{25}^{-y}.C_Y.S_{25}^{+y}.C_Y^{†}.S_{03}^{+y}.C_Y.S_{03}^{-y}.C_Y^{†}$$E20

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

$$V_X=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& -\mathit{sin}{\beta }_2& 0& \mathit{cos}{\beta }_2& 0\\ 0& \mathit{sin}{\beta }_0& 0& 0& 0& \mathit{cos}{\beta }_0\end{array}\right]$$E21

with collision angles

$${\beta }_0={\delta }^2\frac{\partial n_y/\partial x}{n_y^2},{\beta }_2={\delta }^2\frac{\partial n_z/\partial x}{n_z^2},$$E22

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

$$V_Y=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& o\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& \text{cos}{\beta }_3& \text{sin}{\beta }_3& 0& 0\\ 0& 0& 0& 0& 1& 0\\ -\mathit{sin}{\beta }_1& 0& 0& 0& 0& \mathit{cos}{\beta }_1\end{array}\right]$$E23

with collision angles

$${\beta }_1={\delta }^2\frac{\partial n_x/\partial y}{n_x^2},{\beta }_3={\delta }^2\frac{\partial n_z/\partial y}{n_z^2}.$$E24

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

$$\mathbf{Q}\left(t+\mathrm{\Delta }t\right)=V_Y.V_X.{\mathbf{U}}_{\mathbf{Y}}.{\mathbf{U}}_{\mathbf{X}}.\mathbf{Q}\left(t\right)$$E25

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

3.1 Non-unitary external potential operators

Recently, considerable effort has been expended into developing more efficient approximation for handling the evolution operator of a complex Hamiltonian system than the standard Suzuki-Trotter expansion Eqs. (21) and (23) [22]. In particular, the idea [23, 24] has been floated of approximating the full unitary operator by a sum of unitary operators. The actual implementation onto a quantum computer we will leave to another paper, as one of the outcomes of QLA discussed here will be a quantum-inspired highly efficient classical supercomputer algorithm. Moreover, its encoding onto a quantum computer will require error-correcting qubits with long coherence times, something currently out of reach in the noisy qubit regime we are in. Here, we will show the 4 unitary operators needed whose sum yields the sparse non-unitary potential operator $V_X$, Eq. (21). Letting ${\mathbb{I}}_6$ be the $6\times 6$ identity matrix, then it is easily verified that

$$V_X=\frac{1}{2}\sum _{i=1}^4{\mathit{LCU}}_i$$E26

where the first two unitaries are diagonal

$${\mathit{LCU}}_1={\mathbb{I}}_6,{\mathit{LCU}}_2=\mathit{diag}\left(-\mathrm{1,1,1},-1,-1,-1\right)$$E27

and the remaining two unitaries are

$${\mathit{LCU}}_3=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& o\\ 0& \text{cos}{\beta }_0& 0& 0& 0& -sin{\beta }_0\\ 0& 0& \text{cos}{\beta }_2& 0& sin{\beta }_2& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& -sin{\beta }_2& 0& \text{cos}{\beta }_2& 0\\ 0& sin{\beta }_0& 0& 0& 0& \text{cos}{\beta }_0\end{array}\right],$$E28

and

$${\mathit{LCU}}_4=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& o\\ 0& -cos{\beta }_0& 0& 0& 0& sin{\beta }_0\\ 0& 0& -cos{\beta }_2& 0& -sin{\beta }_2& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& -sin{\beta }_2& 0& cos{\beta }_2& 0\\ 0& sin{\beta }_0& 0& 0& 0& cos{\beta }_0\end{array}\right].$$E29

3.2 Conservation of energy

In a fully unitary representation, the norm of Q is a constant of the motion. This is simply the conservation of energy in the electromagnetic field. For fields being a function of $\left(x,y\right)$, we have from Eqs. (11)–(13)

$$\mathrm{ℰ}\left(t\right)=\frac{1}{L^2}{\int }_0^L\mathit{\text{dxdy}}\mathbf{Q}\cdot \mathbf{Q}=\frac{1}{L^2}{\int }_0^L\mathit{\text{dxdy}}\left[{\epsilon }_x{E}_x^2+{\epsilon }_y{E}_y^2+{\epsilon }_z{E}_z^2+{\mu }_o\left(B_x^2+B_y^2+B_z^2\right)\right]$$E30

where the (diagonal) tensor dielectric ${\epsilon }_x=n_x^2,\dots$ and we restrict ourselves to nonmagnetic materials, for simplicity.

4. 2D numerical simulations from the QLA for electromagnetic scattering from 2D dielectric objects

We now present detailed QLA simulations of 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 toward a localized dielectric object of refractive index $n\left(x,y\right)$. The initial pulse has nonzero field components $E_z,B_y$, Figure 1, and scatters from either a localized cylindrical dielectric, Figure 2a, or a conic dielectric object, Figure 2b. These simulations were performed for $\delta =0.1$.

Figure 1.

A 1D electromagnetic pulse with initial fields $-E_z\left(x,t=\right.$ 0)$,B_y\left(x,t=\right.$ 0). 2D simulation grid $L\times L$ with $L=$ 8192. Pulse full-width (in lattice units) $\approx$ 200. Since the Maxwell equations are linear and homogenous, the initial amplitude of the fields is arbitrary.

Figure 2.

(a) The dielectric cylinder, diameter $\approx$ 200, has rapidly increasing boundary dielectric from vacuum $n=1$ to $n_{\mathit{max}}=3$_, whereas (b) the conic dielectric, base $\approx$ 240, has smoothly increasing dielectric from vacuum to conic peak of $n_{\mathit{max}}=$ 3.

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.

4.1 Effects of broken symmetry

When the 1D pulse scatters off the dielectric object with refractive index $n\left(x,y\right)$ the initial electric field spatial dependence $E_z\left(x\right)$ now becomes a function $E_z\left(x,y\right)$, while the initial magnetic field $B_y\left(x\right)$ will become a function $B_y\left(x,y\right)$. Now the scattered field has $\partial B_y/\partial y\ne 0$. Thus, for $\nabla \cdot \mathbf{B}=0$, the scattered field must develop an appropriate $B_x\left(x,y\right)$. This is seen in our QLA simulations, even though the explicit discrete collide-stream algorithm only models asymptotically the Maxwell subset, Eq. (4). $\nabla \cdot \mathbf{D}=0$ and this is exactly conserved in our QLA simulation.

Similarly, for initial $E_y\left(x\right)$ polarization. In this case, the scattered magnetic field $B_z\left(x,y\right)$ satisfies, for our 2D scattering in the x-y plane, $\nabla \cdot \mathbf{B}=0$ exactly and no other magnetic field components are generated. Our discrete QLA recovers this $\nabla \cdot \mathbf{B}=0$ exactly. However, in an attempt to preserve $\nabla \cdot \mathbf{D}=0$, the QLA will generate a nonzero $E_x\left(x,y\right)$ field.

4.2 Scattering from localized 2D dielectric objects with refractive index $n\left(x,y\right)$

Consider a 1D pulse with polarization $E_z\left(x,t\right)$ propagating in a vacuum toward a 2D dielectric scatterer. In Figures 3–6, we consider the time evolution of the resultant scattered $E_z$-field. The initial pulse is followed for a short time while it is propagating in the vacuum to verify that the QLA correctly determines its motion. As part of the pulse interacts with the dielectric object, the pulse speed within the dielectric itself is decreased by the inverse of the refractive index profiles, $n\left(x,y\right)$. The remainder of the 1D pulse propagates undisturbed since it is still propagating within the vacuum.

Figure 3.

The scattered $E_z$ field after 15,000 iterations (i.e., $t=$ 15k). (a) There is an internal reflection at the back of the cylindrical dielectric.

Figure 4.

The scattered $E_z$ field at $t=$ 23.4k. (a) a reflected circular wavefront occurs as that part of the pulse reaches the back-end of the cylindrical dielectric, along with the initial reflected circular wavefront with its $\pi$-changed phase at the front of the vacuum-cylinder boundary. (b) for the conic dielectric, there is an internal reflection from the apex of the cone’s $n_{\mathit{max}}$, which then propagates out of the weakly varying cone edges.

Figure 5.

The scattered $E_z$ field at $t=$ 31.2k. (a) There are multiple reflections/transmissions within the boundaries of the dielectric cylinder. (b) There is only one major reflection from the apex of the cone and which then propagates readily out from the cone edge.

Figure 6.

The scattered $E_z$ field at $t=$ 49.2k. (a) There are multiple reflections/transmissions within the boundaries of the dielectric cylinder. (b) There is only one major reflection from the apex of the cone and which then propagates readily out from the cone edge.

One sees in Figure 3a, a circular-like wavefront reflecting back into the vacuum, with its $E_z$ field $\pi$ out of phase as the reflection is occurring from a low to higher refractive index around the vacuum-cylinder interface. One does not find such a reflected wavefront when the pulse interacts with the conic dielectric, Figure 3b.

At $t=23.4k$, there is a major wavefront emanating from the back of the cylindrical dielectric, Figure 4a. For the conic dielectric, there is a major wavefront reflected from the apex of the conic dielectric, and this propagates out of the cone with a little reflection, Figure 4b.

One clearly sees at $t=31.2k$ that more $E_z$ wavefronts are being created because of the large refractive index gradients at the vacuum-cylinder dielectric boundary, while such gradients are missing from the vacuum-cone interface which leads to no new wavefronts in the scattering of the dielectric cone, Figure 5.

At $t=49.2k$, the complex $E_z$ wavefronts are due to repeated reflections and transmissions from the cylinder dielectric, Figure 5a. However, because of the slowly changing boundaries of the dielectric cone there are no more reflections and one sees only the outgoing wavefront from the pulses’ interaction with the region around the $n_{\mathit{max}}$ of the cone. Since the pulse reaches the apex of the cone before the corresponding pulse hits the backend of the dielectric cylinder, the conic wavefront is further advanced than that of the cylindrical wavefront, Figure 5b.

4.2.1 Auxiliary fields and $\nabla \cdot \mathbf{B}$

For incident $E_z$ polarization and with 2D refractive index $n\left(x,y\right)$, the scattered electromagnetic fields will need to generate a $B_x$ field in order to have $\nabla \cdot \mathbf{B}=0$. In Figure 6a and b, we plot the self-generated $B_x\left(x,y\right)$ field at $t=23.4k$ and $t=49.2k$ for scattering from the dielectric cylinder. It is also found that $\mid \nabla \cdot \mathbf{B}\mid /B_0$ is typically zero everywhere in the spatial lattice except for a very localized region around the vacuum-dielectric boundary layer where the normalized $\mathit{max}\left(|\nabla \cdot \mathbf{B}|\right)/B_0$ reaches around 0.01 at very few isolated grid points (Figure 7).

Figure 7.

The self-consistently generated $B_x\left(x,y\right)$-field after the 1D incident pulse with $B_y=B_y\left(x\right)$ scatters from a local dielectric with refractive index $n\left(x,y\right)$ at times: (a) $t=$ 23.4k, corresponding to Figure 4a for $E_z$, and (b) $t=$ 49.2k, corresponding to Figure 5a for $E_z$.

4.3 Time dependence of $\mathrm{ℰ}\left(t\right)$ on perturbation parameter $\delta$

The discrete total electromagnetic energy $\mathrm{ℰ}\left(t\right)$. Eq. (30), is not constant since our current QLA is not totally unitary. However, the variations in $\mathrm{ℰ}$ decrease significantly as $\delta \to 0$. $\delta$ is a measure of the discrete lattice spacing. The maximal variations occur shortly after the 1D pulse scatters from the 2D dielectric object. For $\delta =0.1$, this occurs around $t=15k$, with variations in the 5th decimal, Eq. (31). However, when one reduces $\delta$ by a factor of $10$ on the same lattice grid, then one recovers the same physics a factor of $10$ later in time since $\delta$ controls the speed of propagation in the vacuum. Thus, the wallclock time of a QLA run is also increased by this factor of $10$. We find, for $\delta =0.01$, that the largest deviation in the total electromagnetic energy is now in the 8th decimal, Eq. (31).

$$\mathit{\delta }\mathbf{=}\mathbf{0.1}:\begin{array}{ccc}\mathit{\text{time}}& \mathrm{ℰ}\left(t\right)\times {10}^{-4}& \\ 0& 1.442426615& \\ 15000& 1.442448\dots & \end{array}\mathit{\delta }\mathbf{=}\mathbf{0.01}:\begin{array}{ccc}\mathit{\text{time}}& \mathrm{ℰ}\left(t\right)\times {10}^{-4}& \\ 0& 1.442426615& \\ 150000& 1.442426620& \end{array}$$E31

For $\delta ={10}^{-3}$, there is variation in $\mathrm{ℰ}$ in the 11th decimal.

4.4 Multiple reflections/transmissions within dielectric cylinder

We now examine the scattered electromagnetic fields—particularly the polarization $E_z\left(x,y\right)$_—within and in the vicinity of the dielectric cylinder. These plots complement the global scattered $E_z\left(x,y\right)$ in Figures 3a,4a, and 5a, but for better resolution we choose $\delta =0.01$ and a slightly different ratio of pulse width to dielectric cylinder diameter. In Figures 8–13, the perspective is looking down from above with the 1D pulse propagating from left to right ($\to$), seen as a dark vertical band. The dielectric cylinder appears as a pink cylinder with the smaller darker pink being the base of the cylinder. The time is expressed in normalized time: $t=t^{\ast }/10$, where $t^{\ast }$ is the QLA time for $\delta =0.01$.

Figure 8.

A view from the $z$-axis of the 1D incident $E_z$ wavefront, with x-y the plane of the page. The vacuum pulse is propagating in the $x$-direction, $\to$ (a) the 1D incident pulse has encountered the localized dielectric cylinder, with both transmission and reflection at the thin vacuum-dielectric boundary layer. The reflected $E_z$ circular wavefront undergoes a $\pi$-phase change. (b) the transmitted $E_z$, within the dielectric, has a lower phase speed and so lags the 1D vacuum pulse.

Figure 9.

As the 1D vacuum part of the wavefront moves past the dielectric cylinder, the two vacuum-dielectric boundary “points” move closer together: (a) at t = 18 k, (b) at t = 22.2 K. the vacuum-reflected wavefront keeps radiating out.

Figure 10.

Wavefronts of $E_z$ at times (a) t = 28.2 k, and (b) t = 32 k around and within the dielectric cylinder after the original 1D pulse has moved past the dielectric. (a) the pinching of the two boundary “points” results in a focussing of $E_z$ and its subsequent spiking at $t=28.2k$. This spike now propagates toward the backend of the dielectric cylinder, (b), and “diffuses,” one should also note the wavefront emanating from the 1D vacuum pulse.

Figure 11.

Wavefronts of $E_z$ at times around and within the dielectric cylinder. At (a) t = 38.4 k the transmitted $E_z$ within the dielectric is radiating outward, with one part reaching the back of the dielectric and resulting in a complex transmission into the vacuum region at the back end of the dielectric, (b) at t = 43.2 k, and a complex reflection back into the dielectric. There is no phase change in the reflected $E_z$.

Figure 12.

Wavefronts of $E_z$ at times (a) t = 47.4 k, and (b) t = 52.2 k around and within the dielectric cylinder. The major vacuum wavefront that is transmitted out of the dielectric now radiates out in the xy-plane. The two boundary contact “points” of the wavefront are now propagating back to the front of the dielectric cylinder, as clearly seen in (a) and (b). These localized wavefronts will have their global wavefronts similar to those shown in Figures 3a,4a, and 5a.

Figure 13.

Wavefronts of $E_z$ at times (a) t = 55.8 k, and (b) t = 60 k around and within the dielectric cylinder.

In Figure 8a, at time $t=7.2k$, a part of the incident pulse has just entered the dielectric cylinder with the transmitted $E_z$ field starting to lag behind the main 1D vacuum pulse since $1<n_{\mathit{cyl}}$. Also, the reflected part of $E_z$ emanates from the two boundary points at the sharp vacuum-dielectric boundary and has undergone a $\pi$-phase change because the incident 1D pulse is propagating from low to high refractive index. In Figure 8b, the slower transmitted $E_z$ wavefront within the dielectric is very evident, as is the reflected part of $E_z$ back into the vacuum.

By $t=18k$, Figure 7b, the 1D pulse has propagated past the dielectric. The $E_z$-field within the dielectric is now being focussed due to its motion toward the backend of the cylinder, with its increasing amplitude but reduced base. As it reaches the backend of the dielectric, part of $E_z$ will be transmitted into the vacuum while the other part will be reflected back into the dielectric but now without any phase change since the pulse is propagating from high to low refractive index.

5. Dissipative classical systems, open quantum systems, and Kraus representation

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 $\mid \psi \left(t\right)⟩$ is given by the unitary evolution from the Schrodinger equation: $\mid \psi \left(t\right)⟩=U\left(t\right)\mid \psi \left(0\right)⟩$ with $U=\mathit{exp}\left[-{\mathit{itH}}_0\right]$ unitary for the Hermitian Hamiltonian $H_0$. The evolution of the density matrix, $\rho =\mid \psi ⟩⟨\psi \mid$, is governed by the corresponding von Neumann equation: $\rho \left(t\right)=U\left(t\right)\rho \left(0\right)U^{†}\left(t\right)$. 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

$$\rho \left(t\right)=\sum _k{K}_k\rho \left(0\right)K_k^{†},\text{with}\sum _k{K}_k^{†}{K}_k=I$$E32

where the only restriction on the set of so-called Kraus matrices $K_k$ is that the sum of $K_k^{†}{K}_k$ is the identity matrix. The evolution of the density matrix, Eq. (32), is no longer unitary for $k\ge 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 nontrivial example is the $1D$ Maxwell equations (without sources) for a homogeneous scalar medium with electrical losses,

$$i\frac{\partial }{\partial t}\left[\begin{array}{c}{E}_y\\ H_z\end{array}\right]=\left[\begin{array}{cc}0& \frac{{\epsilon }^{\ast }}{{\left|\epsilon \right|}^2}{\widehat{p}}_x\\ \frac{1}{{\mu }_0}{\widehat{p}}_x& 0\end{array}\right]\left[\begin{array}{c}{E}_y\\ H_z\end{array}\right]$$E33

with complex permittivity $\epsilon ={\epsilon }_R+i{\epsilon }_I.$${\epsilon }^{\ast }={\epsilon }_R-i{\epsilon }_I.$${\widehat{p}}_x=-i{\partial }_x$ is the momentum operator. Introducing the Dyson map $\mathit{\rho }=\mathit{diag}\left(|\epsilon |/\sqrt{{\epsilon }_R},\sqrt{{\mu }_0}\right)$ into Eq. (33) and after some algebraic manipulations the evolution equation can be written as

$$i\frac{\partial \mathbf{Q}}{\partial t}=\left[v_{\delta }\left({\sigma }_x+\frac{1}{2}{\mathit{\delta \sigma }}_y\right){\widehat{p}}_x-\frac{i}{2}\delta v_{\delta }{\sigma }_x{\widehat{p}}_x\right]\mathbf{Q},$$E34

where the state vector $\mathbf{Q}=\mathit{\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\left(1+{\delta }^2\right)}$ and ${\sigma }_x,{\sigma }_y,{\sigma }_z$ are the Pauli matrices.

5.1 Classical dissipation as a quantum amplitude damping channel

In symbolic form, the Maxwell equations with electric resistive losses, Eq. (34), can be written in the Schrodinger-form

$$i\frac{\partial \mid {\psi }_S⟩}{\partial t}=\left({\widehat{H}}_0\left(\mathbf{r}\right)-i{\widehat{H}}_1\left(\mathbf{r}\right)\right)\mid {\psi }_S⟩$$E35

where the Hamiltonians ${\widehat{H}}_0$ and ${\widehat{H}}_1$ are Hermitian, and the dissipative operator $i{\widehat{H}}_1$ is anti-Hermitian and positive definite. The positive definiteness requirement for the specific case of propagation in a lossy medium translates to

$$\mathit{Im}\left[E_y^{\ast }\frac{\partial H_z}{\partial x}+H_z^{\ast }\frac{\partial E_y}{\partial x}\right]>0,\text{with}{\epsilon }_I>0.$$E36

In general, the dissipative operator ${\widehat{H}}_1$ is relatively simple and models the phenomenological or coarse-graining of the underlying microscopic dissipative processes.

We aim to represent the dissipation in the Schrodinger picture, Eq. (35), as an open quantum system S interacting with its environment $\mathit{Env}$. The full system, $S+\mathit{Env}$, is closed, and hence its time evolution is unitary. Let $\widehat{\mathcal{U}}$ be this unitary operator, and $\rho$ the total density matrix with $\widehat{\mathcal{U}}:\rho \left(0\right)\to \rho \left(t\right)$. We make the usual assumption that the initial total density matrix is separable into the system and into the environmental Hilbert spaces: $\rho \left(0\right)={\rho }_S\left(0\right)\otimes {\rho }_E\left(0\right)$. A quantum operation $E$ on the open system of interest is defined as the map that propagates the open system density in time $t$:

$${\rho }_S\left(t\right)=E\left({\rho }_S\left(0\right)\right).$$E37

But under the conditions of initial separability, the action of the full unitary operators on the total density matrix will yield, after taking the trace over the environment,

$${\rho }_S\left(t\right)={\mathit{Tr}}_E\left(\rho \left(t\right)\right)={\mathit{Tr}}_E\left(\widehat{\mathcal{U}}\rho \left(0\right){\widehat{\mathcal{U}}}^{†}\right)$$E38

Assuming a stationary environment, ${\rho }_E\left(0\right)=\mid a⟩⟨a\mid$, Eq. (38) can be written as

$${\rho }_S\left(t\right)=\sum _{\mu }{\widehat{K}}_{\mu }{\rho }_S\left(0\right){\widehat{K}}_{\mu }^{†}.$$E39

where ${\widehat{K}}_{\mu }=⟨\mu |\widehat{\mathcal{U}}|a⟩$. These operators ${\widehat{K}}_{\mu }$ will form a Kraus representation for the quantum operation E for an open system, Eq. (37), provided the so-called Kraus operators satisfy the extended “unitarity” condition

$$\sum _{\mu }{K}_{\mu }^{†}{K}_{\mu }=I$$E40

Note that the individual Kraus operator need not be unitary. Based on this framework for open quantum systems, we proceed to construct a physical unitary dilation for the combined system-environment by identifying dissipation as an amplitude-damping operation, [21].

Let $d$ be the dimension of the system Hilbert space, and $r$ the dimension of the dissipative Hamiltonian $H_1$, Eq. (35). We require $d\ge 2r$, for optimal results but the dilation technique can be also applied to systems with $d=r$. If the system was quantum mechanical in nature, then there can be a set of $d^2$ Kraus operators at most. The matrix representation of the total unitary dilation evolution operator consists of listing all the Kraus matrices in the first column block. The remaining columns must then be determined, so that $\widehat{\mathcal{U}}$ is unitary. This unitary dilation is equivalent to the Stinespring dilation theorem [25]. The advantage of the Kraus approach is that it avoids the need to actually know the physical properties of the environment.

Returning to the Schrodinger representation of the classical system Eq. (35), one can employ the Trotter-Suzuki expansion to $\mathit{exp}\left[-\mathit{i\delta t}\left({\widehat{H}}_0-i{\widehat{H}}_1\right)\right]$

$$\mid \psi \left(\mathit{\delta t}\right)⟩=\left[e^{-\mathit{i\delta t}{\widehat{H}}_0}\cdot e^{-\mathit{\delta t}{\widehat{H}}_1}+O\left(\delta t^2\right)\right]\mid \psi \left(0\right)⟩.$$E41

Even though $\mathit{exp}\left(-\mathit{\delta t}{\widehat{H}}_1\right)$ is not unitary, ${\widehat{H}}_1$ is Hermitian and can be diagonalized by a unitary transformation $U_1$

$${\widehat{H}}_1={\widehat{U}}_1{\widehat{D}}_1{\widehat{U}}_1^{†}\text{with diagonal}{\widehat{D}}_1=\mathit{diag}\left[{\gamma }_1,\dots ,{\gamma }_r\right],$$E42

where ${\gamma }_i>0$ are the dissipative rate eigenvalues of ${\widehat{H}}_1$. Thus, Eq. (41) becomes

$$\mid \psi \left(\mathit{\delta t}\right)⟩=\left[e^{-\mathit{i\delta t}{\widehat{H}}_0}{\widehat{U}}_1{\widehat{K}}_0{\widehat{U}}_1^{†}+O\left(\delta t^2\right)\right]\mid \psi \left(0\right)⟩,$$E43

where ${\widehat{K}}_0$ is

$${\widehat{K}}_0=\left[\begin{array}{cc}{\widehat{\mathtt{\Gamma }}}_{r\times r}& 0_{r\times r}\\ 0_{\left(d-r\right)\times \left(d-r\right)}& I_{\left(d-r\right)\times \left(d-r\right)}\end{array}\right],\text{with diagonal}{\widehat{\mathtt{\Gamma }}}_{r\times r}=\mathit{diag}\left[e^{-{\gamma }_1\mathit{\delta t}}\dots e^{-{\gamma }_r\mathit{\delta t}}\right].$$E44

The non-unitary ${\widehat{K}}_0$ will be one of our Kraus operators, and it describes the physical dissipation in the open system. We must now introduce a second Kraus operator ${\widehat{K}}_1$, so that ${\widehat{K}}_0^{†}{\widehat{K}}_0+{\widehat{K}}_1^{†}{\widehat{K}}_1=\mathrm{ℐ}$:

$${\widehat{K}}_1=\left[\begin{array}{cc}{0}_{\left(d-r\right)\times \left(d-r\right)}& 0_{\left(d-r\right)\times \left(d-r\right)}\\ \sqrt{I_{r\times r}-{\widehat{\mathtt{\Gamma }}}^2}& 0_{r\times r}\end{array}\right].$$E45

${\widehat{K}}_1$ represents a transition that is not of direct interest. These Kraus operators ${\widehat{K}}_0,{\widehat{K}}_1$ are the multidimensional analogs of the quantum amplitude damping channel [21]: with ${\widehat{K}}_0$ corresponding to the dissipation processes, while ${\widehat{K}}_1$ corresponds to an unwanted quantum transition.

The block structure of the final unitary dilation evolution operator ${\widehat{\mathcal{U}}}_{\mathit{\text{diss}}}$, corresponding to the non-unitary dissipation operator $e^{-\mathit{\delta t}{\widehat{H}}_1}$, consists of column blocks of the Kraus operators ${\left(K_0{K}_1\dots \right)}^T$, and the remaining column blocks are of those matrices required to make ${\widehat{\mathcal{U}}}_{\mathit{\text{diss}}}$ unitary [21]:

$${\widehat{\mathcal{U}}}_{\mathit{\text{diss}}}=\left[\begin{array}{cccc}\widehat{\mathtt{\Gamma }}& 0& 0& -\sqrt{I_{r\times r}-{\widehat{\mathtt{\Gamma }}}^2}\\ 0& I_{\left(d-r\right)\times \left(d-r\right)}& 0& 0\\ 0& 0& I_{\left(d-r\right)\times \left(d-r\right)}& 0\\ \sqrt{I_{r\times r}-{\widehat{\mathtt{\Gamma }}}^2}& 0& 0& \widehat{\mathtt{\Gamma }}\end{array}\right].$$E46

Thus, it can be shown that the evolution of the system $\mid {\psi }_0⟩$ and environment $\mid 0⟩$ is given by

$$\mid 0⟩\mid {\psi }_0⟩=\frac{1}{\sqrt{E_0}}\sum _q{\psi }_{0q}\mid 0⟩\mid q⟩\to \mid 0⟩\otimes e^{-\mathit{i\delta t}{\widehat{H}}_0}{\widehat{U}}_1{\widehat{K}}_0{\widehat{U}}_1^{†}\mid {\psi }_0⟩+\mid 1⟩\otimes {\widehat{K}}_1\mid {\psi }_0⟩,$$E47

where measurement of the first qubit-environment by $\mid 0⟩⟨0\mid \otimes I_{d\times d}$ yields a state analogous to $\mid 0⟩\mid \psi \left(\mathit{\delta t}\right)⟩$. Finally, on taking the trace over the environment will yield the desired system state $\mid \psi \left(\mathit{\delta t}\right)⟩$. The corresponding quantum circuit for Eq. (47) is shown in Figure 14.

Figure 14.

Quantum circuit diagram for Eq. (47).

It is important to highlight that the implementation of the dissipative case is directly overlapping with the QLA framework. The QLA can be used to implement the $exp\left[-\mathit{i\delta t}{\widehat{H}}_0\right]$ part in Eq. (43) as proposed in the previous sections. Specifically, for the lossy medium, the exponential operator of the ${\widehat{H}}_0$ term in Eq. (35) can be easily handled with QLA [12, 13].

6. Conclusions

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 ${\widehat{P}}_0=\mid 0⟩⟨0\mid \otimes I_{d\times d}$ is

$$p\left(0\right)=\sum _{i=1}^r{e}^{-2{\gamma }_i\mathit{\delta t}}{\left|{\psi }_{i0}\right|}^2+\sum _{i=r+1}^d{\left|{\psi }_{i0}\right|}^2\ge 1+\left(e^{-2{\gamma }_{\mathit{max}}\mathit{\delta t}}-1\right)\sum _{i=1}^r{e}^{-2{\gamma }_i\mathit{\delta t}}{\left|{\psi }_{i0}\right|}^2$$E48

The form of the unitary operator ${\widehat{\mathcal{U}}}_{\mathit{\text{diss}}}$, Eq. (42) implies that it can be decomposed into r two-level unitary rotations ${\widehat{R}}_y\left({\theta }_i\right)$ with $cos\left({\theta }_i/2\right)=e^{-{\gamma }_i\mathit{\delta t}}.$ Then, the quantum circuit implementation of ${\widehat{\mathcal{U}}}_{\mathit{\text{diss}}}$ requires $O\left({\mathit{\text{rlog}}}_2^{2}d\right)$ CNOT - and ${\widehat{R}}_y\left({\theta }_i\right)$ quantum gates, so that there is an improvement in the circuit depth of $O\left(r/d^2\right)$. 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.

Acknowledgments

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.