# Quantum Computing Perspective for Electromagnetic Wave Propagation in Cold
Magnetized Plasmas

Efstratios Koukoutsis<EMAIL_ADDRESS>Kyriakos Hizanidis  School of
Electrical and Computer Engineering, National Technical University of Athens,
Zographou 15780, Greece  George Vahala  Department of Physics, William & Mary,
Williamsburg, Virginia 23187, USA  Min Soe  Department of Mathematics and
Physical Sciences, Rogers State University, Claremore, Oklahoma 74017, USA
Linda Vahala  Department of Electrical and Computer Engineering, Old Dominion
University, Norfolk, Virginia 23529, USA  Abhay K. Ram  Plasma Science and
Fusion Center, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA

###### Abstract

The study of electromagnetic wave propagation in magnetized plasmas is of
paramount importance in various fields, including astrophysics, fusion energy,
and communication systems. In thermonuclear fusion experiments where transient
interaction phenomena between electromagnetic waves and plasma can disrupt the
overall confinement, we have to rely on the modern state of the art,
computational tools to delve into the physics of wave propagation in plasma.
However, even those sophisticated computational methods are facing challenges
in terms of memory resources and speed when they are forced to capture all the
physical processes that occur in wave-plasma interaction. Simultaneously, the
rapidly advancing field of quantum technologies has opened up exciting new
frontiers in the computational studies, by promising a minimization on the
computational strain. In this paper we examine a theoretical quantum computing
re-conceptualization of Maxwell equations inside a cold, inhomogeneous,
magnetized plasma that can lead to quantum simulation of electromagnetic wave
propagation and scattering from inhomogeneities. By constructing a quantum
Schrodinger representation of Maxwell equations in plasma that admit unitary
-energy preserving- evolution we formulate a unitary product sequence of
operators that can form the basis of either a Qubit Lattice Algorithm (QLA) or
a pure quantum computing implementation. As an illustration of the power of
QLA, a full-wave simulation of wave-packet scattering from different shaped,
non-dispersive dielectrics is presented. QLAs when they are fully unitary,
they can be directly encoded into a quantum computer, further establishing
their versatility and capabilities but more importantly, indicating the impact
that quantum computers will have in the computational studies of wave
propagation in a fusion plasma.

##  I Introduction

Propagation of electromagnetic waves in thermonuclear fusion plasmas is one of
the most significant fields of research in the pursuit for magnetic fusion. In
magnetic confinement experiments, electromagnetic waves play a vital role in
plasma temperature control, localized non-inductive current drive, heating,
and plasma instability control. Therefore, there is an utmost need for
understanding the physics and mechanics of wave propagation and scattering
inside an inhomogeneous magnetized plasma to enable the optimization for
fusion applications.

While the bedrock for the theoretical and analytical studies of wave
propagation in plasmas has long been established,[1, 2] penetrating into the
complex processes that occur in plasmas and unraveling their physics require a
computational treatment. To that end, taking into consideration the
aforementioned importance of electromagnetic wave propagation in plasmas, a
plethora of computational tools have been developed,[3, 4, 5] ranging from
ray-tracing methods to full-wave simulations along with different domains of
application. Those state-of-the-art algorithmic tools can incorporate mode
conversion, Landau damping, cyclotron resonance damping through the complete
hot plasma dispersion relation as well as linear and quasi-linear loss
mechanisms and collisions.

However, solving the mathematical and physical problem of wave propagation in
an actual fusion device poses a challenge even for the most advanced
supercomputers. This is because the set of partial differential equations
describing a hot plasma is non-linear with a spatial-temporal dependence
inside a complex three-dimensional geometry. Therefore, computational
resources needed in terms of memory and speed are very difficult to be met
with in conventional computational systems.

In an attempt to resolve computational resources issues, new technologies are
emerging such as neuromorphic[6] or reservoir[7] computing, whereas reduced
models are adopted by retaining the physical features that are deemed
essential according to the plasma characteristics. The former have the
disadvantage that are oriented in solving problems associated with machine
learning and AI applications and are not yet directly applicable for
electromagnetic simulations. The latter display limited capability to showcase
relevant physical mechanisms that are present in wave propagation inside a
fusion plasma. For example, in short gradient scale lengths the Geometric
Optics[3] (GO) approach is expected to break down.

With classical computers eventually reaching their limits and fusion research
heavily relying on computational results we motivate a shift in the
traditional computational methods, engaging the modern and uprising quantum
technologies and quantum computing in particular.

Quantum computing is one of those computational pathways that can yield faster
computations than those achieved on a classical computer, [8, 9] the so called
quantum advantage, and has gained significant attention in the plasma physics
community. Considerations on general applications in plasma simulation can be
found in Ref. [10], whereas a fusion oriented review of possible quantum
computing applications is Ref.[11]. In Refs. [12] and [13] the authors exploit
the Quantum Signal Processing (QSP) protocol [14] for simulation of
electrostatic Landau damping and wave propagation in a cold fluid plasma
respectively. In addition, a quantum computing treatment for Vlasov equation
with collisions has been presented in Ref. [15]. Finally, a comprehensive
review on quantum computing applications in plasmas can be found in Ref.[16].

In this paper, we try to address the question whether wave propagation in the
simple fluid plasma model is amendable to quantum computing without tackling
the question of computational advantage over the classical methods. This is
accomplished by establishing the theoretical elements for a quantum computing
implementation of Maxwell equations in a cold, inhomogeneous, magnetized
plasma. Quantum computers are restricted to unitary operations following the
physical laws of closed quantum systems. Thus, the first step towards a
quantum implementation is to reformulate Maxwell equations as a quantum
Schrodinger equation with Hermitian structure. Then, the second challenge is
to decompose the relevant unitary operator of evolution into a product of
simpler 1-qubit and 2-qubit unitary operators that can be encoded efficiently
on a quantum computer. Special focus will be given on the full-wave simulation
results of electromagnetic wave propagation and scattering in a reduced case
of our formulation, for an inhomogeneous, tensorial dielectric without
dispersion, derived with Qubit Lattice Algorithms[17, 18, 19, 20] (QLAs).
Those simulations are implemented on classical supercomputers but can be also
directly transferred to quantum computers, acting as a precursor of QLA
generalization into cold magnetized plasma in the near term future.

This paper is organized as follows. Section II sets up the theoretical
formulation of Maxwell equations as a quantum Schrodinger equation in two
stages. In Sec.II.1 an augmented form of Maxwell equations in magnetized
plasma is presented, serving as a stepping stone for the construction of a
Schrodinger-Maxwell equation with unitary evolution in Sec.II.2. The
importance of initial and boundary conditions is discussed in Sec. II.3.
Decomposition of the established unitary evolution in a product formula of
simple unitary operators based on Trotterization is the main subject of
Sec.III.1. In sections III.2 and III.3 we present the algorithmic scheme of
QLA along with some initial value simulations for scattering of an
electromagnetic wave-packet from two-dimensional (2D) scalar, non-dispersive
inhomogeneous dielectric objects. In particular, we contrast the different
scattering characteristics from a local cylindrical dielectric with strong
gradients in the finite boundary layer between the dielectric and vacuum, with
that scattering from a local conic dielectric with weak boundary layer
gradients in the refractive index. Sec. III.4 presents the quantum
implementation of QLA and the respective scaling with the number of qubits.
Then, a commentary section III.5 follows, containing perspectives on both the
quantum and QLA implementation of the unitary product formula for the general
plasma case. Finally, in Sec.IV we discuss our results along with the next
necessary steps for an actual QLA implementation in the near future.

##  II Quantum representation of Maxwell equations in cold magnetized plasma

For a non-dispersive, tensorial and inhomgeneous medium, Maxwell equations can
be written as a Schrodinger equation with unitary evolution[21]

$i\partialderivative{\boldsymbol{\psi}}{t}=\hat{D}_{\rho}\boldsymbol{\psi},\quad\hat{D}_{\rho}=\hat{D}^{\dagger}_{\rho},\quad\boldsymbol{\psi}(\boldsymbol{r},0)=\boldsymbol{\psi}_{0},$
(1)

under a Dyson transformation $\hat{\rho}$ on the electromagnetic fields
$\boldsymbol{u}=(\boldsymbol{E},\boldsymbol{H})^{T}$, with
$\boldsymbol{\psi}=\hat{\rho}\boldsymbol{u}$. In particular, the Hermitian
operator $\hat{D}_{\rho}$

$\hat{D}_{\rho}=\hat{\rho}\hat{D}\hat{\rho}^{-1}=\hat{\rho}\hat{W}^{-1}(\boldsymbol{r})\hat{M}\hat{\rho}^{-1},$
(2)

with

$\hat{M}=i\begin{bmatrix}0_{3\times 3}&\boldsymbol{\nabla}\times\\\
-\boldsymbol{\nabla}\times&0_{3\times
3}\end{bmatrix},\quad\hat{W}=\begin{bmatrix}\epsilon(\boldsymbol{r})&0_{3\times
3}\\\ 0_{3\times 3}&\mu(\boldsymbol{r})\end{bmatrix}.$ (3)

