Quantum Algorithms in Electromagnetic Propagation Modelling for Telecommunications

Modelling of electromagnetic wave propagation in telecommunications has evolved from empirical models to highly deterministic ray tracing and numerical methods. The extraordinary computational effort of these methods and their wave nature seem ideally suited to be approached by quantum algorithms. We first examine current progress by reviewing recent proposals in the field. We scrutinize potentially advantageous quantum architectures, ranging from mainstream gate-level computers and near-term, mid-scale noisy architectures with limited capabilities, to adiabatic and annealing approaches that are already in commercial use. We analyze the weaknesses and strengths of recent proposals. Beyond the core algorithm, mechanisms to bridge the quantum and classical worlds are of particular interest. Extremely diverse algorithm specifications, from those based on Hamiltonian simulations and emulation of variational optimization to the unconstrained binary formulation, are compared with the use of pure gate-level circuits and known quantum subroutines. We show that the graph Laplacian, given its ability to integrate boundary conditions, is uniquely suited for quantum propagation modelling algorithms rooted in differential numerical methods. Quantum computers could overcome the temporal and spatial limitations of classical methods for larger computational domains and, to some extent, address the problems of dispersion and stability in finite-difference approximations. The ability to express the solution of a problem as an eigenvalue problem turns out to be an advantage in the quantum world, where eigenvalues and eigenvectors are inextricably intertwined with quantum mechanics. In this paper, we identify the most promising techniques and scenarios that hold the greatest potential.


I. INTRODUCTION
Without a solid understanding of electromagnetic (EM) wave propagation, the reliable, efficient, and secure operation of modern telecommunications networks would be impossible.EM propagation modelling enables network engineers to design, optimize, and scale wireless communications systems.This includes the positioning and setup of infrastructure such as antennas, routers and signal amplifiers.Modelling helps determine the best locations for these installations based on factors such as distance, terrain and building materials, all of which affect signal strength and quality.In addition, understanding the propagation of radio waves The associate editor coordinating the review of this manuscript and approving it for publication was Sandra Costanzo .
helps efficiently allocate and manage spectrum to avoid interference between different services and meet regulatory requirements.With the transition to higher frequency communications systems, such as 5G and beyond, understanding electromagnetic propagation becomes even more important, as these higher frequencies behave differently than the lower frequencies used in previous systems.
Empirical propagation models [1], while still commonly used in radio coverage planning, do not allow modelling of the spatial and temporal characteristics of the channel at the required level of detail.Parameters such as signal delay spread and direction of arrival are readily available in advanced deterministic models that incorporate detailed knowledge of the surrounding geometry.Two larger groups of deterministic algorithms are ray-tracing algorithm [2], [3] and full-wave numerical algorithms [4].Modelling EM propagation by ray tracing is based on the principles of geometrical optics and has become very popular in the last decade.However, the high-frequency approximation, in which the rays mimic narrow beams of light, exhibits significant deviations from the diffraction and scattering behavior at radio frequencies.The geometric theory of diffraction [5], [6] improves the accuracy of the modelling to some extent.On the other hand, if the problem is on the order of a few dozen wavelengths, EM propagation modelling based on numerical methods dominates [7].Numerical methods are based on Maxwell's partial differential equations of electrodynamics.Numerical solutions of these equations are a fundamental tool for the development of electronics, communication devices, lasers, antenna systems, and in many other areas with problems of similar magnitude.In the last decade, numerical methods have also been introduced for modelling telecommunications channels where the size of the problem increases to several hundred wavelengths or more [8].
Because of the extraordinary computational requirements, most numerical modelling in telecommunications is limited to two-dimensional geometries with many simplifications.It should be noted that there are other equally significant obstacles to large-scale modelling, such as numerical dispersion and instability, which are consequences of the finitedifference approximation.Even less sophisticated ray-tracing algorithms involve significant computational effort, which limits their applicability to smaller geometries or leads to reducing the dimensionality of the problem or modelling only the dominant propagation paths.
Quantum computing has brought remarkable advances in research and technology, including in the field of telecommunications with so-called quantum-assisted communications [9].Given this progress, the emergence of quantum algorithms for EM propagation modelling was anticipated.However, quite surprisingly, there have not been many proposals to date, even though this field would benefit greatly from the acceleration provided by quantum advantage.For clarity, 'quantum advantage' is a term used to describe the superiority of a quantum algorithm over the best-known classical algorithms, specifically in terms of its ability to solve a given problem more efficiently.
The potential power of quantum systems to simulate and solve complex problems intractable for classical computers was first recognized by Richard Feynman in 1982 [10].Feynman's ideas on quantum simulation laid the foundation for the field of quantum computing and quantum simulation and inspired subsequent research and advancement.The replacement of classical physics by quantum mechanics is considered by many to be the next revolution in computing, although the arguments of sceptics should not be ignored [11].The superposition property of quantum states, where a system of n qubits can be described by 2 n quantum amplitudes, allows massive parallelism on a scale not found in classical computing, while the entanglement of quantum states enables new forms of interactions and operations.The initial theoretical success in applying quantum computing to cryptanalysis was quickly replicated in other fields such as computational chemistry and optimization research.
Here we review current proposals of quantum algorithms in the field of EM propagation modelling.We perform an analysis of the strengths and weaknesses of each algorithm, provide time and space complexity where available, and discuss the preparation of initial states and extraction of results to bridge the quantum and classical worlds.We identify algorithmic features that are the source of the quantum advantage over classical algorithms.Two of the quantum algorithms discussed support ray tracing, while others deal with numerical simulations.
Since the algorithms cover a range of different quantum architectures, from those requiring fault-tolerant computation to variational quantum algorithms and quantum annealers, we give a brief introduction to each architectural paradigm; however, we do not cover the basics of quantum computing and quantum information theory here, as they can be found in numerous sources and textbooks [12].
Exploring quantum computing as an accelerator of simulation tasks is a compelling future research goal.We identify the most promising techniques and scenarios with the greatest potential in the field of EM propagation modelling, including the benefits of graph Laplacians for finite-differential approximation of discretized domains in Hilbert space and the importance of finding equivalent eigenvalue problems.Analyzing current work in this area, we aim to provide important insights, illuminate future prospects and challenges, and offer guidance to researchers and engineers in the field.We believe that overcoming the identified problems will lead to a dramatic improvement in the performance of propagation simulation, at least for some algorithms, while some of the concepts presented here will serve mainly as proofs of principle.
In the following, we first give an overview of propagation modelling in telecommunications in Section II.The quantum algorithms presented are designed to run on different quantum architectures, and the corresponding classes of algorithms are presented in Section III.Subsequently, quantum algorithms for EM propagation modelling are described in sections IV to X, including the definition of the actual problem to be solved, the basic quantum steps, and the space and time complexity, if available.Example problems are presented that were used to validate the algorithm.Section XI is a comparative analysis of the algorithms studied.Strengths and weaknesses are examined, including the practical potential of the algorithms and their quantum advantages over classical algorithms.Future outlook and prospective techniques are discussed.We conclude in Section XII.

II. PROPAGATION MODELLING IN TELECOMMUNICATIONS
Propagation modelling in telecommunications has a long history with numerous references in the literature.Chronologically, modelling techniques have evolved from stochastic 111546 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and empirical models to highly deterministic ray tracing and numerical approaches.We provide here an overview of the main modelling techniques, which is by no means exhaustive.

A. STOCHASTIC AND EMPIRICAL MODELS
We mention the group of stochastic and empirical models only for historical reasons, as they are unlikely to be accelerated by quantum computing.For a general overview of radio propagation models, see COST 273 [13].Stochastic models provide a simplified representation of the complex propagation phenomena that occur in wireless communication systems.Channel behavior is subject to random fluctuations, which are usually modeled by statistical distributions such as Gaussian or lognormal distributions.Empirical models target more diverse propagation environments [14], [15].They rely on extensive measurement campaigns and approximate radio channels by parametric functions and generally include an expression for path loss with an environment-specific exponent [1], [16], [17], [18], [19].Stochastic and empirical models offer a tradeoff between accuracy and computational efficiency.

B. RAY TRACING
In contrast to the empirical approach, ray tracing is a highly deterministic channel modelling.It allows the evaluation of advanced channel properties, such as delay spread or direction of arrival, but at the cost of higher processing overhead.The method is based on the principles of geometrical optics and involves tracing individual rays as they interact with objects and surfaces in the environment.The method assumes that EM waves propagate in a straight line until they encounter an object or interface at which reflection, diffraction, scattering, or transmission occur.This principle makes it possible to track individual rays and determine their paths.Reflection of radio waves occurs when they strike reflective surfaces such as walls, buildings or the ground.The laws of reflection are used to determine the angle of incidence and the angle of reflection, and the energy is redirected accordingly.Diffraction occurs when radio waves are refracted off the edges of objects or obstacles.When radio waves encounter interfaces between different media, such as air-ground or airbuilding interfaces, some of the energy is transmitted through the interface.Finally, scattering is the redirection of radio waves in different directions when they encounter rough or irregular surfaces.
Algorithms from this group refer to the base principle as ray launching [20], ray shooting and bouncing (SBR) [21], pincushion method [2], or more elaborated ray-tube [22] and beam tracing [23], the latter aggregating rays to reduce computational complexity and effectively converging to the second approach known as the method of images [3].In the method of images, a virtual image of the source is created on the opposite side of the reflecting surface.This image behaves as if it is emitting or receiving electromagnetic waves.Gaussian beam tracking [24] builds on the further improvements of image theory.
The common denominator of all ray-tracing algorithms is the high-frequency approximation of propagating waves, where a single ray mimics the behavior of a thin beam of light.The simplification most notably affects the accuracy of diffraction modelling, which for many practical purposes is negligible at optical frequencies [25].Geometrical Theory of Diffraction (GTD) introduces diffracted rays to approximate Maxwell's equations at the edge of two conducting half-planes [5] with discontinuity between the incident and reflected shadow regions.The Uniform Theory of Diffraction (UTD) proposed by Kouyoumjian and Pathak [6] softens the transition between the areas.The UTD was subsequently extended to handle diffraction edges with finite conductivity [26], [27].
Accurate ray tracing in buildings requires proper treatment of diffuse scattering [28], [29], [30], [31].Geometric optics, which is the basis of ray tracing, does not provide a satisfactory solution to this problem.Several approaches [29], [31] attempt to address this shortcoming, but all at the cost of significantly increased computational load and runtime [28].

C. NUMERICAL METHODS
The Finite-Difference Time-Domain (FDTD) method is arguably the simplest numerical technique.In the FDTD method, both the space and time domains are discretized and the field equations are approximated by finite-difference approximations.It was proposed by Kane Yee, who discretized Ampere's and Faraday's laws by second-order central differences in 1966 [4].The FDTD method uses a time-stepping approach in which the EM field equations are solved incrementally in discrete time steps.The Courant-Friedrich-Levy (CFL) constraint on the space-totime ratio is mainly responsible for the high computational cost.Therefore, early applications were limited to electrically small problems.To simulate open boundaries or interfaces between different media, suitable boundary conditions are needed.Commonly used boundary conditions in FDTD simulations include the perfectly matched layer (PML) absorbing boundary condition (ABC) or the Mur absorbing boundary condition.The ABCs usually require multiple layers of specialized nodes in space [32], [33], [34], [35], [36], [37].Neumann boundary condition effectively prevents the reflection of waves by setting the derivative of the field variable at the boundary to zero, mimicking an absorbing effect.Boundary conditions are an active area of research with numerous proposals in recent decades.
The sources of error in finite-difference full-wave methods are well understood and mathematically explained.The high computational cost and the progressive accumulation of errors with increasing propagation distance are the fundamental limitations.The FDTD method is subject to numerical dispersion, which is an error caused by the finite-difference approximations used to represent derivatives.Numerical dispersion can lead to spurious frequency-dependent effects and can alter the propagation characteristics of waves, especially in dispersive media or waveguides.Iterations lead to the accumulation of delay or phase errors, which show up as unphysical phenomena such as anisotropy, broadening and ringing of pulses, inaccurate wave cancellations, and virtual refractions.An exception is the one-dimensional version of the problem, where accurate computation is possible under proper conditions.Techniques such as refining the grid, smaller time steps, and proper implementation of boundary conditions can help reduce some of the errors and improve the accuracy of the simulations.However, a finer computational grid has limited practical value.Time constraints quickly prevent computation in a reasonable time frame.Instability can arise when the chosen time step is too large, causing the numerical solution to diverge or exhibit unphysical oscillations.
Noteworthy numerical method is the generalization of the Transmission Line Matrix (TLM) method from electrical circuit design.Its adaptation to modelling propagation in urban environments can be considered as a full-wave method in the time domain [38].The principle of flows implicitly models wave reflections and diffractions.The high computational costs are slightly reduced in its multi-resolution frequency-domain variant MR-FDPF (Multi-Resolution Frequency Domain Parallel Flow) [39].Simplifications in three-dimensional space by avoiding propagation modes with low impact or a combination of several two-dimensional subproblems [40] lead to approximate solutions that require calibration.
The integral form of Maxwell's equations serves as the basis for a numerical solver with very limited applications to small indoor environments in [41] and [42].The approach does not require explicit boundary conditions.Homogeneous dielectric or conductive geometries are solved by surface integrals [43], while inhomogeneous materials require the use of volume integration [8].The discretization of the equations is based either on the Method of Moments (MoM) [43] or on the hybrid Finite-Element Boundary-Integral (FE-BI) [8].The computation can be accelerated for larger repetitive geometries by Array Decomposition-Fast Multipole Method (AD-FMM) [44].The Volume Electric Field Integral Equation (VEFIE) has been proposed for three-dimensional geometries [41].
The base premise of the parabolic equation method is the paraxial approximation.The wave equation is reduced to a parabolic partial differential equation by assuming that the spatial variations of the wave are much slower compared to the variations in the propagation direction.It was introduced in [45] for a problem of radio propagation in the atmosphere.The modelled geometry should have a preferred propagation direction, with important physical phenomena not occurring at angles greater than 15 degrees to this direction, while the propagation media should be weakly inhomogeneous.Tunnels and other special geometries are best suited for the parabolic equation [46].Wider application is made possible by relaxing the 15-degree constraint [47], [48].
The PSTD (Pseudo-Spectral Time-Domain) method combines the advantages of spectral methods and finite-difference time-domain methods to provide efficient simulations [49].Spectral representation involves expanding the fields into a series of basis functions, such as Fourier series, Chebyshev polynomials, or other orthogonal functions.The choice of basis functions depends on the specific problem and the geometry of the computational domain.Spatial derivatives are evaluated exactly for at least two points per wavelength.The temporal derivatives still contain the second-order error term.Perfectly matched layer ABC is required to prevent the Fast Fourier Transform (FFT) periodic behavior to generate additional errors.However, only geometries consisting entirely of dielectrics can be modelled with PSTD, since continuity of the tangential fields is required, which is clearly violated on metal surfaces.

D. HYBRID METHODS
Hybrid methods are a common compromise between fast but less accurate ray tracing and computationally intensive numerical methods.Numerical treatment is used in the areas with complex discontinuities or in areas where the size of the details is close to the simulated wavelength, while ray tracing is used in the rest of the geometry.The transition between areas or domains is a major challenge.Moreover, user interaction is required to define the domain partitioning [41], [50], [51].Due to the complexity involved, hybrid approaches in three dimensions are rare [52], [53].