In Eq.(3) the $\hat{M}$ operator is the Maxwell curl operator and the
Hermitian, positive definite $\hat{W}$ matrix represents the constitutive
relations of the medium. The explicit form of the Dyson map $\hat{\rho}$
depends on the structure of the material matrix $\hat{W}$:
$\hat{\rho}=\sqrt{\hat{W}}$.

On the other hand, the cold magnetized plasma as a dielectric medium is
characterized by dispersion. This translates into a frequency dependent
permittivity matrix $\tilde{\epsilon}(\omega)$ in the frequency domain.
Following the Stix notation[1],

$\tilde{\epsilon}(\omega)=\begin{bmatrix}S&-iD&0\\\ iD&S&0\\\
0&0&P\end{bmatrix}$ (4)

with

$\displaystyle S=$
$\displaystyle\epsilon_{0}\Big{(}1-\sum_{j=i,e}\frac{\omega^{2}_{pj}}{\omega^{2}-\omega_{cj}^{2}}\Big{)}$
$\displaystyle D=$
$\displaystyle\epsilon_{0}\sum_{j=i,e}\frac{\omega_{cj}\omega^{2}_{pj}}{\omega(\omega^{2}-\omega_{cj}^{2})}$
(5) $\displaystyle P=$
$\displaystyle\epsilon_{0}\Big{(}1-\sum_{j=i,e}\frac{\omega^{2}_{pj}}{\omega^{2}}\Big{)}.$

The definition of elements (II) in the Stix permittivity tensor is taken for a
two-species, ions (i) and electrons (e), plasma with inhomogeneous plasma
frequency
$\omega^{2}_{pj}(\bf{r})=\frac{n_{j}(\bf{r})q^{2}_{j}}{m_{j}\epsilon_{0}}$ and
cyclotron frequency $\omega_{cj}=\frac{q_{j}B_{0}}{m_{j}}$. The homogeneous
magnetic field $B_{0}$ is along the $z$ axis and $m_{j}$, $q_{j}$ are the mass
and charge of the $j$-species respectively. $n_{j}(\bf{r})$ is the $j^{th}$
species density.

###  II.1 Maxwell equations in temporal domain

In contrast to the optical response case, the temporal domain transformation
of $\tilde{\epsilon}(\omega)$ is expressed through a convolution integral. As
a result, the temporal domain, constitutive relations for a cold magnetized
plasma are

$\boldsymbol{d}=\hat{W}_{0}\boldsymbol{u}+\frac{1}{2\pi}\int_{0}^{t}\int_{-\infty}^{\infty}(\tilde{\epsilon}(\omega)-\epsilon_{0}I_{3\times
3})e^{-i\omega(t-\tau)}\boldsymbol{E}(\boldsymbol{r},\tau)d\,\omega d\,\tau,$
(6)

with $\boldsymbol{d}=(\boldsymbol{D},\boldsymbol{B})^{T}$. The matrix
$\hat{W}_{0}$ represents the optical response, as in Eq.(3), but now only that
of the vacuum. Evaluation of the inner integral term in Eq. (6) requires the
Plemelj formula[1] to yield

$\boldsymbol{d}=\hat{W}_{0}\boldsymbol{u}+\int_{0}^{t}\hat{K}(t-\tau)\boldsymbol{E}(\boldsymbol{r},\tau)d\,\tau,$
(7)

with the inhomogeneous susceptibility kernel $\hat{K}(t)$

$\hat{K}(t)=\epsilon_{0}\sum_{j=i,e}\begin{bmatrix}\frac{\omega^{2}_{pj}}{\omega_{cj}}\sin{\omega_{cj}t}&\frac{\omega^{2}_{pj}}{\omega_{cj}}(\cos{\omega_{cj}t}-1)&0\\\
\frac{\omega^{2}_{pj}}{\omega_{cj}}(1-\cos{\omega_{cj}t})&\frac{\omega^{2}_{pj}}{\omega_{cj}}\sin{\omega_{cj}t}&0\\\
0&0&\omega^{2}_{pj}t\end{bmatrix}.$ (8)

From the expressions (7) and (8), Maxwell equations for a cold magnetized
plasma now take the form

$i\partialderivative{\boldsymbol{u}}{t}=W_{0}^{-1}\hat{M}\boldsymbol{u}-i\int_{0}^{t}\partialderivative{\hat{G}(t-\tau)}{t}\boldsymbol{u}(\boldsymbol{r},\tau)d\,\tau$
(9)

where

$\partialderivative{\hat{G}(t)}{t}=\begin{bmatrix}\frac{1}{\epsilon_{0}}\partialderivative{\hat{K}}{t}&0_{3\times
3}\\\ 0_{3\times 3}&0_{3\times
3}\end{bmatrix},\quad\frac{1}{\epsilon_{0}}\partialderivative{\hat{K}}{t}=\sum_{j=i,e}\omega^{2}_{pj}(\bf{r})\begin{bmatrix}\cos{\omega_{cj}t}&-\sin{\omega_{cj}t}&0\\\
\sin{\omega_{cj}t}&\cos{\omega_{cj}t}&0\\\ 0&0&1\end{bmatrix}.$ (10)

###  II.2 Schrodinger representation

Returning back to $\tilde{\epsilon}(\omega)$ in Eq. (4), its Hermitian
structure ensures that the conductivity current does not produce dissipation
inside the plasma, i.e the cold magnetized plasma is a lossless dispersive
dielectric. Hence, it is possible to construct a Schrodinger representation of
Maxwell equations (9) that admit unitary evolution corresponding to
electromagnetic energy conservation. Such mathematical representations of
Maxwell equations for lossless dispersive media are well studied in the
literature[22, 23].

Defining the total conductivity current density $\boldsymbol{J}_{c}$ as

$\boldsymbol{J}_{c}=\int_{0}^{t}\partialderivative{\hat{K}}{t}\boldsymbol{E}(\boldsymbol{r},\tau)d\,\tau=\boldsymbol{J}_{ce}+\boldsymbol{J}_{ci},$
(11)

we exploit the rotational symmetry of $\partialderivative{\hat{K}}{t}$ in
Eq.(10) to reformulate Maxwell equations (9) as

$\displaystyle i\partialderivative{\boldsymbol{E}}{t}$
$\displaystyle=\frac{i}{\epsilon_{0}}\boldsymbol{\nabla}\times\boldsymbol{H}-\frac{i}{\epsilon_{0}}\boldsymbol{J}_{c},$
$\displaystyle i\partialderivative{\boldsymbol{H}}{t}$
$\displaystyle=-\frac{i}{\mu_{0}}\boldsymbol{\nabla}\times\boldsymbol{E},$
(12) $\displaystyle i\partialderivative{\boldsymbol{J}_{cj}}{t}$
$\displaystyle=i\epsilon_{0}\omega^{2}_{pj}(\boldsymbol{r})\boldsymbol{E}+\omega_{cj}\hat{S}_{z}\boldsymbol{J}_{cj},\quad
j=i,e.$

The set of equations (II.2) represent the augmented Maxwell system which self-
consistently describes the behaviour of electromagnetic fields inside a cold
magnetoplasma. We point out that Eq.(II.2) is the basis for FDTD
simulations,[24] but for a stationary plasma. The Hermitian matrix
$\hat{S}_{z}$,

$\hat{S}_{z}=\begin{bmatrix}0&-i&0\\\ i&0&0\\\ 0&0&0\end{bmatrix}$ (13)

represents the projection of spin-1 onto the $z$-axis.

To obtain an explicit Schrodinger representation of Eq.(II.2) we apply a Dyson
transformation[21],

$\hat{\rho}=diag(\epsilon^{1/2}_{0}I_{3\times 3},\mu^{1/2}_{0}I_{3\times
3},\frac{1}{\epsilon_{0}^{1/2}\omega_{pi}}I_{3\times
3},\frac{1}{\epsilon_{0}^{1/2}\omega_{pe}}I_{3\times 3})$ (14)

resulting in

$i\partialderivative{t}\begin{bmatrix}\epsilon_{0}^{1/2}\boldsymbol{E}\\\
\mu_{0}^{1/2}\boldsymbol{H}\\\
\frac{1}{e_{0}^{1/2}\omega_{pi}}\boldsymbol{J}_{ci}\\\
\frac{1}{e_{0}^{1/2}\omega_{pe}}\boldsymbol{J}_{ce}\end{bmatrix}=\begin{bmatrix}0_{3\times
3}&ic\boldsymbol{\curl}&-i\omega_{pi}&-i\omega_{pe}\\\
-ic\boldsymbol{\curl}&0_{3\times 3}&0_{3\times 3}&0_{3\times 3}\\\
i\omega_{pi}&0_{3\times 3}&\omega_{ci}\hat{S}_{z}&0_{3\times 3}\\\
i\omega_{pe}&0_{3\times 3}&0_{3\times
3}&\omega_{ce}\hat{S}_{z}\end{bmatrix}\begin{bmatrix}\epsilon_{0}^{1/2}\boldsymbol{E}\\\
\mu_{0}^{1/2}\boldsymbol{H}\\\
\frac{1}{e_{0}^{1/2}\omega_{pi}}\boldsymbol{J}_{ci}\\\
\frac{1}{e_{0}^{1/2}\omega_{pe}}\boldsymbol{J}_{ce}\end{bmatrix}\Leftrightarrow
i\partialderivative{\boldsymbol{\psi}}{t}=\hat{D}\boldsymbol{\psi}.$ (15)