III. QUANTUM COMPUTING
The principles of quantum mechanics allow quantum computers to perform certain types of computations much more efficiently than classical computers.Here we give a brief introduction to quantum computing focusing on the different architectural approaches.
The Schrödinger equation provides the quantum mechanical framework and underlies the principles of quantum computing.In quantum computing, a qubit can be in a state of 0 or 1, or any combination of the two.This property arises from the fact that the Schrödinger equation is a linear differential equation, and if we have two solutions to the equation, any linear combination or superposition of these solutions is also a solution.Two or more qubits can be entangled, i.e., their states are connected.When two qubits are entangled, their common state is described by a single wave function.The wave-like behavior of quantum states leads to interference effects that are central to many quantum algorithms.In describing the time evolution of a quantum system, the Schrödinger equation has the form where h is Planck's constant and H is a fixed Hermitian operator known as the Hamiltonian of a closed system.The eigenvalues of H determine the energy levels of the system, with the lowest energy state called the ground state.Just as classical computers use logic gates to perform operations on bits, quantum computers use quantum gates to perform operations on qubits.Unlike classical gates, however, quantum gates are reversible and operate on qubits while they are in a superposition of states.Quantum gates are unitary transformations, and these transformations evolve the quantum state of the system according to (1), where the connection between discrete time unitary operators and the continuous time system Hamiltonian is Not all approaches to quantum computation are gate-based.Adiabatic quantum computing (AQC) can be considered a subclass or model of quantum computing.It is based on the principles of quantum mechanics but operates under a different paradigm.AQC is particularly useful for optimization problems.However, there are open questions about its universality, i.e., its ability to efficiently solve certain computational problems.The core of AQC is the adiabatic theorem, which guarantees that a system initially in its ground state will tend to stay in the lowest energy state if the Hamiltonian of the system is changed slowly enough [54].Suppose the Hamiltonian H 0 has an easily set up ground state and the Hamiltonian H P is such that its ground state is a solution to our problem.A purely adiabatic calculation can change the system Hamiltonian with time The system remains in the ground state while the measurement at time T yields a solution to the optimization problem.Quantum algorithms can be classified according to various criteria, such as the type of the problem or the complexity that can be achieved.Some quantum algorithms are designed for adiabatic quantum computers or quantum annealing systems, while others are designed for gate-based quantum computers.Moreover, some quantum algorithms are hybrid algorithms that incorporate a significant amount of classical computation in addition to quantum computation.Others, such as Shor's or Grover's algorithms, are pure quantum algorithms.Since quantum algorithms in the field of EM propagation modelling are not only gate-level algorithms, but also hybrid algorithms and algorithms developed for annealers, a brief overview of both of these approaches follows.

A. VARIATIONAL QUANTUM ALGORITHMS
In the short term, noisy quantum computers of limited scale are more likely and are considered by many researchers to be better suited for the machines of the near future than faulttolerant architectures.The time for which a quantum system remains coherent and the maximum number of consecutive unitary operations of expected real quantum computers would prohibit the use of quantum algorithms requiring thousands of qubits, deep circuits, and noise-free operation.Variational quantum algorithms (VQAs) are a class of quantum algorithms that leverage the concept of variational optimization to solve computational problems while using a limited number of qubits and having limited circuit depth [55].The central idea of VQAs is to express a quantum problem as an optimization problem.The goal is to find the minimum of an objective function, usually called a cost or energy function, by varying a set of parameters in a quantum circuit.The circuit is referred to as the ansatz.Therefore, the final algorithm is a hybrid of classical and quantum parts.
Given the widespread knowledge of machine learning algorithms (ML), one can explain variational quantum algorithms by analogy.VQA and ML have several things in common.In machine learning, models are trained on data using various learning algorithms to adjust the parameters of the model.VQAs train the parameters of the ansatz quantum circuit by iteratively optimizing them using classical optimization methods.The ansatz should have enough variational freedom to explore different parts of the Hilbert space.This allows the optimization algorithm to search for a wide range of possible solutions and increase the likelihood of finding the optimal solution.The circuit should be parameterized, where the parameters are the variables to be optimized.In both VQAs and machine learning algorithms, a cost function or performance metric is evaluated.In machine learning, the cost function can be based on metrics such as accuracy, loss functions, or likelihoods.In VQAs, the cost function represents the energy of the problem-specific quantum system, which is evaluated based on the measurement outcomes.
Notable representatives of VQAs include the Variational Quantum Eigensolver (VQE), a type of variational quantum algorithm for estimating the lowest eigenvalue and corresponding eigenvector of a given Hamiltonian, the Quantum Approximate Optimization Algorithm (QAOA), specifically designed for solving combinatorial optimization problems, the Quantum Neural Networks (QNN), which combine principles of quantum computing and classical neural networks, the Quantum Principal Component Analysis (QPCA), which is inspired by classical principal component analysis, the Quantum Boltzmann Machine (QBM), a quantum analog of classical Boltzmann machines, and others.We give an overview of the VQE as it is used to solve the waveguide mode problem presented later.

1) VARIATIONAL QUANTUM EIGENSOLVER
The VQE is a specific type of VQA designed to estimate the lowest eigenvalue and the corresponding eigenvector of a given Hamiltonian.The VQE begins by encoding the problem of interest into a Hamiltonian operator and prepares the ansatz in a trial state that approximates the eigenvector of the Hamiltonian.The choice of ansatz depends on the problem and the available resources.At each iteration, the VQE executes the parameterized circuit on a quantum computer and measures the expectation value of the observable that matches the Hamiltonian.The expectation value represents the average outcome of the measurements and provides an estimate of the energy associated with the trial state.The average value of the observable M is defined as E(M ) = ⟨ψ|M |ψ⟩, where |ψ⟩ is the state of the ansatz circuit.In practice, a large number of repeated algorithm runs is required to measure the expectation value.
To find the set of parameters for the ansatz that minimizes the energy expectation value, the VQE employs classical optimization algorithms, such as gradient-based methods.The optimizer adjusts the parameters of the ansatz iteratively based on the measured energy, aiming to converge to the lowest possible energy value.When the optimization is complete, the set of optimized parameters corresponds to an approximation of the eigenvector of the Hamiltonian with the lowest eigenvalue.
The ansatz quantum states should have some properties that help the optimizer find a solution, e.g., some overlap or similarity with the true ground state of the quantum system under study, which provides a good starting point for optimization and speeds up convergence.The initial state should be relatively easy to prepare.In some cases, domainspecific knowledge or intuition about the problem can be used to design the initial state.For example, if the problem has known symmetries or specific characteristics, the initial state can be tailored to exploit these properties, leading to improved efficiency.The problem of the optimal initial quantum state can be a challenging task and often requires a combination of theoretical analysis and trial-and-error.General ansatz topologies have been proposed, one of which, the hardware-efficient ansatz (HEA), will be presented in the section on waveguide modes callculation, which includes a practical example of how to encode the problem as an eigenvalue problem and define the corresponding measurement observable.

B. QUANTUM ANNEALING ALGORITHMS
Quantum Annealing (QA) was developed to solve combinatorial optimization problems by exploiting the adiabatic evolution of a quantum system.QA not only exploits the aforementioned adiabatic theorem, but also builds on the phenomenon of quantum tunneling, which allows particles to move between states even when separated by energy barriers.As with VAQ, the goal of QA is to find the lowest-energy state of the quantum system that corresponds to the solution of the problem.However, the way such a state is achieved is very different.
In quantum annealing, the quantum system is first initialized in a simple initial state, typically a superposition state, and then the Hamiltonian of the system is slowly changed toward the one that encodes the optimization problem.The goal is to find the ground state of the targeted Hamiltonian that corresponds to the optimal solution of the optimization problem.The main difference from the purely adiabatic calculation is that quantum annealing typically involves thermal effects and does not necessarily operate at the ground state of the system.The system can explore higher energy states before returning to a lower energy state, potentially avoiding getting stuck in local minima.Tunneling allows the quantum system to transition directly from one state to another by passing through the energy barrier rather than overcoming it.On the other hand, the transition from a local minimum in classical annealing algorithms can be very computationally intensive.
In the final step, a measurement is made on the qubits that, due to the nature of quantum measurement, collapses the superposition of states into a single state that corresponds to a possible solution to the problem.If the quantum annealing process was successful, this is the global minimum or a good approximation of it.
While it has been shown that quantum annealing can provide a speed advantage over classical algorithms for some problems, including certain optimization problems, it is still an open question whether quantum annealing can solve NP problems in polynomial time.

1) QUADRATIC UNCONSTRAINED BINARY OPTIMIZATION
An important class of problems that can be solved with QA are problems expressed as quadratic unconstrained binary optimization (QUBO).QUBO allows the mathematical formulation of a wide range of optimization problems, including binary integer programming, quadratic programming, and combinatorial optimization.In QUBO, the goal is to find the binary vector x ∈ B N that minimizes the objective function which is a quadratic function of the binary variables.The objective function is usually represented as a real-valued upper triangular matrix Q, where the elements Q ij represent the strength of the interaction between the binary variables x i and x j .
To solve a QUBO problem with quantum annealing, a conversion to the Ising model is usually performed.In the Ising model, the so-called spin variables can take the values −1 or 1.The models are closely related and equivalent in computational power.The conversion from QUBO to Ising is a simple transformation s = 2x − 1, where s is the Ising spin variable and x is the QUBO variable.The Ising formulation allows a straightforward mapping of the objective function to the Hamiltonian of a physical system.The procedure involves summing the tensor products of weighted Pauli-Z matrices acting on the individual qubits in a Hilbert space of N qubits [56].The system is then evolved over time using a quantum annealing schedule.D-Wave Systems is particularly known for their work in the field of quantum annealing [57].The company has produced a series of quantum annealing processors for solving optimization problems.These processors operate using a network of superconducting quantum bits.

IV. QUANTUM TRANSMISSION LINE MATRIX
First, we take a look at the pure quantum gate-level algorithm proposed in the context of EM propagation modelling.Like FDTD, the TLM method is an explicit time-domain method that appeared a few years after FDTD in 1971 [58], [59].Based on the transmission lines circuit model, it evaluates the propagating EM fields in analogy to the currents and voltages in electrical networks.The underlying concept is the discretized Huygens' model of wave propagation and scattering.The computational domain is represented by a mesh of orthogonal transmission lines that interconnect at the scattering nodes.There are different types of scattering nodes, e.g., the symmetrical condensed node (SCN), the symmetrical super-condensed node (SSCN), the hybrid node, etc.The electromagnetic field is represented by pulse signals propagating through the network of transmission lines.When these signals hit a node, they scatter into the other transmission lines connected to that node.The solution is stepwise in time, calculating the scattering process at each node for one time step before moving to the next.The scattered signals are related to the field quantities of interest, such as electric and magnetic fields.For example, in the version of the transverse electric wave problem (TE Z ), the relations of the quantities are where I , V , E, and H are current, voltage, electric field, and magnetic field, respectively.The scattering matrix describes the reflected signal along the transmission lines for a given node and a single time step.In free space, where the scattering event is lossless, the four incident voltages are mapped to the four reflected voltages according to where the subscript k denotes a time step, the superscript r the reflected values, and the superscript i the incident values.The reflected voltages then become incident voltages in the adjacent nodes.For given x and y coordinates of nodes, the propagation rules apply.In a more comprehensive model, boundary conditions, complex geometries, and different materials would have to be considered, which may change the scattering matrix and the resulting equations.Similar to the FDTD method, the TLM is sensitive to the propagation direction and signal frequency, while harmonic solutions require the use of the Fourier transform.The TLM method can be very computationally intensive due to the fine spatial discretization required for high frequency simulations.Therefore, it could benefit from quantum acceleration.

A. PROBLEM BEING SOLVED
In [60], the TLM method was implemented as a gate-level quantum algorithm (QTLM).The proposal relies heavily on the earlier Hilbert space formulation of the TLM [61], [62], in which the scattering process at a node is represented by a unitary operator acting on the signal vector.By representing the time-stepping process as a sequence of such operations, the TLM can be interpreted as a dynamical system evolving in a quantum space.The steps of the mesh calculation are implemented in a straightforward way.Each quantum basis amplitude corresponds to a single field value of the internode connection.Although the Hilbert space formulation can be conveniently mapped to unitary operations, decomposed by the authors into a product of universal gates, we cannot directly observe the state vector.The proposed QTLM avoids reading the field values by addressing a variant of the design problem in electromagnetics, where the task is to select an electromagnetic structure with the required electromagnetic properties corresponding to a given electromagnetic response.Such a formulation does not need actual amplitude extraction.We know in advance what the result should be, but we do not know for which structure.

B. FIELD STATE SPACE
Field amplitudes are computed in an M × N mesh of nodes, each of which having four connections to its neighbors.Four incident fields per node yield two qubits in addition to the log 2 M + log 2 N qubits needed to represent the field values of a single structure.Therefore, the field Hilbert space is described by the orthonormal basis vector set |m, n, i⟩, where m and n are the node coordinates and i is the transmission line index of the node.
For a given electromagnetic structure, one must combine scattering and connection operations and repeat them as many times as there are simulation steps.The connection operator can be decomposed into a product of qubit permutations in the x and y directions.An efficient implementation of the connection operator requires the ordering of adjacent states based on Gray code encoding to obtain consecutive addresses that differ only in one qubit value.

C. STRUCTURE STATE SPACE
In the case presented in [60], each node can be either a free space or an ideal conductor.Although two different scattering operators can be implemented, the authors choose to use the free-space scattering operator for all nodes and invert the behavior in conducting nodes with the superimposed operator.The free space scattering operator S F is given by the matrix in (6), while the scattering operator for a perfect conductor is given by the 4 × 4 matrix S C = −I .To invert a free space scattering and simulate an ideal conducting node, S F should simply be multiplied by −S F .Note that the scattering operator depends on the material of the node.
To evaluate the EM field for different structures simultaneously, a number of structure qubits must be defined, forming a so-called structure Hilbert space.Each structure qubit represents a single element of a real world, which can be either present or absent and serves as a control qubit for the above scattering operator.The number of structures grows exponentially with the number of elements that make up these structures.In the case presented here, where each element can be either a free space or an ideal conductor, the number of distinct structures doubles with each additional building block, resulting in 2 B possibilities, where B is the number of building blocks and the number of structure qubits.In the extreme case, when each node can be an independent object, we need MN structure qubits.Before they can be used in the controlled scattering gates, the qubits must be brought into an ideal superposition state by a Hadamard transform to allow simultaneous field evaluation of all configurations.Fig. 1 shows the overall architecture of the QTLM before the result is extracted.

D. RESULT EXTRACTION
To find out the structure with the preferred response, the authors propose to first measure the field qubits with a carefully designed observable O.One eigenvector of the observable, e.g., |e 1 ⟩, should correspond to the expected field solution.The other 4MN − 1 eigenvectors are orthogonal to the first eigenvector and have distinct nondegenerate eigenvalues If the measured eigenvalue matches the first eigenvalue belonging to the desired electromagnetic response, the structure qubits collapse to a structure description with high probability that it is the structure that led to the above response.Measuring the structure qubits in the computational basis can then possibly provide the solution.If a different eigenvalue is measured, the result is discarded.It should be noted that the correctness of the obtained solution should be verified by a classical TLM.It is assumed that this can be checked quickly, but the number of structures to be checked is too high to be handled by classical computers.The measurement step is shown in Fig. 2.
The authors acknowledge that the probability of measuring the required eigenvalue and obtaining the correct In [63], nonlinear quantum mechanics [64] and its profound implications for computational theory are discussed in the context of improving QTLM.Basically, the approach would amplify the amplitudes of the desired field states.However, nonlinear mechanics is highly speculative and physically infeasible at the current state of quantum architectures.

V. QUANTUM SIMULATION OF WAVE EQUATION
EM propagation can be described by a single wave equation using either an electric or a magnetic field, rather than a set of coupled differential equations as is the case with the FDTD algorithm.In a source-free region, the electric field equation is of the form where c is the speed of light in free space and ∇ 2 is a differential operator over a vector field, i.e., a vector Laplacian.In addition to the propagation of electromagnetic waves, the general wave equation has many important applications in physics, engineering, and other fields since it relates the spatial and temporal variations of a wave to its frequency and velocity.

A. MAPPING TO A QUANTUM SYSTEM
Costa et al. [65] base their solution of a quantum wave equation simulation on the fundamental time evolution of a quantum system, which naturally relates the wave function of a system to its energy in the form of the Schrödinger equation ( 1).In the proposed solution, the finite domain of a simulated space is discretized on a cubic lattice grid with spacing a.The grid is then represented as a graph with a set of points and edges describing the neighborhood relation 111552 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
between the grid points.The Laplacian matrix of such a graph is defined as L = D − A, where D is the degree matrix and A is the adjacency matrix.The entry in the i-th row and j-th column of the adjacency matrix is 1 if there is an edge between points i and j, and 0 otherwise.The usefulness of the Laplacian matrix lies in its property that in the limiting case, when the lattice spacing a approaches 0, the −L/a 2 approximates ∇ 2 .The discretized wave equation is then Here is used to represent a column vector of field values, and the wave propagation velocity is taken to be 1.
To evaluate (10) by evolving quantum system described by (1), the following quantum system is proposed.Based on the lattice graph, let the matrix B be w j if j is an edge with i as source, − √ w j if j is an edge with i as sink, 0 otherwise.(11) where the weights w j are set to 1, except for the weights due to the boundary conditions, which will be detailed shortly.
Setting the Hamiltonian of a quantum system to and the Hilbert space associated with the quantum system to a direct sum of vertex and edge spaces H = H V ⊕ H E , we can show that BB † = L and, by taking the second derivative, the amplitudes of a subspace H V in behave according to (10).
To solve practical scattering problems, one has to include different boundary conditions into the computational algorithm.Usually, Dirichlet and Neumann boundary conditions are used.While the former specifies the value of the solution at the boundary of the scattering region, the latter specifies the derivative.In simulating the wave equation, only changes to the graph Laplacian are required to satisfy these two conditions.By removing vertices within the scattering region, the Neumann boundary condition is realized.On the other hand, weighted self-loops should be added to the boundary vertices to implement the Dirichlet boundary condition, where the weights are the number of missing edges to the removed interior vertices.