It should be noted that we have switched from using the Riemann-Silberstein-
Weber[25] field representation to the vacuum field representation, and the
plasma inhomogeneity is now thrust into the source terms
$\boldsymbol{J}_{ci},\boldsymbol{J}_{ce}$ through the species plasma
frequencies $\omega_{pj}(\bf{r})$. Additionally, Eq.(15) can be easily
extended to incorporate different ions species by adding the respective ion-
species current components in the stave vector $\boldsymbol{\psi}$. In
realistic fusion experiments there will be hydrogen, deuterium and tritium
ions, so their contribution must be included in Eq.(15) for a complete
description of the total inhomogeneity profiles.

Under suitable Dirichlet boundary conditions the operator $\hat{D}$ in the
Schrodinger-Maxwell Eq.(15) is Hermitian. As a result, the evolution operator
$\hat{\mathcal{U}}=e^{-it\hat{D}}$ is unitary and corresponds to the
conservation of an extended electromagnetic energy $E(t)$ through the inner
product,

$E(t)=\innerproduct{\boldsymbol{\psi}}{\boldsymbol{\psi}}=\int_{\Omega}\Big{(}\epsilon_{0}\absolutevalue{\boldsymbol{E}}^{2}+\frac{\absolutevalue{\boldsymbol{B}}^{2}}{\mu_{0}}\Big{)}d\,\boldsymbol{r}+\int_{\Omega}\Big{(}\frac{\absolutevalue{\boldsymbol{J}_{ci}}^{2}}{\epsilon_{0}\omega^{2}_{pi}(\bf{r})}+\frac{\absolutevalue{\boldsymbol{J}_{ce}}^{2}}{\epsilon_{0}\omega^{2}_{pe}(\bf{r})}\Big{)}d\,\boldsymbol{r}=E(0)=\int_{\Omega}\Big{(}\epsilon_{0}\absolutevalue{\boldsymbol{E}_{0}}^{2}+\frac{\absolutevalue{\boldsymbol{B}_{0}}^{2}}{\mu_{0}}\Big{)}d\,\boldsymbol{r},\quad\Omega\subset\mathbb{R}^{3}.$
(16)

The extended electromagnetic energy Eq.(16) consists of two terms. The first
term is the standard electromagnetic energy in a vacuum whereas the second
term reflects the energy associated with the cold plasma response.

A subtlety related with the extended electromagnetic energy (16) is the
smoothness of $E(t)$ because of the Laplace Transform in Eq.(6). As a result,
even for resonant frequencies $\omega=\omega_{cj}$ we obtain a bounded
dispersive electromagnetic energy $E_{disp}(t)\leq E(0)$. Thus, it is possible
to quantify the resonant energization for each plasma population without
considering resonant wave-particle interactions or pertubative approximations
for the RF field.

###  II.3 Initial and boundary conditions

In this section we will restate our problem comparing the imposed mathematical
conditions with the ones in a plasma fusion device.

The plasma as a dielectric is considered to be confined inside a volume
$\Omega\subset\mathbb{R}^{3}$ with a boundary surface $\partial\Omega$. By
selecting the boundary condition

$\boldsymbol{n}\times\boldsymbol{E}=0,\quad\text{on {\hbox{\partial\Omega}}},$
(17)

the ”Hamiltonian operator” $\hat{D}$ in the Maxwell-Schrodinger equation (15)
is Hermitian so the standard quantum-mechanical analogies are present. In
fusion devices, the plasma is confined by a vacuum vessel at which the Perfect
Electric Conductor (PEC) boundary condition (17) no longer holds due to
electromagnetic losses in the walls. Alteration of the PEC boundary condition
results in the non-Hermicity of the operator $\hat{D}$ and subsequently, a
break in the unitary evolution. In this case, the quantum simulation of the
dynamics becomes troublesome. A resort has been proposed in Ref.[26] where
instead of the quantum simulation of the Maxwell dynamics, the linear system
of equations is solved through quantum singular value decomposition as a
boundary value problem. This approach could run into some difficulties as one
moves to 2D and 3D plasma wave propagation. Alternatively, one could resort to
some dilation by embedding the subsystem into a higher dimensional Hilbert
space and thereby recover unitarity within this higher dimensional space.

A different kind of boundary condition arises when one considers the
inhomogeneous electron and ion density profiles,
$\omega_{pi,e}(\boldsymbol{r})$. QLA is an initial value algorithm, i.e, no
internal boundary conditions are imposed in the evolution of the scattered
fields even though there exists various density structures within the plasma
with both sharp and slowly varying gradients (see e.g, Fig. 1).

For completeness, one could eventually introduce into the set of equations
(II.2) the effect of an antenna by coupling the Faraday equation with a
monochromatic oscillator[13]
$\boldsymbol{Q}(\boldsymbol{r},t)=\boldsymbol{Q}_{a}(\boldsymbol{r}_{a})e^{-i\omega_{a}t}$
with frequency $\omega_{a}$. The subscript $a$ denotes the antenna-related
quantities. In that way, the Faraday equation in (15) is augmented by

$\displaystyle i\partialderivative{(\mu_{0}^{1/2}\boldsymbol{H})}{t}$
$\displaystyle=-ic\boldsymbol{\curl}(\epsilon_{0}^{1/2}\boldsymbol{E})+\beta_{\boldsymbol{r},\boldsymbol{r}_{a}}\boldsymbol{Q}$
(18) $\displaystyle i\partialderivative{\boldsymbol{Q}}{t}$
$\displaystyle=\beta_{\boldsymbol{r},\boldsymbol{r}_{a}}(\mu_{0}^{1/2}\boldsymbol{H})+\omega_{a}\boldsymbol{Q},$

where
$\beta_{\boldsymbol{r},\boldsymbol{r}_{a}}=\beta\delta_{\boldsymbol{r},\boldsymbol{r}_{a}}$,
$\delta_{\boldsymbol{r},\boldsymbol{r}_{a}}$ is the Kronecker symbol and
$\beta$ is the coupling strength between the antenna emitted wave and the
magnetic field.

Finally we turn our attention to the initial conditions. The initial state
vector of Eq. (15) is

$\boldsymbol{\psi}(\boldsymbol{r},0)=\boldsymbol{\psi}_{0}=\begin{bmatrix}\epsilon_{0}^{1/2}\boldsymbol{E}_{0}\\\
\mu_{0}^{1/2}\boldsymbol{H}_{0}\\\ 0\\\ 0\end{bmatrix}.$ (19)

Inclusion of the antenna coupling Eq. (18) adds to the initial state
$\boldsymbol{\psi}_{0}$ the term
$\boldsymbol{Q}(\boldsymbol{r},0)=\boldsymbol{Q}_{a}$. The selection of the
initial vacuum electromagnetic filed profiles is dictated by the satisfaction
of the divergence set of Maxwell equations.

$\boldsymbol{\divergence}\boldsymbol{D}_{0}=\boldsymbol{\divergence}\boldsymbol{E}_{0}=0,\quad\boldsymbol{\divergence}\boldsymbol{B}_{0}=0.$
(20)

In that way, the divergence Maxwell equations are guaranteed to be satisfied
for $t>0$ along with $\boldsymbol{\divergence}\boldsymbol{J}_{cj}=0$ from the
charge continuity equation in the continuum limit. From current discrete
simulation 2D QLA runs[18, 20], it appears that divergence cleaning is not
required as QLA divergence errors are spatially localized and do not
accumulate.

##  III Connection with quantum computing and Qubit Lattice Algorithms

Application of QLA or any other quantum protocol for simulation of
electromagnetic wave propagation in a cold inhomogeneous magnetized plasma
requires a decomposition of the $\hat{D}$ operator in Eq.(15) into simpler
matrices,

$\hat{D}=\hat{D}_{vac}+\sum_{j=i,e}\hat{D}_{\omega_{pj}}+\hat{D}_{\omega_{cj}},$
(21)

with

$\displaystyle\hat{D}_{vac}$ $\displaystyle=-\frac{c}{2}(I_{2\times
2}+\hat{\sigma}_{z})\otimes\hat{\sigma}_{y}\otimes\boldsymbol{\curl}$ (22)
$\displaystyle\hat{D}_{\omega_{pi}}$
$\displaystyle=\frac{1}{2}\hat{\sigma}_{y}\otimes(I_{2\times
2}+\hat{\sigma}_{z})\otimes\omega_{pi}$ (23)
$\displaystyle\hat{D}_{\omega_{pe}}$
$\displaystyle=\frac{1}{2}(\hat{\sigma}_{x}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{x})\otimes\omega_{pe}$
(24) $\displaystyle\hat{D}_{\omega_{ci}}$
$\displaystyle=\frac{1}{4}(I_{2\times 2}-\hat{\sigma}_{z})\otimes(I_{2\times
2}+\hat{\sigma}_{z})\otimes\omega_{ci}\hat{S}_{z}$ (25)
$\displaystyle{D}_{\omega_{ce}}$ $\displaystyle=\frac{1}{4}(I_{2\times
2}-\hat{\sigma}_{z})\otimes(I_{2\times
2}-\hat{\sigma}_{z})\otimes\omega_{ce}\hat{S}_{z}.$ (26)