B. COMPLEXITY
The authors do not provide a gate-level implementation of the time evolution e −iHt , but estimate time and space complexity based on the properties of general sparse Hamiltonians [66]

C. INITIAL STATE AND RESULT EXTRACTION
The initial state [ V , E ] can be prepared in polynomial time if the quantum representation x w(x)|x⟩ of a discretized twice-differentiable initial wave packet w(x − ct) can be prepared in polynomial time.In the general case, the authors discuss a scenario where the initial state is prepared using the quantum algorithm for solving a system of linear equations [67].
The final quantum state encodes the intensity of the wave as the lattice field values (T ) and the values associated with the edges B −1 ∂ (T )/∂t.No sophisticated procedure is provided for measuring the field values, so each measurement can only provide a sample of the distribution of state probabilities, which is proportional to the square of the field amplitude.

VI. VQA FOR WAVEGUIDE MODES
The propagation of EM waves in a waveguide exhibits diverse field patterns and unique cutoff frequencies.Knowledge of these characteristics is necessary for efficient transmission, filtering and other applications.Propagation modes depend on both the geometry of the waveguide and the operating frequency.In addition to the fundamental mode, higher modes with more complex field patterns and typically higher cutoff frequencies are also of interest.

A. PROBLEM SPECIFICATION
In [68], the finite difference method is proposed for the calculation of propagation modes in a hollow metallic waveguide using VQE, a variant of VQA.Transverse electric (TE) and transverse magnetic (TM) modes are considered as two separate problems.The Helmholtz equations describe the propagation of EM waves in waveguides and are used to calculate the modes or solutions for the electric and magnetic fields within the waveguide.They are derived from Maxwell's equations and take the form of partial differential equations.2D scalar wave equations are where ∇ 2 is the Laplacian operator, k is the wavenumber, and E and H are the electric and magnetic field components, respectively.The solutions of the Helmholtz equations yield the modes or eigenfunctions of the waveguide, which represent the various possible configurations in which electromagnetic waves can propagate within the waveguide.

B. VQE FORMULATION
The finite difference method starts with a discretization of the rectangular cross-section of a waveguide with a uniform shifted grid.The Laplacian operator in the Helmholtz equations of the TE and TM waves simplifies to the second-order derivative approximated by the central finite difference.Since the x-and y-directions can be considered as two orthogonal problems, the quantum state space can be a tensor product of the x-and y-Hilbert space, with the waveguide axis coincident with the z-direction, as shown in Fig. 3.The Helmholtz equations can then be transformed into observables for the TM mode and for the TE mode, where M {x,y},{D,N } are Laplacian matrices for the x and y directions and either the Dirichlet boundary condition or the Neumann boundary condition are considered for the TM and TE modes, respectively.The matrices have similarity to the Laplacian matrix used by Costa et.al [65] for the lattice in the quantum wave equation simulation.
The propagation modes are then the result of an eigenvalue problem where M corresponds to either the M TE or the M TM , λ is an eigenvalue and ν is an eigenvector.The lowest propagation mode corresponds to the lowest eigenvalue.To solve the problem with the variational quantum eigensolver as presented earlier, a suitable cost function must be defined.The basic value of the parameterized ansatz |ψ(θ)⟩ would iteratively converge towards the eigenvalue of the fundamental mode.A penalty term is needed in the cost function to guide the variational quantum algorithm to higher order eigenvalues where β i should be any constant satisfying β i > λ k − λ i , and θ i is the value of the parameter θ that approximately minimizes the cost function of lower mode i.The optimization step in the proposed algorithm is performed by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimizer, which is based on gradient descent.

C. PARAMETERIZED ANSATZ
The parameterized ansatz uses a combination of a single-qubit rotations and two-qubit entanglement gates.The rotations can be adjusted by parameters, allowing for variation and optimization during the algorithm.The circuit is known as hardware-efficient ansatz (HEA).In this case, n layers were used, with an initial state |0⟩ ⊗n as shown in Fig. 4, where n = n x + n y .The HEA was specifically designed to address the limitations and constraints of current quantum hardware, such as limited qubit connectivity and the inability to execute certain types of quantum gates.These gates are chosen to be intrinsic to the quantum hardware, meaning they can be implemented directly without the need for complex decompositions or additional gate operations.

D. COMPLEXITY
The decomposition of M into simpler observables for implementation can, in principle, be done using the Pauli basis.However, the number of terms in such a decomposition increases exponentially with n.Therefore, the authors express M as a linear combination of simple Hamiltonians [69] with a constant number of terms.Although the overall complexity of the number of quantum gates is not specified, the fact that only a constant number of expectation values need to be evaluated independently of n in each iteration step should provide the quantum advantage over the classical algorithm for sufficiently large n.

E. EXAMPLE
The example waveguide presented in [65] is 15 mm wide and 10 mm high with 16 and 8 discretization points, respectively, i.e., n x = 4 and n y = 3 qubits.The quantum part of the VQA was simulated using the Qiskit simulator from IBM [70].The cutoff frequency agrees well with the classical numerical solution and the analytical solution, with an error below 0.001%.

VII. QUANTUM RAY TRACING
Ray tracing seems to be less suitable for quantum computing when it comes to computing field, but the inherent parallelism of geometric subproblems can still benefit from the quantum advantage, as shown in [71].

A. PROBLEM BEING SOLVED
In general, ray tracing consists of a systematic enumeration of ray paths between known transmitters and receivers based on the position mirroring technique or some other pathfinding technique.Common to all approaches is the problem of determining visibility between two points in a model of the world described by a collection of geometric primitives.
The algorithm must establish whether the direction of the ray is occluded and which geometric primitive intersects the ray path before the path is reflected, scattered, or diffracted.
A similar problem is solved in rendering algorithms where ray tracing is used to color a pixel instead of computing an electric field.
If we were to naively check every object in the scene for intersection with a ray, we would have to perform an intersection test for each object, leading to a time complexity of O(N ), where N is the number of objects.The proposal in [71] is a very basic attempt to reduce this complexity to O( √ N ) by using the well-known Grover's quantum search algorithm for the ray occlusion and intersection problems.Note that O(N ) can be improved ideally to O(log N ) in the classical case by using a spatial acceleration technique such as a bounding volume hierarchy (BVH) [72].However, the time required to create the BVH in the first place is typically O(N log N ).
The occlusion problem in [71] is extremely simplified and therefore serves only as a proof of principle.The rays are assumed to be parallel to the z-axis, which corresponds to the viewing plane z = 0 in the orthographic projection.The geometric primitives all have the same distance from the viewing plane, they do not overlap, and they have a rectangular shape.All coordinates are treated as positive integers.is only one possible solution, namely the primitive, which can block the path of a given ray.It should be noted that knowing the number of solutions is important for choosing the number of iterations in the basic quantum search algorithm.

B. GROVER SEARCH
The indices of N = 2 n geometric primitives are encoded as computational basis states |ψ⟩ of a state space formed by n qubits.The occlusion search algorithm should define an oracle required in the Grover algorithm that adds a negative phase to the solution states.In our case, the solution state is the one that represents the intersected object for a given ray.Note that it is easy to check the intersection between a ray and an object, i.e., to compute the oracle, but time-consuming to check all the intersections individually.Quantum parallelism allows simultaneous intersection tests for all geometric primitives by bringing the states |ψ⟩ into a uniform superposition.Measuring such a state space without the Grover algorithm in the standard basis would result in the state collapsing into any of the basis states with equal probability.The Grover oracle, followed by a diffusion operator, is able to amplify the amplitude of the solution state, which increases the chances of measuring the identification of the intersected object.
The Grover search algorithm begins with an equally weighted superposition of states |ψ⟩.The iteration step of the algorithm is (2|ψ⟩⟨ψ|−I )O, where I is the identity matrix and O is our oracle.If the vector |ψ⟩ is expressed in the plane spanned by two orthonormal vectors, the first of which is the normalized sum of all computational states that do not represent our solution, and the other is a normalized sum of all solutions to the search problem, it can be shown that the oracle operator reflects the current state on such a plane with the reference to the nonsolutions basis vector while the diffusion operator, i.e., 2|ψ⟩⟨ψ| − I , reflects the current state around |ψ⟩ on the same plane.Thus, the Grover step rotates the quantum state to the basis vector formed by all solutions.Measuring the first n qubits will then lead to one of the solutions with high probability.The probability of finding the solution is brought close to 1 by choosing the optimal number of iterations.This depends on how close ⌈ 4 π √ N /t⌉ is to the integer value, where t is the known number of solutions, i.e., one in the nonoverlapping occlusion test.
The oracle implementation requires additional qubits to test the limits of a primitive encoded in |ψ⟩.In addition to some extra qubits required for intermediate results and reversibility of the computation, ⌈log 2 (max(max x , max y ))⌉ qubits of register |c⟩ are needed to store the minimum and maximum coordinates of a primitive in both dimensions.The computation of the occlusion condition is performed sequentially, with the register |c⟩ being uncomputed between steps, and the boundaries of the primitives must be hard coded in the unitary transforms used in the oracle calculation.

C. SINGLE SOLUTION INTERSECTION TESTING
The proposed hybrid algorithm runs the quantum code for each ray until an intersection is found or the iteration limit is reached.The intersection is verified by the classical code.The algorithm has been successfully tested on IBM's quantum simulator [70].The example scene with 8 primitives required 195 quantum gates acting on 15 qubits.However, running the examples on a real IBM Q Network machine was not successful, except for the trivial examples, because the length of the longest gate sequence of the proposed circuits, e.g., 68 for the above scenario, is too long for the current real quantum machine.Even for the extremely trivial 4-primitive example, the results were so noisy that the probability of correct reading was only 0.53.

D. MULTIPLE SOLUTION EXTENSION
The algorithm has been extended for the case where primitives are at different depths and may overlap from the ray perspective.In addition to extending the oracle calculation to handle z-coordinates in boundary extraction and intersection checking, the unknown number of solutions must also be handled appropriately in the classical part of the hybrid algorithm.Not knowing the number of intersections in advance poses a problem for the basic quantum search algorithm, since this information is needed to determine the optimal number of iterations and to bring the solution amplitude close to 1.In the proposed extension, the exponential search [73] and the quantum algorithm for finding the minimum [74] were combined with the hybrid algorithm.In the proposed algorithm, the oracle circuit should be adjusted at each iteration.

VIII. QUANTUM SHORTEST RAY-PATH
In [75] the imaging method, a variant of ray tracing, is approached with a quantum annealer.Given a number of interactions N and a set of M interaction objects, the problem to be solved is to find the dominant propagation path, i.e., the path that experiences exactly N reflections or diffractions and has minimal propagation loss (PL).Furthermore, the goal is to find multiple ray-paths in ascending order with respect to PL.The number of combinations to be searched for the above single path problem is M (M − 1) N , which determines the complexity of the classical algorithm.
In order to formulate the task as a QUBO problem, further simplifications have been made which make the applicability of the solution to the real problem difficult and, similar to the earlier quantum search of occlusions, make the solution only of theoretical interest.The first assumption is that the radio wave is scattered uniformly in space after each interaction, either reflection or diffraction.Ignoring Snell's law for reflections and Keller's cone attenuation for diffraction is a major departure from reality, and it is not possible to see how the QUBO formulation can be improved in this respect.Moreover, the path lengths of the rays between the interaction points are assumed to correspond to the distances between the objects, which may be justified for small and distant objects, but the phase shift information is certainly lost.Making all the above changes, the problem of finding a shortest path with a minimum PL given N interactions becomes a version of the problem of finding the shortest path in a graph.Let the binary variable x i,k encode the solution path by taking the value 1 if the object m i is the k-th interaction point, and 0 otherwise.The shortest path in QUBO form is then the one that minimizes the distance expression where and δ i,j is the Kronecker delta and d i,j is the distance between objects m i and m j .In addition, two constraints should be included in the QUBO formulation (20) a) interactions with multiple objects at the same time are forbidden and b) the interactions cannot be successive with the same object.Both constraints must be in QUBO form.For example, the absence of path loops in the latter constraint has the form The expression would be 0 if the constraint is satisfied, and a positive value otherwise, causing the solver to choose a path without loops.With this formulation, QA would search for a single dominant path on which N interactions occur.Searching for other paths in ascending order with respect to path length requires a sequential approach, where the objective function is modified to exclude paths already found.Therefore, the final Hamiltonian coding QUBO for the quantum shortest raypath (QSRP) encodes a real-valued matrix Q (4) composed of the weighted sum of the distance expression, two path shape constraints, and the expression to exclude already found paths with shorter lengths.
The authors compare simulated annealing in software for a problem with M = 9 objects and N = 3 interactions with the algorithm running on Fujitsu Digital Annealer [76].Apart from finding that the physical realization of the annealer is faster than the simulated annealing, no further complexity analysis was performed in comparison to classical algorithms.The digital annealer operates probabilistically and cannot always find the optimum, which was the case for the 9th and 10th paths in descending order in the example presented.It should be noted that the digital annealer is a classical device that does not benefit from quantum effects such as superposition and quantum tunneling.It is a specially designed classical hardware to accelerate annealing simulations.

IX. QUANTUM FINITE-ELEMENT METHOD
The finite element method (FEM) is widely used in EM wave propagation problems as well as in many other fields.The method is a powerful computational tool for solving a variety of engineering and scientific problems.In telecommunications, FEM can predict the behavior of electromagnetic waves in the high-frequency range and can be used to analyze wave propagation in waveguides, optical fibers, and other guiding structures.FEM simulations provide insight into propagation modes, dispersion, and losses that are critical to the development of modern communications systems and optical devices.The method helps predict antenna radiation patterns, electromagnetic fields of microwave circuits, electromagnetic interference, and it used to solve many other problems.
In the FEM, the continuous partial differential equations (PDEs) describing the problem are discretized into a system of algebraic equations.This is usually done by applying variational principles or Galerkin methods to obtain a set of equations that relate the unknowns, such as nodal values or element-based field representations, to the known parameters of the problem.These equations usually take the form of a sparse linear system, represented as where K is a coefficient matrix, e is the vector of unknowns, and c is the right-hand side vector.Solving this linear system gives the solution of the original problem [7].It is important to note that the HHL algorithm is highly sensitive to errors and noise in quantum systems.Implementing and executing the HHL algorithm on current quantum hardware remains a major challenge due to the requirements for high coherence and precise gate operations.Further, the exponential speedup in solving a FEM problem has been questioned when predetermined accuracy is required [78].HHL has some limitations that prevent a straightforward solution of (24).First, in order to find x that solves Ax = b with HHL, the matrix A must be Hermitian and the vector b should be represented as a quantum state |b⟩.For the HHL algorithm to achieve exponential speedup, the matrix A must not only be sparse and prepared efficiently, but also have certain properties, such as a small condition number, since the running time is proportional to the square of this number.Further, the state |b⟩ must also be efficiently preparable.The solution A −1  |b⟩ is encoded in the final quantum state, which cannot be easily extracted.To achieve exponential speedup, the output of the HHL algorithm should be the expectation value ⟨x|M |x⟩ for a user-defined operator M specific to the problem at hand.
The HHL algorithm, shown in Fig. 5, starts with a vector b encoded as quantum state |b⟩.By applying Quantum Phase Estimation (QPE) to a unitary operator e iAt for a superposition of different times, the HHL algorithm generates a quantum state that can be viewed as a decomposition of |b⟩ in the eigenbasis |u i ⟩ of A and extended by the eigenvalue phase information, i.e., i β i |u i ⟩|λ i ⟩, where β i are the decomposition coefficients and λ i are the corresponding eigenvalues.A controlled rotation based on these eigenvalues is performed next on an ancillary qubit originally set to 0. The operation roughly translates |λ i ⟩ into Cλ −1 i |λ i ⟩, where C is a normalization factor, and effectively maps the inverse eigenvalues of the matrix to the amplitudes of the quantum state.The final step is to uncompute the quantum state to its original state before phase estimation by applying the inverse of QPE.The ancillary qubit is then measured, and if the result is 1, the final state of the remaining system is proportional to A −1  |b⟩, which is the encoded solution of the system of linear equations.If the measurement of the ancillary qubit yields a value of 0, the entire process must be repeated.Note that as part of the phase estimation step, i.e., the operator U in Fig. 5, a Hamiltonian simulation is required, which may not have efficient implementation.
The generalization of HHL to arbitrary problem specifications is given in [79].The HHL does not give instructions on how to efficiently prepare the state |b⟩.Clader et.al [79] propose a modified initial state that is entangled with the ancilla qubit and can be prepared efficiently if the computation of the individual amplitude and phase components of |b⟩ is efficient.Next, they extend the original HHL with two additional ancilla qubits and show how to extract the entangled solution |x⟩ for three different problems without a post-selection step using the ancilla qubit measurement: overlapping a solution with an arbitrary vector, computing moments of the solution ⟨x|x n |x⟩, and extracting individual values of the solution ⟨j||x⟩.Note that accessing the entire solution in exponentially large space would negate any speed advantage.Furthermore, since the Hamiltonian simulation in the original HHL scales linearly with the condition number, they propose a preconditioning procedure.Instead of solving Ax = b, they solve the modified system MAx = M b.By finding MA with a low condition number, they preserve the speedup for a larger number of problems.A class of preconditioners known as sparse approximate inverse preconditioners has been proposed and integrated into the algorithm for general linear systems.