For simplicity let us assume that all quantities are only $x$-dependent,
rendering our model 1D. The inclusion of $y$\- and $z$-dependence is
straightforward, following the usual alternate direction iteration (ADI)
Cartesian procedures with no extraneous couplings of the respective quantum
operators. Then, the curl operator in Eq.(22) reads

$\boldsymbol{\curl}=\hat{S}_{x}\hat{p}_{x},\quad\hat{S}_{x}=\begin{bmatrix}0&0&0\\\
0&0&-i\\\ 0&i&0\end{bmatrix},\quad\hat{p}_{x}=-i\partialderivative{x}.$ (27)

###  III.1 Trotter Product Evolution Approximation

Trotterizing the total unitary evolution $e^{-i\delta t\hat{D}}$ whose
components are given in Eqs.(21)-(26) we obtain

$\boldsymbol{\psi}(\boldsymbol{r},\delta t)=e^{-i\delta
t\hat{D}_{vac}}\prod_{j=i,e}e^{-i\delta t\hat{D}_{\omega_{pj}}}e^{-i\delta
t\hat{D}_{\omega_{cj}}}\boldsymbol{\psi}_{0}+\textit{O}(\delta t^{2}).$ (28)

Each of the exponential operators in Eq.(28) can be written as a product of
unitary operators based on the their tensor-fold Pauli structure.
Specifically, we have the following diagonalization relations for the
$\hat{\sigma}_{y},\hat{\sigma}_{x},\hat{S}_{x},\hat{S}_{z}$ matrices

$\displaystyle\hat{\sigma}_{x}=\hat{H}\hat{\sigma}_{z}\hat{H},\quad$
$\displaystyle\hat{\sigma}_{y}=\hat{H}_{y}\hat{\sigma}_{z}\hat{H}_{y},$ (29)
$\displaystyle\hat{S}_{x}=\hat{H}^{(x)}_{y}\hat{\sigma}^{(x)}_{z}\hat{H}^{(x)}_{y},\quad$
$\displaystyle\hat{S}_{z}=\hat{H}^{(z)}_{y}\hat{\sigma}^{(z)}_{z}\hat{H}^{(z)}_{y},$

where $\hat{H}$ is the unitary Hadamard gate, $\hat{H}_{y}$ is the unitary
variant of Hadamard gate that diagonalizes $\hat{\sigma}_{y}$ whereas the
unitary set of matrices $\hat{H}^{(x)}_{y},\hat{H}^{(z)}_{y}$ and Hermitian
$\hat{\sigma}^{(x)}_{z},\hat{\sigma}^{(z)}_{z}$ are the three-dimensional
extensions of $\hat{H}_{y}$ and $\hat{\sigma}_{z}$ for $x$ and $z$ axes
respectively:

$\displaystyle\hat{H}$ $\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix},\quad\hat{H}_{y}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-i\\\
i&-1\end{bmatrix},\quad\hat{H}^{(x)}_{y}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&0&0\\\ 0&1&-i\\\
0&i&-1\end{bmatrix},$ (30) $\displaystyle H^{(z)}_{y}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-i&0\\\ i&-1&0\\\
0&0&1\end{bmatrix},\quad\hat{\sigma}^{(x)}_{z}$
$\displaystyle=\begin{bmatrix}0&0&0\\\ 0&1&0\\\
0&0&-1\end{bmatrix},\quad\hat{\sigma}^{(x)}_{z}$
$\displaystyle=\begin{bmatrix}1&0&0\\\ 0&-1&0\\\ 0&0&0\end{bmatrix}.$

This enables us to express the unitary exponential of operators (22)-(26)
using the identities

$\displaystyle e^{-i\delta
t\hat{V_{1}}\hat{A}\hat{V_{1}}^{\dagger}\otimes\hat{V_{2}}\hat{B}\hat{V_{2}}}=(\hat{V}_{1}\otimes\hat{V}_{2})e^{-i\delta
t\hat{A}\otimes\hat{B}}(\hat{V}^{\dagger}_{1}\otimes\hat{V}^{\dagger}_{2}),$
(31) $\displaystyle e^{-i\delta tI_{2\times 2}\otimes\hat{A}}=I_{2\times
2}\otimes e^{-i\delta t\hat{A}},$ (32) $\displaystyle
e^{-i\frac{\theta}{2}\hat{\sigma}_{i}\otimes\hat{A}}=I_{2\times
2}\otimes\cos{(\hat{A}\theta/2)}-i\hat{\sigma}_{i}\sin{(\hat{A}\theta/2)}.$
(33)

Therefore, the exponential operator $e^{-i\delta t\hat{D}_{vac}}$ can be
written

$e^{-i\delta t\hat{D}_{vac}}=\hat{C}_{vac}\hat{S}\hat{C}_{vac}$ (34)

where the unitary collision operator $\hat{C}_{vac}$ has the form

$\hat{C}_{vac}=I_{2\times 2}\otimes\hat{H}_{y}\otimes\hat{H}^{(x)}_{y},$ (35)

and the streaming operator in $x$, with $\delta x=c\delta t$:

$\hat{S}=\exp{i(I_{2\times
2}+\hat{\sigma}_{z})\otimes\hat{\sigma}_{z}\otimes\hat{\sigma}^{(x)}_{z}\delta
x\hat{p}_{x}/2}.$ (36)

Similarly, we express the rest of the operators in the Trotterized evolution
Eq.(28) as follows

$e^{-i\delta
t\hat{D}_{\omega_{pi}}}=\hat{C}_{\omega_{pi}}(\hat{\mathcal{R}}^{(pi)}_{z}\otimes
I_{3\times 3})\hat{C}_{\omega_{pi}},$ (37)

where $\theta_{pi}=\omega_{pi}\delta t$, $\hat{C}_{\omega_{pi}}$ is the
collision operator

$\hat{C}_{\omega_{pi}}=\hat{H}_{y}\otimes I_{2\times 2}\otimes I_{3\times 3}$
(38)

and the $\hat{\mathcal{R}}^{(pi)}_{z}$ operator is defined through identity
(33) which in principle represents a functional $\hat{R}_{i}(\cdot)$
rotations,

$\hat{\mathcal{R}}^{(pi)}_{z}=[\hat{R}_{z}(\theta_{pi})\otimes I_{2\times
2}]\hat{R}_{z}(\hat{\sigma}_{z}\theta_{pi}).$ (39)

For $e^{-i\delta t\hat{D}_{\omega_{pe}}}$ we obtain

$e^{-i\delta
t\hat{D}_{\omega_{pe}}}=\hat{C}^{(1)}_{\omega_{pe}}(\hat{\mathcal{R}}_{z}^{(pe)}\otimes
I_{3\times
3})\hat{C}^{(1)}_{\omega_{pe}}\hat{C}^{(2)}_{\omega_{pe}}(\hat{\mathcal{R}}_{z}^{(pe)}\otimes
I_{3\times 3})\hat{C}^{(2)}_{\omega_{pe}}$ (40)

with

$\displaystyle\hat{C}^{(1)}_{\omega_{pe}}$
$\displaystyle=\hat{H}\otimes\hat{H}_{y}\otimes I_{3\times 3},$ (41)
$\displaystyle\hat{C}^{(2)}_{\omega_{pe}}$
$\displaystyle=\hat{H}_{y}\otimes\hat{H}\otimes I_{3\times 3},$ (42)
$\displaystyle\hat{\mathcal{R}}_{z}^{(pe)}$
$\displaystyle=\hat{R}_{z}(\hat{\sigma}_{z}\theta_{pe}).$ (43)

We now move to the terms containing the cyclotron angle $\theta_{cj}$,

$\displaystyle e^{-i\delta t\hat{D}_{\omega_{ci}}}$
$\displaystyle=\hat{C}_{\omega_{ci}}[I_{4\times
4}\otimes\hat{R}_{z}^{(z)}(\theta_{ci}/2)][I_{2\times
2}\otimes\hat{R}_{z}(\hat{\sigma}^{(z)}_{z}\theta_{ci}/2)]$ (44)
$\displaystyle\times\hat{\mathcal{R}}^{(1),(ci)\dagger}_{z}\hat{\mathcal{R}}^{(2),(ci)\dagger}_{z}\hat{C}_{\omega_{ci}},$

with

$\hat{C}_{\omega_{ci}}=I_{2\times 2}\otimes I_{2\times
2}\otimes\hat{H}^{(z)}_{y}$ (45)