B. QFEM ALGORITHM
Research papers on quantum FEM are numerous, but few have presented an example that can be considered close to a propagation topic.Worth mentioning is [78], which deals with the solution of a Poisson equation used in problems of electrostatics.The importance of this work lies in the conclusion that in order to achieve predetermined accuracy ϵ the quantum FEM can only achieve a polynomial speedup.Moreover, in case of low dimensionality and sufficiently smooth solution, the running time of the quantum algorithm can be even worse than that of the classical algorithm.The exponential speedup does not vanish if the quantum state is the desired result and one does not evaluate a property of the state by measurements.
In [79], an EM scattering cross section is considered.The authors show that the solution using a preconditioner has quadratically better complexity in the condition number than the original HHL, since post-selection on the ancilla qubit is omitted.Note, however, that [78] claims that [79] does not fully incorporate the accuracy parameter in the run time analysis.
The QFEM variant in [80] addresses the necessary determination of hyperparameters for the application of the original HHL.The transformation of a nonhermitian matrix into a Hermitian one has already been proposed in [77], in [80] a procedure for integer powers of 2 of rows and columns in the augmented finite element matrix A is proposed.In addition to using the complex conjugate transpose of K , the transformation requires the diagonal matrix where λ min is the smallest positive eigenvalue of K .The vector b that enters the HHL algorithm should then be b = [c † , c T 0, 0, . . ., 0, 0].
The size L of a register containing the phase information is given by where λ max is the largest positive eigenvalue of A. The appropriate value of the parameter t of the Hamiltonian simulation where λmin = ⌈(t min λ max 2 L−1 )/π⌉ (30) and After completing the HHL on the initial state, which consists of one ancilla qubit, L working qubits, and ⌈log 2 M ⌉ + 1 I/O qubits, where M is the number of rows of K , the I/O register is proportional to the EM field.The solution of a 2D electrostatic problem with 200 cells was verified using Cirq [81], a quantum simulator supported by Google.

X. QUANTUM METHOD OF MOMENTS
The Method of Moments (MoM), also known as the Boundary Element Method (BEM), is a numerical technique that focuses on solving integral equations derived from the underlying physics of the problem.In MoM, the problem domain is divided into surfaces or boundaries, and the integral equations are formulated based on the electric and magnetic field relationships at these boundaries.The MoM discretizes the integral equations into a system of linear equations, similar to the FEM discretization of partial differential equations.The integral equations relate the unknown quantities, such as surface or volume currents, to the known field quantities or sources.These equations are solved to obtain the unknown current distribution or other field quantities.MoM is particularly suitable for solving problems with open regions or unbounded domains, such as radiation and scattering problems.It has already been applied to propagation problems, e.g., in small indoor environments [46], [47].Since MoM involves solving a set of linear equations as part of the numerical solution process, the quantum HHL algorithm can be used in a similar way as for FEM.
A Quantum MoM (QMoM) algorithm applied to an EM problem can be found in [82], where an example consisting of electrostatic integral equations for a 2D cross section of a rectangular two-conductor transmission line using Cirq is presented.On the other hand, a variant of QMoM in [83] uses preconditioning to reduce condition number and dense Hamiltonian simulation [84] while the inverse of A is computed using the Fourier transform as proposed in [67].The method is applied to problems of relative surface current density and EM radar cross section.However, both publications provide very few details on the examples.

XI. DISCUSSION
In the following, we discuss the main features of the presented algorithms, including their practical potential and the achievable speedup compared to classical algorithms.Tables 1 and 2 summarize classes of the algorithms, their quantum advantages, strengths, and weaknesses.Since the preparation of initial states and the extraction of results can affect the algorithm performance, Table 3 gives an overview of the I/O techniques as well as the asymptotic behavior of the algorithms in time and space.From what has been presented, it is clear that EM propagation modelling can benefit greatly from quantum advantage, however, there are a number of obstacles and opportunities that should be explored in the future.

A. ALGORITHM CLASSES
The main finding of the survey is that there are conceptually very different ways to perform EM propagation modelling using quantum computers, from mainstream gate-level algorithms and quantum circuits, based on the assumption of the availability of fault-tolerant architectures (QTLM, QWE, QRT, QFEM, QMoM), to variational quantum algorithms relying on the near-term and less ideal quantum computers with limited capabilities (VQAWM), adiabatic and annealing approaches (QSRP).For the latter, quantum advantage is still the subject of ongoing research, but real machines with several thousand qubits already exist.While D-Wave's quantum annealers have demonstrated quantum effects, the extent to which they can provide quantum speedup over classical computers is not yet known.The problem-specific nature of quantum annealing also means that its advantage only applies to certain types of problems, particularly certain combinatorial optimization problems, and may not VOLUME 11, 2023 111559 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

be as broadly applicable as other classes of quantum
On the other hand, VQAs offer a hybrid approach that can potentially outperform both classical algorithms and pure quantum circuits.VQAs can be designed to leverage the specific capabilities and limitations of different quantum hardware, making them applicable to a wide range of quantum computing technologies.Although the quantum advantage has not yet been conclusively demonstrated for many problems, ongoing research and development efforts are aimed at identifying problem areas where VQAs can provide significant speedups or better solutions.
Last but not least, pure quantum gate circuits have advantages in scenarios where specific quantum gates and operations must be applied precisely.Pure gate-level algorithms are often used when exact unitary transformations or specific quantum algorithms such as the Grover algorithm are required (QRT).
All these approaches have different strengths and weaknesses.The choice of architecture depends on the nature of the problem, the available quantum resources, and the specific computational requirements.In Fig. 6, various classes are shown schematically, including the reviewed algorithms and the most common algorithm building blocks.

B. QUANTUM ADVANTAGE
The studied algorithms can be applied to very different problems in EM propagation modelling and are hardly comparable in terms of speedup.In particular, QRT and QSRP are of limited value because the problem simplifications are simply too extreme, making their extension to practical cases difficult to imagine.Moreover, while QRT is conceptually interesting, it is not practical: there are classical algorithms that work faster than searching through all objects.On the other hand, QSRP does not guarantee a solution even if one exists.While QSRP QUBO is well-defined, QUBO problems in general can be difficult to solve for large domains, and the performance of quantum annealing on QUBO problems is still an active area of research.
QTLM is actually a variant of classical TLM that singles out the structure with the desired EM response, since reconstructing the actual field values would negate the advantage of field computation in exponentially growing state space.A similar problem arises with HHL, which is part of QFEM.The HHL algorithm is required to output an expectation value for a problem-specific observable, rather than the actual solution to the linear system problem, to achieve the exponential speedup.
QWE, which exploits the similarity between the simulated wave equation and the fundamental time evolution of a quantum system, is an excellent solution but does not provide an efficient method for extracting the results.
On the other hand, VQAWM does not suffer from complexity penalty due to algorithm input or output, as explained below.However, the speedup over the classical algorithm is only guaranteed by an efficient decomposition of the problem matrix into simpler observables regardless of the problem size.