and operators
$\hat{R}_{z}^{(z)}(\theta_{ci}/2),\hat{\mathcal{R}}^{(1),(ci)}_{z},\hat{\mathcal{R}}^{(2),(ci)}_{z}$
representing $z$-rotation based on the $3\times 3$ $\hat{\sigma}^{(z)}_{z}$
matrix and functional $z$-rotations respectively,

$\displaystyle\hat{R}_{z}^{(z)}(\theta_{ci}/2)$
$\displaystyle=e^{-i\frac{\theta_{ci}}{4}\hat{\sigma}_{z}^{(z)}},$ (46)
$\displaystyle\hat{\mathcal{R}}^{(1),(ci)\dagger}_{z}$
$\displaystyle=e^{i\frac{\theta_{ci}}{4}\hat{\sigma}_{z}\otimes I_{2\times
2}\otimes\hat{\sigma}_{z}^{(z)}},$ (47)
$\displaystyle\hat{\mathcal{R}}^{(2),(ci)\dagger}_{z}$
$\displaystyle=e^{i\frac{\theta_{ci}}{4}\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}^{(z)}}.$
(48)

Finally,

$\displaystyle e^{-i\delta t\hat{D}_{\omega_{ce}}}$
$\displaystyle=\hat{C}_{\omega_{ce}}[I_{4\times
4}\otimes\hat{R}_{z}^{(z)}(\theta_{ce}/2)][I_{2\times
2}\otimes\hat{R}^{\dagger}_{z}(\hat{\sigma}^{(z)}_{z}\theta_{ce}/2)]$ (49)
$\displaystyle\times\hat{\mathcal{R}}^{(1),(ce)\dagger}_{z}\hat{\mathcal{R}}^{(2),(ce)}_{z}\hat{C}_{\omega_{ci}}.$

It is important to note that after we have made the somewhat standard leading-
order Trotterized approximation to the total unitary evolution operator in
Eq.(15), the evaluations of all the operators in Eqs.(34)-(49) are exact and
no further approximations have been made.

Consequently, the fully unitary evolution sequence reads

$\displaystyle\boldsymbol{\psi}(\boldsymbol{r},\delta
t)=\hat{C}_{vac}\hat{S}\hat{C}_{vac}\hat{C}_{\omega_{pi}}(\hat{\mathcal{R}}^{(pi)}_{z}\otimes
I_{3\times
3})\hat{C}_{\omega_{pi}}\hat{C}^{(1)}_{\omega_{pe}}(\hat{\mathcal{R}}_{z}^{(pe)}\otimes
I_{3\times
3})\hat{C}^{(1)}_{\omega_{pe}}\hat{C}^{(2)}_{\omega_{pe}}(\hat{\mathcal{R}}_{z}^{(pe)}\otimes
I_{3\times 3})\hat{C}^{(2)}_{\omega_{pe}}\hat{C}_{\omega_{ci}}[I_{4\times
4}\otimes\hat{R}_{z}^{(z)}(\theta_{ci}/2)]$ (50)
$\displaystyle\times[I_{2\times
2}\otimes\hat{R}_{z}(\hat{\sigma}^{(z)}_{z}\theta_{ci}/2)]\hat{\mathcal{R}}^{(1),(ci)\dagger}_{z}\hat{\mathcal{R}}^{(2),(ci)\dagger}_{z}\hat{C}_{\omega_{ci}}\hat{C}_{\omega_{ce}}[I_{4\times
4}\otimes\hat{R}_{z}^{(z)}(\theta_{ce}/2)][I_{2\times
2}\otimes\hat{R}^{\dagger}_{z}(\hat{\sigma}^{(z)}_{z}\theta_{ce}/2)]\hat{\mathcal{R}}^{(1),(ce)\dagger}_{z}\hat{\mathcal{R}}^{(2),(ce)}_{z}\hat{C}_{\omega_{ci}}\boldsymbol{\psi}_{0}.$

A QLA discretization of this unitary sequence (50) in the spatial domain
should lead to a desired quantum algorithm for simulation of wave propagation
in magnetized plasma with arbitrary density profile.

###  III.2 Example: QLA for scattering from 2D scalar non-dispersive
dielectric objects

To highlight the connection between the unitary product formula (50) with a
QLA sequence we briefly present the algorithmic scheme for a 2D $x-y$
scattering of a wave-packet from a scalar but non-dispersive localized
inhomogeneities with refractive index $n=n(x,y)$, as displayed in Fig.1.

(a)

(b)

Figure 1: Two different inhomogeneity refractive index profiles $1\leq
n(x,y)\leq 2$ and the electric field $E_{z0}(x)$ of the incident wave-packet.
The cylinder dielectric has strong spatial gradient near the vacuum-dielectric
interface, while the conic dielectric has very weak spatial gradients. In
Fig.1a these two profiles are shown superimposed. In Fig.1b the conic
dielectric is shown together with the incident wave-packet of arbitrary
normalization.

QLA’s were first developed in the late 1990’s to solve the 1D Schrodinger
equation using unitary collision and streaming operators acting on some qubit
basis[27, 28]. QLA recovered the Schrodinger equation in the continuum limit
to second order in the spatial lattice grid spacing $\delta$.

In our reduced case of non-dispersive dielectric, QLA is a discrete
representation of unitary representation of Maxwell equations (1), which, at a
mesoscopic level, uses an appropriately chosen interleaved sequence of three
non-commuting operators. Two of the operators are unitary collision and
streaming operators – the collision operator entangles the on-site qubits and
the streaming operator propagates the entangled state through the lattice. The
gradients in the medium constitutive properties are included via a third
operator referred to as a potential operator.

For 2D $x-y$ scattering of electromagnetic for a scalar dielectric state
vector that evolves unitarily is

$\boldsymbol{q}=\begin{bmatrix}nE_{x}\\\ nE_{y}\\\ nE_{z}\\\
\mu_{0}^{1/2}H_{x}\\\ \mu_{0}^{1/2}H_{y}\\\
\mu_{0}^{1/2}H_{z}\end{bmatrix}=\begin{bmatrix}q_{0}\\\ q_{1}\\\ q_{2}\\\
q_{3}\\\ q_{4}\\\ q_{5}\end{bmatrix}.$ (51)

In (diagonal) tensor dielectric media one would simply have $q_{0}\rightarrow
n_{x}E_{x}$, $q_{1}\rightarrow n_{y}E_{y}$, $q_{2}\rightarrow n_{z}E_{z}$.

The decomposition of the electromagnetic Schrodinger equation (1) in Cartesian
components is

$\displaystyle\partialderivative{q_{0}}{t}=\frac{1}{n}\partialderivative{q_{5}}{y},\quad\partialderivative{q_{1}}{t}=\frac{1}{n}\partialderivative{q_{5}}{x},\quad\partialderivative{q_{2}}{t}=\frac{1}{n}\Big{[}\partialderivative{q_{4}}{x}-\partialderivative{q_{3}}{y}\Big{]},$
(52)
$\displaystyle\partialderivative{q_{3}}{t}=\partialderivative{(q_{2}/n)}{y},\quad\partialderivative{q_{4}}{t}=\partialderivative{(q_{2}/n)}{x},$
$\displaystyle\partialderivative{q_{5}}{t}=-\partialderivative{(q_{1}/n)}{x}+\partialderivative{(q_{0}/n)}{n_{y}}.$

For the discrete QLA, using ADI, the unitary collision operators in the x and
y directions are

$\hat{C}_{X}=\begin{bmatrix}1&0&0&0&0&0\\\
0&\cos{\theta_{0}}&0&0&0&-\sin{\theta_{0}}\\\
0&0&\cos{\theta_{0}}&0&-\sin{\theta_{0}}&0\\\ 0&0&0&1&0&0\\\
0&0&\sin{\theta_{0}}&0&\cos{\theta_{0}}&0\\\
0&\sin{\theta_{0}}&0&0&0&\cos{\theta_{0}}\end{bmatrix},$ (53)
$\hat{C}_{Y}=\begin{bmatrix}\cos{\theta_{0}}&0&0&0&0&\sin{\theta_{0}}\\\
0&1&0&0&0&0\\\ 0&0&\cos{\theta_{0}}&\sin{\theta_{0}}&0&0\\\
0&0&-\sin{\theta_{0}}&\cos{\theta_{0}}&0&0\\\ 0&0&0&0&1&0\\\
-\sin{\theta_{0}}&0&0&0&0&\cos{\theta_{0}}\end{bmatrix}.$ (54)

with collision angle $\theta_{0}=\delta/4n$. The form of $\hat{C}_{X}$ can be
readily discerned from the coupling of the $\partialderivative{t}$ with
$\partialderivative{x}$ derivatives in (52): $q_{1}-q_{5}$, and $q_{2}-q_{4}$,
as well as the respective collision angle. Similarly for the unitary matrix
$\hat{C}_{Y}$.

We now define the unitary streaming operator $\hat{S}_{ij}$ which shifts the
amplitudes $\\{q_{i},q_{j}\\}$ one lattice unit, either in the $x$ or in the y
direction, while leaving all the other amplitudes unaffected. Then the
collide-stream sequence along each direction is,

$\displaystyle\hat{U}_{X}$
$\displaystyle=\hat{S}^{+x}_{25}\hat{C}^{\dagger}_{X}\hat{S}^{-x}_{25}\hat{C}_{X}\hat{S}^{-x}_{14}\hat{C}^{\dagger}_{X}\hat{S}^{+x}_{14}\hat{C}_{X}\,.\,\hat{S}^{-x}_{25}\hat{C}_{X}\hat{S}^{+x}_{25}\hat{C}^{\dagger}_{X}\hat{S}^{+x}_{14}\hat{C}_{X}\hat{S}^{-x}_{14}\hat{C}^{\dagger}_{X}$
(55) $\displaystyle\hat{U}_{Y}$
$\displaystyle=\hat{S}^{+y}_{25}\hat{C}^{\dagger}_{Y}\hat{S}^{-y}_{25}\hat{C}_{Y}\hat{S}^{-y}_{03}\hat{C}^{\dagger}_{Y}\hat{S}^{+y}_{03}\hat{C}_{Y}\
\,.\,\hat{S}^{-y}_{25}\hat{C}_{Y}\hat{S}^{+y}_{25}\hat{C}^{\dagger}_{Y}\hat{S}^{+y}_{03}\hat{C}_{Y}\hat{S}^{-y}_{03}\hat{C}^{\dagger}_{Y}.$

It should be noted that the first set of four collide-stream operators in
$\hat{U}_{X}$ and $\hat{U}_{Y}$ would yield (52) to first order in $\delta$.

The terms in (52), containing the derivatives of the refractive index, are
recovered through the following potential operators