C. HANDLING INPUT AND OUTPUT
Quantum algorithms generally perform well when they execute operations in Hilbert space, i.e., quantum state space.111560 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Weaknesses occur at the interface between the classical and quantum worlds, either at the stage of preparing an input or at the stage of obtaining the result.Quantum systems can have a huge number of basis states with respect to the number of input qubits, which cannot even be stored by classical computers except for smaller problems, and preparation of an arbitrary state requires control of the amplitudes and phases of each basis state (QTLM, QWE, QFEM, QMoM).
As the size of the quantum system increases, the number of parameters required to specify the state grows exponentially, making it increasingly difficult to accomplish the task.On the other hand, the preparation of the initial quantum state, which involves superposition and entanglement between qubits in the ground state (QSRP), is much less demanding.
For Grover-based algorithms such as QRT, setting up the oracle could easily overshadow the run time complexity.Note that the oracle for QRT should encode geometric primitives including their boundaries.On the contrary, VQAs and especially VQAWM have the advantage of using efficient ansatz.Moreover, the final set of optimized parameters in VQAs already corresponds to an approximation to the eigenvector of the problem Hamiltonian, i.e., the solution.
The input/output difficulties disappear when the quantum algorithm is used as a subroutine of another quantum algorithm or only the reduced result in terms of an expectation value is needed (HHL).Reading a complete set of quantum state amplitudes is time consuming, requires numerous repetitions of the entire algorithm, and can easily negate the speedup potential (QTLM, QWE).

D. PROSPECTIVE TECHNIQUES AND FUTURE DIRECTIONS
The algorithms studied differ in the details given for the actual implementation.For example, the use of Hamiltonian simulation to describe the algorithm is elegant but generally difficult to implement (QWE, VQAWM, HHL, QFEM, QMoM).Numerical methods based on the Schrödinger equation are often used to simulate the time evolution of a quantum system governed by a Hamiltonian.Methods can vary depending on the specific problem and level of complexity.Common techniques include matrix exponentiation, diagonalization, and Trotterization.Running a Hamiltonian simulation on a real quantum computer can be challenging due to the limitations of current quantum hardware, such as limited coherence times, high error rates, and limited qubit connectivity.One way to address the problem is to use universal quantum circuits for sparse and dense Hamiltonian simulations, which have been proposed in [66], [84], [85], and [86].
Modelling propagation by numerical methods typically requires partial differentiation in discretized domains.Graph Laplacians are extremely useful in this context, since they are related to the first-order approximation of the continuous Laplacian operator, except for a multiplicative constant, and can also be used to approximate higher-order derivatives.The matrix form can be used either as part of an observable (WQAWM) or a system Hamiltonian (QWE).The square matrix encodes the local relationships between neighboring nodes in the graph, e.g., points of the discretized space in propagation modelling algorithms.Moreover, the matrix is always symmetric, and its eigenvalues are nonnegative.The number of zero eigenvalues of the graph Laplacian is equal to the number of connected components in the graph and various boundary conditions can be integrated into the matrix.
Most algorithms involve some form of the eigenvalue and eigenvector decomposition problem, either explicitly (HHL, VQAWM) or as a fundamental behavior of any quantum measurement.The manipulation of eigenvectors is central to the design and success of many quantum algorithms.The eigenvalues of the Hamiltonian operator represent the possible energy levels of a quantum system.Understanding and manipulating these energy levels and the corresponding eigenvectors are fundamental to studying the dynamics and behavior of quantum systems.Eigenstates are often used as starting points for quantum algorithms.Many quantum algorithms, such as the VQE and the QPE, revolve around finding certain eigenvalues and eigenvectors of interest.These algorithms exploit the properties of eigenvalues and eigenvectors to solve optimization problems, simulate quantum systems, or estimate properties such as energies or phases.
The exponentially growing quantum state space could be useful for other open problems in EM propagation modelling.Finer spatial and temporal discretization helps improve the stability of the FDTD method.Smaller steps ensure that fast fluctuations in the electromagnetic field are properly resolved, which prevents numerical instability.However, finer discretization is associated with higher computational costs.It is often necessary to perform convergence studies to determine an appropriate level of discretization that represents a good tradeoff between accuracy and computational cost.This increase in computational resources could potentially be handled more efficiently in quantum space.Another example is the parabolic equation method for radio propagation problems with a preferred propagation direction, which could potentially benefit from a quantum approach since there are some similarities between the parabolic and Schrödinger equations.

XII. CONCLUSION
Although the last decade has seen the introduction of numerical methods in classical telecommunication modelling that extend beyond the limit of a few hundred wavelengths, the extraordinary computational effort and numerical artifacts have limited the proposals mainly to two-dimensional geometries with simplifications.Quantum computing seems ideally suited for numerical wave-based methods simply because of the underlying physical implementations of quantum information processing.The wave nature of the Schrödinger equation is a central feature of quantum mechanics and plays a key role in understanding the behavior of particles at the microscopic level.
In this paper, we have highlighted possible techniques for a broader transition to quantum algorithms in the field of EM propagation modelling.In this context, it is shown that graph Laplacians are ideal for finite-differential approximation of discretized domains in Hilbert space, where they can also handle boundary conditions.Next, finding equivalent eigenvalue problems proves beneficial for quantum implementation, since eigenvalues and eigenvectors are fundamentally intertwined with the principles of quantum mechanics.Algorithms can be well defined without having to go to the level of quantum gates and circuits.System Hamiltonians are particularly useful, although implementation can be difficult.To exploit the principles of quantum annealing, a problem specification in terms of a quadratic binary unconstrained optimization is sufficient.In general, exponential or quadratic speedup is expected for the quantum kernels of the algorithm, while input preparation and result extraction remain problematic, as is the case with many other quantum algorithms.Among the most useful building blocks for propagation modelling problems are the well-known quantum phase estimation, Fourier transform, and unconstrained Grover search, to name a few.

FIGURE 1 .
FIGURE 1.The architecture of QTLM before the extraction of result.

FIGURE 2 .
FIGURE 2. The extraction of the structure with the required EM response in QTLM.structure may be very low.Repeating the calculation until the correct solution is found could negate the quantum advantage.No case study is presented in this regard.The number of CNOT gates required for the connection operator scales with O(M +N ) in addition to O(M log M +N log N ) of TOFFOLI gates.The scattering operator, on the other hand, requires O(B log(MN )) CNOT gates, O(B log(MN )) TOFFOLI gates, and O(B) single qubit gates.The total number of qubits required is log 2 M + log 2 N + 2 + B.In[63], nonlinear quantum mechanics[64] and its profound implications for computational theory are discussed in the context of improving QTLM.Basically, the approach would amplify the amplitudes of the desired field states.However, nonlinear mechanics is highly speculative and physically infeasible at the current state of quantum architectures.
. The time complexity without logarithmic factors is approximated to Õ(tD 2 /a), where D represents the dimensionality of the problem.The space complexity in terms of required gates is estimated as D log(l/a), where l is the side length of a simulated cubic domain.The number of qubits required is log 2 [(1 + D)(l/a) D ], since the cardinality of the vertex and edge spaces is |V | = (l/a) D and |E| = D(l/a) D , respectively.

FIGURE 3 .
FIGURE 3. Alignment of a rectangular waveguide in space.

FIGURE 4 .
FIGURE 4. Quantum circuit of the hardware-efficient ansatz used in VQAWM.

A
. HARROW-HASSIDIM-LLOYD ALGORITHM Several variants of Quantum FEM (QFEM) are based on the well-known Harrow-Hassidim-Lloyd (HHL) quantum algorithm for solving sets of linear equations, its extensions and generalizations.HHL was proposed by Harrow et al. in 2009 [77].It can potentially offer an exponential speedup compared to classical algorithms for certain types of problems.

FIGURE 5 .
FIGURE 5.The architecture of the Harrow-Hassidim-Lloyd algorithm for solving systems of linear equations.

FIGURE 6 .
FIGURE 6.Quantum architectures, reviewed algorithms and the most common building blocks with potential for EM wave propagation modelling.

TABLE 1 .
Reviewed Algorithms potential in EM propagation modelling.

TABLE 2 .
Reviewed Algorithms potential in EM propagation modelling -continued.

TABLE 3 .
I/O techniques and asymptotic behavior.