$\hat{V}_{X}=\begin{bmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&-\sin{\beta_{0}}&0&\cos{\beta_{0}}&0\\\
0&\sin{\beta_{0}}&0&0&0&\cos{\beta_{0}}\end{bmatrix}$ (56)

and

$\hat{V}_{Y}=\begin{bmatrix}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\
0&0&\cos{\beta_{1}}&\sin{\beta_{1}}&0&0\\\
-\sin{\beta_{1}}&0&0&0&0&\cos{\beta_{1}}\end{bmatrix}.$ (57)

The angles $\theta_{0}$, $\beta_{0}$, and $\beta_{1}$, that appearing in
matrices (53), (54), (56), and (57) are chosen so that the discretized system
reproduces (52) to order $\textit{O}(\delta^{2})$.

The evolution of the state vector $\boldsymbol{q}$ from time $t$ to
$t+\Delta{t}$ is given by,

$\boldsymbol{q}(t+\Delta{t})=\hat{V}_{Y}\hat{V}_{X}\hat{U}_{Y}\hat{U}_{X}\boldsymbol{q}(t).$
(58)

Note that the external potential operators $\hat{V}_{X},\hat{V}_{Y}$, as given
above, are not unitary.

A detailed analysis of the QLA for the more general case of a bi-axial medium
along with simulation results for scattering of Gaussian pulses can be found
in Ref. [20].

###  III.3 QLA simulation results

The electromagnetic or electrostatic structures that propagate in plasmas are
generally not in the form of plane waves. Rather, they are wave-packets that
are localized in space such as spatially confined beams and of finite duration
in time. The interaction of the inhomogeneity plasma profile with the envelope
of the carrier wave, as well as with the individual components that a
spatially confined beam consists of, will lead to complex electromagnetic
structures that will affect the current densities in the dispersive plasma.
More importantly, those transport effects correspond to energy transfer from
the initial electromagnetic fields to the current density fields and can be
explicitly measured due to Eq.(16) which describes the extended
electromagnetic energy. Hence, examination of wave packet propagation in
plasmas is extremely important in realistic fusion experiments.

However, before tackling the propagation of such wave-packets in plasma which
is extremely complex it is instructive to investigate the behavior in a
simpler non-dispersive scalar medium to verify that our framework works
properly. The shape of the inhomogeneities, depicted in Fig.1a, can be related
to cylindrical filaments or smooth localized concentrations of plasma density.

In all simulations, the total energy is conserved to the seventh significant
digit.

(a)

(b)

Figure 2: QLA scattering simulation of $z$-component of an electromagnetic
pulse, $E_{z0}$ off a dielectric inhomogeneity in the shape of a cone
(Fig.2a), versus a cylindrical dielectric (Fig.2b). The perspective is looking
down the z-axis onto the x-y plane. The full-wave simulation for the wave-
cylinder encounter reveals strong initial reflection phenomena whereas the
reflection is very weak in the cone case. This differentiation in the wave
behavior is directly related to the steepness of the inhomogeneity gradient.
The weak reflected wave from the cone corresponds to asymptotic WKB type of
solution.

The initial electromagnetic wave-packet
$\boldsymbol{u}_{0}=(E_{z0}(x),-B_{y0}(x))^{T}$ is a Gaussian envelope with
internal oscillations, Fig.1b. The wave-packet propagates in the
$x$-direction, from a vacuum $n=1$ towards a localized dielectric
inhomogeneous object with $n_{max}(x,y)=2$. This polarization satisfies the
initial divergence conditions. As the 1D vacuum wave-packet interacts with the
2D refractive index of the dielectric. the $B_{y}$ field now becomes 2D, with
$B_{y}(x,y,t)$. This self-consistently generates a $B_{x}(x,y,t)$ so that
$\nabla\cdot\bf{B}=0$ as well as a 2D $E_{z}(x,y,t)$. Throughout the QLA
scattering simulation, $\nabla\cdot\bf{B}$ is monitored and is non-zero in
very small isolated spatial regions with some time instants in which
$max_{x,y}|\nabla\cdot\bf{B}/B_{0}|\leq 0.006$. $\nabla\cdot\bf{D}$ is
identically zero throughout the simulation. [For initiial
$E_{y0}(x)$-polarization, 2D QLA simulations retain $\nabla\cdot\bf{B}=0$
identically zero for all time.]

In Fig.2, the wave-packet has interacted with the dielectric object. The
viewpoint is looking down from the $z-$axis onto the $x-y$ plane. The apex of
the cone is seen as a white dot, while the interior of the dielectric cylinder
is in a somewhat darker color than the surrounding vacuum. In the case of a
dielectric cone, Fig.2a, there is a mild slowing down of that part of the
packet that is around the apex of the cone - since the phase velocity is
reduced to $c/n(x,y)$. But more importantly, one does not see any reflected
part of the packet from the slowly varying boundary region between vacuum and
dielectric. Basically the propagation is WKB-like. On the other hand there are
immediate reflection fronts emitted back into the vacuum from the interaction
of the wave-packet’s oscillation peaks with the steep refractive index
gradients in the boundary region of vacuum and cylinder dielectric, Fig.2b.
There is also considerable retardation in the oscillation peaks within the
dielectric cylinder as the refractive index away from the boundaries are
$n=2$.

As mentioned earlier, the transmitted component of the initial wave-packet
propagates into the respective dielectrics with phase velocity

$v_{ph}=\frac{c}{n(x,y)}$ (59)

because there is no dispersion in the media. However, the wave crests and the
envelope along the $y$-direction possess different phase velocities during
their propagation in the dielectric resulting to a lag between the interior
and outer wave components. Ultimately, both dielectrics act as a focusing lens
for the transmitted wave inside them. This latter behavior is clearly depicted
in Fig.3.

(a)

(b)

Figure 3: The propagation of the transmitted wave within the conical and
cylindrical dielectrics. The wave propagation is now distorted because the
initial wave crests along the $y-axis$ ”see” different effective length of
material. In both cases, Fig.3a, 3b, a focusing phenomenon is observed towards
the exit point of the transmitted wave to vacuum.

As the focused transmitted wave-front within the dielectric approaches the
vacuum boundary, the sudden change in the cylindrical dielectric object
produces a secondary internal reflection that propagates back inside the
cylinder. For the cone case, the smooth transition between the different
regions contributes a negligible secondary reflection. Those secondary
reflections, along with the secondary propagating wave-fronts in the vacuum
region are presented in Fig.4.

(a)

(b)

Figure 4: The absence of internal reflections from the conical dielectric
Fig.4a versus the internal reflections from the cylindrical dielectric Fig.4b.
Similar to the behavior of the primary reflections in Fig.2 the inhomogeneity
gradient of the dielectrics plays a pivotal role on the strength of the
internal reflection.

The back and forth succession from Fig.4 to Fig.2 through higher order
internal reflections in the cylindrical dielectric results in a radiating
temporal pattern. It should be noted that QLA is an initial value solver
giving the temporal (and transient) evolution of the scattered field without
the introduction of any internal boundary conditions to handle vacuum-
dielectric effects.

We will revisit the ability of QLA to generate the proper fields to fulfill
the divergence conditions in Sec.III.5, but for the plasma case.

###  III.4 Quantum encoding

To implement the QLA evolution (58) onto a quantum computer we must express
the participating matrices into elementary quantum gates acting on a set of
qubits. We will use two qubit registers. The first encodes the amplitude
dimensionality of the state vector $\boldsymbol{q}$ in (51), hence containing
$n_{i}=3$ qubits with $\\{\ket{i}\\}$ basis. The second register labels the
spatial $x-y$ discretization. For the two-dimensional lattice with $N$ nodes
and a discretization step $\delta$ in both directions, we will need
$n_{p}=\log_{2}N$ qubits with basis $\\{\ket{p}\\}$. Therefore, a total number
of $n_{total}=n_{p}+3$ qubits are required for the complete description of the
state $\boldsymbol{q}$.

The qubit encoding of the initial condition state vector $\boldsymbol{q}_{0}$
is,

$\ket{\boldsymbol{q}_{0}}=\sum_{p=0}^{M}\sum_{i=0}^{5}q_{0ip}\ket{i}\ket{p},$
(60)

with $M\leq N$ and amplitudes $q_{ip}$ characterize the $i$-component of the
state vector $\boldsymbol{q}$ in the lattice site $p$. normalized to the
square root of the initial (constant) electromagnetic energy so that
$\sum_{i}\absolutevalue{q_{ip}}^{2}=1$.

The collision operators $\hat{C}_{X}$ and $\hat{C}_{Y}$ can be each decomposed
into a product of two, two-level rotation matrices acting on the $\ket{i}$
register of the initial state (60). As a result, quantum implementation of the
collision operators requires $\textit{O}(9)$ CNOTs and single-qubit gates. On
the other hand, taking into consideration the quantum circuit implementation
of the streaming operators $\hat{S}^{+x}$ and $\hat{S}^{-x}$ in Ref. [21],
they can be decomposed into $\textit{O}(n^{3}_{px})$, $\textit{O}(n^{3}_{py})$
CNOTs and single-qubit gates respectively. Consequently, a quantum
implementation of the $\hat{U}_{Y}\hat{U}_{X}$ part in the QLA evolution (58)
is achieved, to leading order, within $\textit{O}(8n^{3}_{p})$ CNOTs and
single-qubit gates.

We refrain from a detailed description of the quantum implementation of the
non-unitary potential operators $\hat{V}_{X},\hat{V}_{Y}$ because it is not
relevant for the QLA sequence for a cold magnetized plasma which is fully
unitary. However, non-unitary operators can be handled using the LCU
method[29]. We direct the reader to Ref. [21] for a detailed discussion on the
quantum circuit implementation of these QLA non-unitary operators.

###  III.5 Discussion

The fully unitary product structure of Eq.(50) not only suggests that it can
be the building block of a QLA simulation but it is also directly encodable
onto a quantum computer. All the unitary collision $\hat{C}$’s operators are
in tensor product of elementary single-qubit gates like the Hadamard gate
$\hat{H}$, $\hat{H}_{y}=\hat{\sigma}_{z}\hat{R}_{x}(\pi/2)$ and the
$\hat{H}_{y}^{(z)},\hat{H}_{y}^{(x)}$ gates which can be easily implemented
within simple, two-qubit gates. As far as the unitary rotation operators are
concerned they are all diagonal and can be decomposed into simpler two-level
$z$-rotations or in worst case scenario, directly implemented within
$\textit{O}(2^{n})$ CNOTs and single-qubit gates [30], where $n=\log_{2}{N}$
is the number of qubits required for the quantum description of the state (see
(60)) in an $N$-node spatial discretization of the $x$-axis.

Comparing the Schrodinger representation of Maxwell equations for
inhomogeneous non-dispersive media Eq.(1) with Eq.(22) for the magnetized
plasma, it is seems that the latter supports more complexity due to the
dimensionality of the state vector $\boldsymbol{\psi}$. But, in contrast with
the optical case where the respective spatial displacement operator interferes
with the inhomogeneity of the refractive index (see Eq.(2)) the respective
exponential operator $e^{-i\hat{D}_{vacc}}$ in Eq.(22) is explicitly
decomposed without implicit dependence on the inhomogeneity plasma profile
which is reflected in the plasma frequencies. As a consequence, the expected
QLA will be free of the non-unitary potential operators such those introduced
in Eqs.(46),(47) resulting in a fully unitary product sequence similar to
Eq.(45).

Subsequently, the QLA sequence of $\hat{U}_{X}$ in Eq.(55) can be immediately
employed to calculate the term $e^{-i\delta t\hat{D}_{vac}}$ in the
Trotterized evolution approximation of $e^{-i\delta t\hat{D}}$,

$\displaystyle e^{-i\delta t\hat{D}}$ $\displaystyle=e^{-i\delta
t\hat{D}_{disp}}e^{-i\delta t\hat{D}_{vac}}+\textit{O}(\delta t^{2})$ (61)
$\displaystyle=e^{-i\delta t\hat{D}_{disp}}\hat{U}_{X}^{vac}+\textit{O}(\delta
t^{2}).$

Implementation of the dispersive part $e^{-i\delta t\hat{D}_{disp}}$, where
$\hat{D}_{disp}=\sum_{j=i,e}\hat{D}_{\omega_{pj}}+\hat{D}_{\omega_{cj}}$ can
be performed in parallel with the QLA. The main advantage of this
approximation is that we can decide whether to classically compute the
$\hat{U}_{X}^{vac}\boldsymbol{\psi}_{0}$, store the information and proceed
with a follow up quantum computation for the $e^{-i\delta t\hat{D}_{disp}}$
term resulting in a hybrid computation, or purely quantum computing the whole
sequence based on the quantum encoding of QLA as described in Sec.III.4.

In addition, both the QLA and its quantum implementation derived from unitary
evolution sequence (50) conserve the extended electromagnetic energy (16) and
the divergence conditions. Thus, no ”approximate” physics take place in our
full-wave scheme and the examined electromagnetic structures can be extended
beyond the usual plane-wave or monochromatic wave approximations as indicated
with the QLA simulations of wave-packet scattering from cylindrical and
conical dielectrics. The physical background of a QLA simulation is further
highlighted when applied for the plasma case. Assuming an initial X-wave
polarization $\boldsymbol{E}_{0}=E_{y}(k_{x}x)\hat{\boldsymbol{y}}$ the
scattering from a two dimensional $x-y$ plasma inhomogeneity will generate the
electromagnetic fields
$\boldsymbol{E}=E_{x}(k_{x}x,k_{y}y,\omega_{X}t)\hat{\boldsymbol{x}}+E_{y}(k_{x}x,k_{y}y,\omega_{X}t)\hat{\boldsymbol{y}}$
and $\boldsymbol{B}=B_{z}(k_{x}x,k_{y}y,\omega_{X}t)\hat{\boldsymbol{z}}$ but
most importantly will produce the conductivity current density
$\boldsymbol{J}_{cj}=J_{xcj}(k_{x}x,k_{y}y,\omega_{X}t)\hat{\boldsymbol{x}}+J_{ycj}(k_{x}x,k_{y}y,\omega_{X}t)\hat{\boldsymbol{y}}$
to satisfy
$\divergence\boldsymbol{E}=\divergence\boldsymbol{B}=\divergence\boldsymbol{J}_{cj}=0$.

Given the fact that the QLA scales linearly with the number of processors and
its quantum variant is probably expected to scale as $\textit{O}(n^{k}),\,k>2$
(see Sec.III.4), it is evident that our considerations pose a strong
alternative to the cost-inefficient FDTD methods, particularly in 2D and 3D.

On the other hand, it may be necessary to further manipulate the evolution
sequence (50) for an optimized QLA to be produced[31]. Therefore, considerable
research is required before applying the QLA for simulation of wave
propagation into a plasma characterized by fusion-reactor parameters. We also
reiterate that in applications of QLA to nonlinear spinor Bose-Einstein
condensates, the QLA produced an algorithm that was ideally parallelized to
all available cores on a classical supercomputer (over $750,000$ cores on the
now-retired IBM Blue Gene/$Mira$ supercomputer at Argonne).

##  IV Conclusions

The two main contribution of this paper are: (1) the analytical formulation of
Maxwell equations in a magnetized plasma, Eq.(15), as a Schrodinger equation,
and (2) a fully unitary QLA representation of this augmented Schrodinger
equation

The augmented Schrodinger representation has advantages over the standard
Helmholtz formulation[32, 33] both in the regularity of the spatial derivative
of the fields as well as in the construction of formal solutions. The
Hermitian structure of the full operator $\hat{D}$ permits a normal mode
decomposition of the solution in terms of the eigenfunctions
$\boldsymbol{\phi}(\boldsymbol{r},\lambda)$ of $\hat{D}$ operator with
$\lambda$ being the respective eigenvalues. This is very important in cases
where the inhomogeneous plasma profile does not possess a simple symmetry. In
addition, the unitary evolution of Eq.(15) explicitly preserves an extended
electromagnetic energy integral (16) beyond the usual Landau and Brillouin
approximations[34].

While various quantum schemes can be devised for the solution of the augmented
Schrodinger equation, we are currently pursuing the QLA scheme. For wave
propagation in a cold magnetized plasma, an appropriate QLA sequence of
unitary collision-streaming operators is determined in terms of 2-qubit gates.
The second part is based on the quantum representation of Maxwell equations in
which the energy preserving evolution is given by the unitary product formula
(50). This decomposition is deemed suitable for construction of a fully
unitary QLA, which no longer require the introduction of potential operators,
and their subsequent quantum encoding.

To benchmark the capabilities of QLA we present here the two dimensional
scattering of a wave-packet from either a cylindrical or a conical scalar,
inhomogeneous non-dispersive dielectrics. For the conic dielectric there are
weak spatial gradients in the layer connecting the vacuum to the dielectric.
As a result, there is negligible reflection at the first encounter of the wave
packet with the dielectric and then following the interaction with the steep
cone apex there is no internal reflections within the dielectric. This results
in a simple scattered field from the cone. However, for the cylindrical
dielectric, there is a sharp (but continuous) gradient in the layer connecting
the dielectric to the vacuum. The initial value QLA simulations (with no
internal boundary conditions being imposed at the dielectric-vacuum interface)
yield an immediate reflected wave front from the first interaction of the wave
packet with the dielectric followed by subsequent reflection/transmission of
the wave packet at the dielectric-vacuum layer. This leads to quite complex
interference in the scattered fields. Even though the final QLA is not fully
unitary due to the introduction of non-unitary potential operators, these
operators can be rewritten in the form of a linear combination of unitary
operators - and thus permit quantum computations, albeit requiring error
correcting qubits for the time evolution of the scattered field. Moreover, QLA
is ideally parallelized on classical supercomputers and seems to yield
alternate algorithms for the solution of classical problems.

We are now exploring QLA simulations of the wave propagation in a cold
magnetized (dispersive) plasma, exploiting the QLA operator splitting
approach. While only the x-dependent fully unitary QLA is presented here, the
use of the Alternating Direction Implicit (ADI) scheme will permit extensions
to fully 3D simulations.

###### Acknowledgements.

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 are however 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. This
research was partially supported by Department of Energy grants DE-SC0021647,
DE-FG0291ER-54109, DE-SC0021651, DE-SC0021857, and DE-SC0021653. 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.

## References

  * Stix [1992] T. H. Stix,  _Waves in plasmas_ (American Institute of Physics, 1992).
  * Swanson [2003] D. G. Swanson,  _Plasma Waves_ (Institute of Physics Publishing, 2003).
  * Friedland and Bernstein [1980] L. Friedland and I. B. Bernstein, “General geometric optics formalism in plasmas,” IEEE Trans. Plasma Sci. 8, 90–95 (1980).
  * Tsironis [2013] C. Tsironis, “On the Simplification of the Modeling of Electron-Cyclotron Wave Propagation in Thermonuclear Fusion Plasmas,” Prog. Electromagn. Res. B 47, 37–61 (2013).
  * Lau  _et al._ [2018] C. Lau, E. F. Jaeger, N. Bertelli, L. A. Berry, D. L. Green, M. Murakami, J. M. Park, R. I. Pinsker, and R. Prater, “AORSA full wave calculations of helicon waves in DIII-D and ITER,” Nucl. Fusion 58, 066004 (2018).
  * Hasler and Marr [2013] J. Hasler and H. Marr, “Finding a roadmap to achieve large neuromorphic hardware systems,” Front. Neurosci. 7 (2013), 10.3389/fnins.2013.00118.
  * G.Tanaka  _et al._ [2019] G.Tanaka, T. Yamane, J. B. Heroux, R.Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, _et al._ , “Strong Quantum Computational Advantage Using a Superconducting Quantum Processor,” Phys. Rev. Lett. 127, 180501 (2021).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, _et al._ , “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Dodin and Startsev [2021] I. Y. Dodin and E. A. Startsev, “On applications of quantum computing to plasma simulations,” Phys. Plasmas 28, 092101 (2021).
  * Joseph  _et al._ [2023] I. Joseph, Y. Shi, M. D. Porter, A. R. Castelli, V. I. Geyko, F. R. Graziani, S. B. Libby, and J. L. DuBois, “Quantum computing for fusion energy science applications,” Phys. Plasmas 30, 010501 (2023).
  * Engel, Smith, and Parker [2019] A. Engel, G. Smith, and S. E. Parker, “Quantum algorithm for the vlasov equation,” Phys. Rev. A 100, 062315 (2019).
  * Novikau, Startsev, and Dodin [2022] I. Novikau, E. A. Startsev, and I. Y. Dodin, “Quantum signal processing for simulating cold plasma waves,” Phys. Rev. A 105, 062444 (2022).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, “Optimal hamiltonian simulation by quantum signal processing,” Phys. Rev. Lett. 118, 010501 (2017).
  * Ameri  _et al._ [2023] A. Ameri, E. Ye, P. Cappellaro, H. Krovi, and N. F. Loureiro, “Quantum algorithm for the linear vlasov equation with collisions,” Phys. Rev. A 107, 062412 (2023).
  * Óscar Amaro and Cruz [2023] Óscar Amaro and D. Cruz, “A living review of quantum computing for plasma physics,” (2023), arXiv:2302.00001 [physics.plasm-ph] .
  * Vahala  _et al._ [2020] G. Vahala, L. Vahala, M. Soe, and A. K. Ram, “Unitary quantum lattice simulations for maxwell equations in vacuum and in dielectric media,” J. Plasma Phys. 86, 905860518 (2020).
  * Vahala  _et al._ [2022] G. Vahala, J. Hawthorne, L. Vahala, A. K. Ram, and M. Soe, “Quantum lattice representation for the curl equations of maxwell equations,” Radiat. Eff. Defects Solids 177, 85–94 (2022).
  * Ram  _et al._ [2021] A. K. Ram, G. Vahala, L. Vahala, and M. Soe, “Reflection and transmission of electromagnetic pulses at a planar dielectric interface: Theory and quantum lattice simulations,” AIP Advance 11, 105116 (2021).
  * Vahala  _et al._ [2023] G. Vahala, M. Soe, L. Vahala, A. K. Ram, E. Koukoutsis, and K. Hizanidis, “Qubit lattice algorithm simulations of maxwell’s equations for scattering from anisotropic dielectric objects,” Comput. Fluids 266, 106039 (2023).
  * Koukoutsis  _et al._ [2023] E. Koukoutsis, K. Hizanidis, A. K. Ram, and G. Vahala, “Dyson maps and unitary evolution for Maxwell equations in tensor dielectric media,” Phys. Rev. A 107, 042215 (2023).
  * Silveirinha [2015] M. G. Silveirinha, “Chern invariants for continuous media,” Phys. Rev. B 92, 125153 (2015).
  * Cassier, Joly, and Kachanovska [2017] M. Cassier, P. Joly, and M. Kachanovska, “Mathematical models for dispersive electromagnetic waves: An overview,” Comput. Math. with Appl. 74, 2792–2830 (2017).
  * Lee and Kalluri [1999] J. H. Lee and D. K. Kalluri, “Three-dimensional fdtd simulation of electromagnetic wave transformation in a dynamic inhomogeneous magnetized plasma,” IEEE Trans. Antennas Propag. 47, 1146–1151 (1999).
  * Khan [2005] S. A. Khan, “An Exact Matrix Representation of Maxwell’s Equations,” Phys. Scr. 71, 440 (2005).
  * Novikau, Dodin, and Startsev [2023] I. Novikau, I. Dodin, and E. Startsev, “Simulation of linear non-hermitian boundary-value problems with quantum singular-value transformation,” Phys. Rev. Appl. 19, 054012 (2023).
  * Boghosian and Taylor [1998] B. M. Boghosian and W. Taylor, “Simulating quantum mechanics on a quantum computer,” Physica D 120, 30–42 (1998).
  * Yepez [2002a] J. Yepez, “An efficient and accurate quantum lattice-gas model for the many-body schrodinger wave equation,” Comput. Phys. Commun. 146, 280–294 (2002a).
  * Childs, Kothari, and Somma [2017] A. M. Childs, R. Kothari, and R. D. Somma, “Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision,” SIAM J. Comput. 46, 1920–1950 (2017).
  * Bullock and Markov [2004] S. S. Bullock and I. L. Markov, “Asymptotically optimal circuits for arbitrary n-qubit diagonal comutations,” Quantum Inf. Comput. 4, 27––47 (2004).
  * Yepez [2002b] J. Yepez, “An efficient and accurate quantum algorithm for the dirac equation,” (2002b), arXiv:quant-ph/0210093 [quant-ph] .
  * Ram and Hizanidis [2016] A. K. Ram and K. Hizanidis, “Scattering of radio frequency waves by cylindrical density filaments in tokamak plasmas,” Phys. Plasmas 23, 022504 (2016).
  * Ram, Hizanidis, and Kominis [2013] A. K. Ram, K. Hizanidis, and Y. Kominis, “Scattering of radio frequency waves by blobs in tokamak plasmas,” Phys. Plasmas 20, 056110 (2013).
  * Jackson [1998] J. D. Jackson,  _Classical Electrodynamics_ (Wiley, 1998).