# Rapid Design of Wide-Area Heterogeneous Electromagnetic Metasurfaces beyond
the Unit-Cell Approximation

Krupali D. Donda 1 and Ravi S. Hegde 1,*

###### Abstract

We propose a novel numerical approach for the optimal design of wide-area
heterogeneous electromagnetic metasurfaces beyond the conventionally used
“unit-cell” approximation. The proposed method exploits the combination of
Rigorous Coupled Wave Analysis (RCWA) and global optimization techniques (two
evolutionary algorithms namely the Genetic Algorithm (GA) and a modified form
of the Artificial Bee Colony (ABC with memetic search phase method) are
considered). As a specific example, we consider the design of beam deflectors
using all-dielectric nanoantennae for operation in the visible wavelength
region; beam deflectors can serve as building blocks for other more
complicated devices like metalenses. Compared to previous reports using local
optimization approaches our approach improves device efficiency; transmission
efficiency is especially improved for wide deflection angle beam deflectors.
The ABC method with memetic search phase is also an improvement over the more
commonly used GA as it reaches similar efficiency levels with upto 35%
reduction in computation time. The method described here is of interest for
the rapid design of a wide variety of electromagnetic metasurfaces
irrespective of their operational wavelength.

Vol. x, y–z, 2014

111Received date222* Corresponding author: Ravi Hegde
(hegder@iitgn.ac.in).3331 Department of Electrical Engineering, Indian
Institute of Technology-Gandhinagar, Gandhinagar, Gujarat, India 382355.

Krupali D. Donda and Ravi S. Hegde

##  1 Introduction

The metasurface [1, 2, 3, 4, 5, 6] is the two-dimensional analogue of the
metamaterial. It is a spatially heterogeneous array of nanoscale resonant
elements (called meta-atoms) that can, in general, alter the amplitude, phase,
spectrum and polarization values of an incident wave-front in a very short
propagation distance and with sub-wavelength resolution in the transverse
plane [6]. Being a structured surface, it can be fabricated easily in relation
to the metamaterial. Additionally, it reduces the insertion losses and
provides compactness.

The metasurface is a frequency agnostic concept and is finding applications
across the electromagnetic spectrum in different tasks. In the optical and
infrared frequency regions, they are of interest in connection with integrated
photonics. Although this concept [7, 8] has been initially explored in
connection with plasmonic nanoantenna [3] attention has now turned towards
all-dielectric metasurfaces [9, 10, 11, 12, 13]. Such high refractive index
dielectric nanoantenna arrays can also be considered as two dimensional high-
index contrast subwavelength diffraction gratings; various optical wavefront
manipulation possibilities have been demonstrated with these so called hcta
(high contrast transmit arrays) [14, 15, 16, 17]. In the microwave region also
low-profile transmitarrays are of interest for a variety of applications [18].
Metasurfaces can achieve phase and polarization control simultaneously and are
of interest in achieving millimeter wave beam-shaping lenses [19, 18] and
carpet-cloaks [20].

The constituent elements of a metasurface, the meta-atoms, are subwavelength
resonators while the transverse extent of the metasurface can be several
orders of magnitude larger than the operating wavelength (in other words,
useful metasurfaces will be electrical large in the transverse plane). For
heterogeneous metasurfaces, this means that the number of free parameters in
the design of the metasurface can exceed $10^{9}$ [21]. This structure is not
easily amenable for analysis, synthesis and optimization tasks.
Conventionally, the so-called unit-cell approximation [22] has been adopted
[21] in the design of heterogeneous metasurfaces, whereby each constituent
meta-atom is designed as if it was part of an infinite periodic lattice. This
unit-cell approach has also been extensively used in reflectarray design in
the microwave domain and it has been reported that it leads to discrepancies
when experiments are compared with the simulated designs [22]. The unit cell
approximation has several limitations, chiefly, it tends to be inadequate in
the presence of strong phase gradients, interactions between neighboring
pillars, and oblique incidence angles. These constraints become important for
instance in the case of a focusing lens [19, 18, 21, 23].

In this paper, we propose an approach whereby spatial order is applied on an
extended cell, which includes several individual resonator elements. This
technique is similar to the Extended Local Periodicity (ELP) approach [22]
proposed for reflectarray design. Specifically, we focus our attention to a
beam deflector element which is a commonly used structure for benchmarking
metasurfaces design strategies [23, 21, 24] as the transmission efficiencies
can be compared across designs. The beam deflector element changes the
direction of an normally incident plane wave by giving it a predesigned
inclination with respect to the normal. Although a simple element, the beam
deflector, in addition to being an important element in itself, can be
combined to produced metasurfaces with more advanced functionality like high
numerical aperture and multi-wavelength focusing lenses and holograms [21, 24,
25, 26]. The beam deflector elements considered in this paper do not use the
local phase approach [24] which is based on the unit-cell approximation but
move beyond this by considering an extended unit-cell consisting of a larger
number of nanoantenna.

The main contribution of this article is that it presents a systematic
investigation of global optimization methods for the design of the extended
unit-cells. In comparison to the beam deflector synthesis method reported by
Byrnes et. al. [21] which use local optimization methods, our method which
relies on global optimization methods like the Genetic Algorithm and
Artificial Bee Colony show significant improvement in efficiency by not
getting stuch at local optima. In comparison to the [24] which has studied the
application of Genetic Algorithm (GA) to the beam deflector design problem, we
have explored unit-cells which do not restrict the design to a rectangular
lattice and cylindrical elements. Additionally, our Artificial Bee Colony
method is seen to outperform GA method in speed by cutting the convergence
time by 35%.

The paper is organized as follows: following this introduction, in section 2,
we describe the beam deflector geometry and the optimization algorithms in
detail. In section 3, we discuss how to phrase the optimization problem in
terms of the GA and Memetic ABC algorithms and suggest how to best choose the
hyperparameters of these algorithms. The performance figures for the designs
and comparisons with previously published reports are presented in section 4
before concluding the paper in section 5.

##  2 Specification of the optimization problem

As discussed earlier, the beam deflector geometry is a canonical element in
the design of a wide variety of metasurfaces. The system considered here can
be directly compared with the geometry reported in Byrnes et. al [21] where a
wide-area focusing lens was designed using beam-deflectors as motifs.
Specifically, this required beam deflectors with deflection angles varying in
the range 20 to 70 degrees. Figure 1 shows the geometry of such a beam-
deflector with rectangular shaped extended unit-cells. Here TiO2 nanoantennae
are arranged in a hexagonal lattice within the unit-cells on a fused silica
substrate. We have kept the overall dimensions of the extended unit-cells for
any particular deflection angle exactly equal to those in Byrnes et. al [21]
so that our designs can serve as a drop-in replacements with higher
efficiencies.

Figure 1: (a) Three dimensional schematic of a beam deflector made of extended
unit-cells containing TiO2 nano elliptical antennae on fused silica substrate
(b) Design parameters of the extended unit-cell. The yellow line shows the
boundaries of the extended unit-cell.

A plane wave is normally incident from the substrate side and upon leaving the
beam deflector emerges at an angle ($\theta$ in the $x-z$ plane) that is
decided by the design (vice-versa a light wave incident at this angle from the
free-space side emerges at normal incidence). The period of the beam deflector
unit-cell in the $x$ direction is called the grating period (G. P.) and the
period in the $y$ direction is called the lateral period (L. P.). These
periods are given by:

$\displaystyle G.P.=\lambda_{0}/sin\theta,$ (1) $\displaystyle L.P.=\Phi
f\tan{\theta},$

where, $\lambda_{0}$ is free space wavelength of source light and $\theta$ is
angle of incidence or the desired deflection angle. $\Phi$ is angular width of
unit cell in $y$ direction and is decided by the lateral width and the desired
focal length. Throughout the paper we consider the lateral period as 400 nm as
in [21]. As seen in equation 2, a beam deflector with a larger deflection
angle requires a larger period.

Within each extended unit-cell, we consider elliptically shaped nanoantennae
each with freely choosable values of the following 5 parameters: major and
minor axes $r_{max}$ and $r_{min}$ respectively, the $x$ and $y$ coordinate of
the center of the ellipse in the cell coordinate system $x_{0}$ and $y_{0}$
respectively, and the orientation angle of the major axis $\alpha$. In each
unit cell the number of nanoantennae $N$ can be calculated by:

$N\leq G.P./d_{min},$ (2)

where $d_{min}$ is minimum diameter of nano structure and can be decided based
on the fabrication constraints. Here we impose the condition that the minimum
feature size of the entire shape should not be lower than 100 nm.

The optimization problem can now be stated as:

$\displaystyle\underset{s}{\text{maximize}}$ $\displaystyle
0.5(\eta^{1}(s)+\eta^{2}(s))$ (3) subject to $\displaystyle f_{i}(s)\leq
b_{i},\;i=1,\ldots,m,$

where $s$ is a $5N$ sized vector describing the extended unit cell ( i.e.
$s=[r_{max}^{i},r_{min}^{i},x_{0}^{i},y_{0}^{i},\alpha^{i}]$ and $i=1\cdots
N$), $\eta^{1}$ and $\eta^{2}$ describe the first order beam efficiency of the
designed grating for $x$ and $y$ polarized normally incident beams
respectively (the first order is made to coincide with the desired angle
$\theta$), and the various functions $f_{i}$ describe constraints imposed on
the values that the vector $s$ or its elements can take (for instance, these
constraints could require that the ratio of the major and the minor axes be
limited to 5). More generally, fabrication techniques impose several
constraints on realizable designs and these could be incorporated as shown
above. In order to calculate the first order efficiency of a particular
extended unit-cell, we use the rigorous coupled wave analysis (RCWA) approach
[27] using the open-source S4 application. S4 combines the S-matrix approach
with the RCWA method. Full-wave simulations are performed on the finally
obtained optimal vectors to verify the obtained efficiency figures.

##  3 Optimization via Genetic Algorithm and Me-ABC Algorithm

Due to the relatively large number of free parameters and stringent
fabrication related constraints required to design efficient metasurface
scanning over the full parameter space is not feasible. Stochastic
optimization methods are more appropriate to improve the efficiency and are
expected to perform better than simpler gradient descent based approaches. In
order to optimize the polarization averaged first order diffraction efficiency
we employed two techniques: Genetic algorithm and hybrid ABC algorithm with
memetic search phase.

###  3.1 Genetic Algorithm

Genetic Search is a very popular method used in electromagnetic problems
having a large multidimensional unknown search space [28, 29]. We have
employed genetic algorithm with three common parameters: selection, crossover
and mutation. Additionally, we have set upper and lower bounds for the radii
and position of nanoantennae in the beam deflector. Radius values of ellipses
are limited between 50 nm and 125 nm and the minimum distance between two
antennae centres is set to 100 nm.

Figure 2: Selection of hyperparameters for properly tuning of GA algorithm via
grid-search. (a-c) diffraction efficiency for 70 degree beam deflector as a
function of mutation and crossover probability for population sizes of 100
(a), 200 (b) and 300 (c).

The elliptical nanoantennae represent individuals in a generation. Each
generation has some number of families $N_{f}$ represented by beam deflectors
with $N_{i}$ number of individuals. The number of individuals in a family is
given by equation 2. The evolution usually starts from a population of
randomly generated individuals that constitute potential metasurface designs.
We are targeting to maximize the polarization averaged first order diffraction
efficiency. For each generation, we choose a fraction $f_{t}$ of the best
performing beam deflectors and a fraction $f_{d}$ of individuals randomly out
of remaining $1-f_{t}$ fraction of the lesser performers from the population
based on diffraction efficiency. These individuals then act as parents for
next generation so that we promote genetic diversity in addition to fitness.
The Crossover step applied next is the process of taking more than one parent
solutions and producing a child solution from them. The crossover operator
given by $P_{c}$ modifies the children generation by randomly mixing couples.
After crossover, We mutate radii, position and angle of nano-ellipses by some
amount given by mutation probability $P_{m}$. The pseudocode for the GA used
here is shown below.

Algorithm 1 Genetic Search Method

1:Generate initial population of size $N_{f}$;

2:Evaluate initial population according to fitness criteria;

3:while Termination criteria meets do

4: Select the best fit individuals for reproduction;

5: Perform crossover operation given by probability $P_{c}$;

6: Perform mutation operation given by probability $P_{m}$;

7: Evaluate fitness of new individuals;

8: Replace the worst individuals of population by best new individuals;

9:Report the best family achieved

We have used traditional Grid search method for tuning of hyper parameters of
the GA such as $f_{t},f_{d},P_{c},P_{m}$ and $N_{f}$. It was found that values
of 0.15 and 0.10 was best suited for $f_{t}$ and $f_{d}$ respectively. The
procedure followed in the grid search method is demonstrated in figure 2. It
shows how different population size, mutation rates and crossover rates affect
the convergence of the efficiency function. Population diversity is crucial to
the genetic algorithm’s ability to continue fruitful exploration without
getting stuck at local maxima. The size of the population $N_{f}$ dictates the
available diversity. While a large value will lead to good diversity it will
also increase the computational time. We observe that 300 families in each
generation is sufficient to produce good fitness. The crossover rate $P_{c}$
and mutation rate $P_{m}$ are both continuous variables within [0,1]. These
two values are important for controlling the balance between exploration and
exploitation. Table 1 summarizes a good set of values of all the
hyperparameters for the GA used in this problem.

###  3.2 ABC algorithm with Memetic Search Phase

Artificial bee colony (ABC) optimization algorithm is a relatively new
population based probabilistic approach for global optimization based on the
concept of swarm intelligence [30]. ABC has outperformed many other nature
inspired optimization algorithms prompting our interest in it for metasurface
design. Memetic algorithms (MA) is another growing area of research in
evolutionary computation. It is inspired by Darwin’s principle of evolution
and Dawkin’s notation of a ”meme” [31]. The ABC algorithm has achieved
excellent results when solving continuous and combinatorial optimization
problems [32]. To achieve the benefits of both algorithm, ABC is hybridized
with Memetic Algorithm [33, 34, 35].

In the ABC algorithm [30], the search of a parameter space is accomplished by
a set of honey bees called the swarm containing three types of bees: employed
bees, onlooker bees, and scout bees. Consider an objective function $f_{i}(s)$
which evaluates some performance metric of an extended unit-cell represented
by $s=[r_{max}^{i},r_{min}^{i},x_{0}^{i},y_{0}^{i},\alpha^{i}]$ (here $s$ is a
vector of $5N$ geometrical parameters and $N$ is number of ellipses on a beam
deflector). Each beam deflector is represented by a food source $s$ in the
swarm and it is generated as follows [33]:

$s_{ij}=s_{minj}+rand[0,1](s_{maxj}-s_{minj}),$ (4)

where function rand[0,1] returns a random value between 0 and 1, $s_{maxj}$
and $s_{minj}$ represents maximum and minimum limit of the candidate solution
$s_{i}$. The swarm will now navigate this parameter space and converge to an
optimal geometry. The parameter space is considered a region of space where
food sources are located; those regions of this parameter space that exhibit a
high value of the function $f$ are considered to be richer food sources. The
solution search equation of the original ABC algorithm is significantly
influenced by a random quantity which helps in exploration at the expense of
exploitation of the search space; there is a significant chance of skipping
true solutions due to the large step sizes that are often used. In order to
balance between diversity and convergence capability of the ABC (in other
words between exploration and exploitation), a new local search phase is
integrated with the basic ABC to exploit the search space identified by the
best individual in the swarm. The addition of a memetic search phase to the
ABC algorithm results in the the memetic ABC (MeABC) algorithm wherein the
step size required to update the best solution is controlled by a Golden
Section Search (GSS) approach. The MeABC has an enhanced exploitation
capability in comparison to the bare ABC Artificial bee colony algorithm [33,
34, 35].

In employed bee step, current solutions are changed by employed bees based on
their individual experience. If the fitness of latest solution is better than
previous solution, then the employed bees update their position to new
solution. The employed bees in ABC algorithm use the following equation in
order to improve self solution [34]:

$v_{ij}=s_{ij}+\Phi_{ij}(s_{ij}-s_{kj}),k\neq i,$ (5)

where $v_{ij}$ is the updated food, $s_{i}$ is the old food, $s_{kj}$ is a
random food from hive. Here $k\in{I,2,...,s_{N}}$ and $j\in{I,2,...,N}$ are
haphazardly chosen indices and $\Phi_{ij}$ is a random number in range
$[-1,1]$. Onlooker bees in the hive expect information of fresh solutions and
their position. In the next step onlooker bees inspect the available
information and pick a solution with a probability given by

$P_{i}=\dfrac{fit_{i}}{\sum_{j=1}^{n}fit_{i}},$ (6)

where $fit_{i}$ is $i^{th}$ solution in the swarm. If the position of a food
source is not updated for a given cycle it is considered to be abandoned. In
the scout bee step, the bee whose food source has been deserted becomes a
scout bee and the deserted food source is replaced by a haphazardly chosen
food source within the search space. Scout bees are agents for global food
search; they replace a food source by another randomly chosen food source
which is generated by the equation

$s_{ij}=s_{minj}+rand[0,1](s_{maxj}-s_{minj}),j\in{1,2\dots N},$ (7)

where $s_{minj}$ and $s_{maxj}$ are the bounds of $s_{ij}$ in the jth
direction.

Memetic phase is designed on the golden grid/section search criteria.
Basically we choose a negative value of a (that is some value on the left of
the X-axis for that dimension of the best food) and equal but positive value
of b (that is some value on the right of the X-axis) [31]. The values a and b
dictate how wide the memetic search should be. The current best food lies
exactly in the middle initially. The memetic search starts by determining
whether a better food is on the Left of the current best food (i.e. more
towards a) or to the right (i.e. more towards b). Similar to a binary search
algorithm it updates the values of a and b and gradually zeroes in on the
better solution around the current best. If it finds a better food source, it
updates the best food else the search fails. This process is carried out until
the difference between a and b (absolute value of the difference) is greater
than the set value of epsilon; epsilon being the stopping criteria of memetic
search phase. In MeABC algorithm, ABC algorithm behaves as a local search
algorithm in which only the best individual of the current swarm updates
itself in its neighbourhood while in memetic search phase the step size
required to update the best individual in the current swarm is controlled by
the golden section search (GSS) approach. GSS processes the interval [a=-0.75,
b=0.75] and generates two intermediate points:

$\displaystyle F_{1}=b-(b-a)\times\Psi,$ (8) $\displaystyle
F_{2}=a+(b-a)\times\Psi,$

where $\Psi$ is the golden ratio. Memetic ABC algorithm has three steps
similar to the ABC algorithm and one more step, the memetic phase, is added
for updating the location of an individual. It changes position [34] given by
the equation

$s_{ij}=s_{i}+\Phi_{i}(s_{i}-s_{k})+\Psi(s_{best}-s_{i}),$ (9)

where $\Phi_{i}$ is a random number in the interval $[0,D]$ and $D$ is a
positive constant. The pseudo code for the Memetic ABC algorithm that we have
employed is given below.

Algorithm 2 MeABC Algorithm

1:Generate an initial population of food sources using equation 2;

2:Evaluate initial population according to fitness criteria;

3:while Termination criteria meets do

4: Deploy employed bee searches to find new food searches in the neighbourhood
using equation 5;

5: Calculate Probability P for each food source using equation 6;

6: Send onlooker bee to food source depending upon P;

7: Evaluate fitness of each new food source;

8: if any employed bee becomes scout bee then

9: Send scout bee to a randomly generated food source;

10: Employ memetic search phase;

11: Memorize the best food source achieved so far;

12:Report the best food source achieved

From equation 4, it can be observed that the step size consists of a random
component $\Phi$ and thus a proper balance is not possible manually [34]
($\Phi$ is random component that decides direction and step size of an
individual). Memetic search phase (MSP) improves the exploitation capability
[36] considerably. We have used MSP to fine tune the value of $\Phi$
dynamically and iteratively using the Golden Section Search strategy. The
range of $\Phi_{ij}$ is set to $[a,b]$ where $a=-0.75$ and $b=0.75$; $a$ and
$b$ dictate how wide the memetic search should be. To tune the hyper
parameters of MeABC algorithm like the swarm population and numbers of each
kind of bees, we used a traditional grid search method. The hyperparameter
values that were found to yield good results in terms of efficiency are
summarized in Table 1.

Table 1: Hyperparameters chosen for the GA and MeABC method  GA Parameters | Me-ABC Parameter
---|---
Parameters |  Value |  Parameters |  Value
No. of Families |  300 |  No. of Food Sources |  25
No. of Individuals in a Family |  2-4 |  No. of Onlooker Bees |  25
Mutation Probability, $P_{r}$ |  0.07 |  No. of Employed Bees |  25
Crossover Probability, $P_{c}$ |  0.90 |  No. of Scout Bees |  1
No. of Generations |  50 |  No. of iterations before bee is tired |  200

##  4 Results and discussion

The above design methodology was applied to the design of a set of beam
deflecting metasurfaces operating at 580 nm wavelength. The extended unit
cells are rectangular shaped but within it the fill fraction is maximized by
adopting a hexagonal grid. The elliptically shaped nanoantenna are made of
$TiO_{2}$ and sit atop a fused silica substrate. The refractive index of
$n(TiO_{2})=2.37$ at $580nm$ and the nanoantenna height is kept at $h=550nm$.

Figure 3: (a) Comparison of the polarization-averaged first order diffraction
efficiencies for beam deflectors with bend angles ranging from 20 to 70
degrees designed via Gradient Ascent, Genetic and MeABC Approach. Top view of
the optimal beam deflectors at deflection angles of 30 degree (b) and 50
degrees (c) (red lines delineate the extended unit-cell). Full-wave simulation
based far field intensity plots (d) and (e) for the gratings (b) and (c)
respectively.  Table 2: Optimal geometrical parameters of the extended unit-
cells and their first order polarization-averaged diffraction efficiencies [%]
at $580nm$ wavelength

| E1 | E2 | E3 | E4 |
---|---|---|---|---|---
D.A. | Rx | Ry | $\theta_{r}$ | Rx | Ry | $\theta_{r}$ | Rx | Ry | $\theta_{r}$ | Rx | Ry | $\theta_{r}$ | Eff.
20∘ | 69.15 | 83.17 | -9.6 | 50 | 50 | 0 | 50 | 50 | 0 | 162.2 | 89.14 | 0.01 | 83.8%
30∘ | 111.14 | 85.63 | 175.7 | 50.04 | 50.06 | 175.16 | 95.99 | 134.80 | 88.66 | - | - | - | 79%
40∘ | 93.25 | 101.76 | 78.36 | 50 | 50 | 0 | - | - | - | - | - | - | 78.6%
50∘ | 58.09 | 77.46 | 179.99 | 50 | 50 | 0 | - | - | - | - | - | - | 81%
60∘ | 106.21 | 79.77 | -179.96 | 50 | 50 | 0 | - | - | - | - | - | - | 71.8%
70∘ | 103 | 89.86 | -177.76 | 50 | 50 | 0 | - | - | - | - | - | - | 65%

Figure 3 (a) shows comparison of the polarization averaged efficiencies
obtainable with all three methods viz: the local steepest gradient ascent
[21], the genetic algorithm, and ABC with memetic phase search. Firstly, note
that all three curves indeed follow the well known fact that efficiencies are
nearly constant up to deflection angles of about 50 degrees but start to decay
rapidly afterwards. While all three methods give nearly equal efficiencies in
the initial angular range, the designs obtained with the global optimization
methods significantly outperform at steeper angles.

Figure 4: Comparison of convergence time needed to arrive at the optimal beam
deflectors for gradient ascent, GA and MeABC methods. Three different
deflection angles have been considered.

Comparison of the convergence times using all three methods for three
different beam deflectors is shown in figure 4. All runs were made on machine
with the following configurations: Processor - Intel (R) Xeon(R) CPU E5-2650
<EMAIL_ADDRESS>RAM - 32 GB; and, System Type - Windows 64 bit operating system.
It is evident that the gradient ascent method gets stuck early on in a local
optima and never improves the overall efficiency. The global optimization
techniques clearly avoid this problem. Furthermore, the time taken by the
global search is not too large in comparison to the local methods. While the
efficiencies achieved by GA and MeABC are almost similar, the MeABC approach
converges to the final geometry in approximately $35\%$ time compared to the
GA approach.

##  5 Conclusion

We have presented two global optimization techniques based on nature inspired
algorithms for the rapid design of all-dielectric metasurfaces and have
applied them to the case of beam deflectors. Compared to local optimization
methods like the gradient ascent algorithm, Genetic and Me-ABC algorithms
provides larger deflection efficiencies. Up to $15\%$ efficiency improvement
is achieved for higher deflection angles. The MeABC method proposed by us is
significantly faster than previously proposed gradient ascent algorithm [21].
It also outperformed the GA in terms of computation time.

The authors acknowledge support from the Department of Science and Technology,
Govt. of India through the Extramural grant SB/S3/EECE/0200/2015.

## References

  * [1] Christopher L Holloway, Edward F Kuester, Joshua a. Gordon, John O’Hara, Jim Booth, and David R Smith.  An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials.  IEEE Antennas and Propagation Magazine, 54(2):10–35, 2012.
  * [2] Patrice Genevet and Federico Capasso.  Holographic optical metasurfaces: a review of current progress.  Reports on progress in physics. Physical Society (Great Britain), 78(2):24401, 2015.
  * [3] Nikolay I Zheludev and Yuri S Kivshar.  From metamaterials to metadevices.  Nature Materials, 11(11):917–924, 2012.
  * [4] D Lin, P Fan, E Hasman, and M L Brongersma.  Dielectric gradient metasurface optical elements.  Science, 345(6194):298–302, 2014.
  * [5] Nina Meinzer, William L Barnes, and Ian R Hooper.  Plasmonic meta-atoms and metasurfaces.  Nature Photonics, 8(12):889–898, 2014.
  * [6] Alexander E Minovich, Andrey E Miroshnichenko, Anton Y Bykov, Tatiana V Murzina, Dragomir N Neshev, and Yuri S Kivshar.  Functional and nonlinear optical metasurfaces.  Laser and Photonics Reviews, 9(2):195–213, 2015.
  * [7] Nanfang Yu, Patrice Genevet, Mikhail a Kats, Francesco Aieta, Jean-Philippe Tetienne, Federico Capasso, and Zeno Gaburro.  Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction.  Science, 334(October):333–337, 2011.
  * [8] Nanfang Yu and Federico Capasso.  Flat optics with designer metasurfaces.  Nature Materials, 13(2):139–150, 2014.
  * [9] Jacob B Khurgin.  How to deal with the loss in plasmonics and metamaterials.  Nature Nanotechnology, 10(1):2–6, 2015.
  * [10] Francesco Aieta, Patrice Genevet, Mikhail a. Kats, Nanfang Yu, Romain Blanchard, Zeno Gaburro, and Federico Capasso.  Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces.  Nano Letters, 12(9):4932–4936, 2012.
  * [11] Salvatore Campione, Lorena I Basilio, Larry K Warne, and Michael B Sinclair.  Tailoring dielectric resonator geometries for directional scattering and Huygens’ metasurfaces.  Optics Express, 23(3):2293, 2015.
  * [12] Manuel Decker, Isabelle Staude, Matthias Falkner, Jason Dominguez, Dragomir N Neshev, Igal Brener, Thomas Pertsch, and Yuri S Kivshar.  High-Efficiency Dielectric Huygens’ Surfaces.  Advanced Optical Materials, 3(6):813–820, 2015.
  * [13] Isabelle Staude, Andrey E Miroshnichenko, Manuel Decker, Nche T Fofang, Sheng Liu, Edward Gonzales, Jason Dominguez, Ting Shan Luk, Dragomir N Neshev, Igal Brener, and Yuri Kivshar.  Tailoring directional scattering through magnetic and electric resonances in subwavelength silicon nanodisks.  ACS Nano, 7(9):7824–7832, 2013.
  * [14] a Arbabi, Y Horie, M Bagheri, and a Faraon.  Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission.  Nat Nanotechnol, 10(11):937–943, 2015.
  * [15] Amir Arbabi, Yu Horie, Alexander J Ball, Mahmood Bagheri, and Andrei Faraon.  Subwavelength-thick Lenses with High Numerical Apertures and Large Efficiency Based on High Contrast Transmitarrays.  Nature communications, 6(May):1–10, 2015.
  * [16] Ehsan Arbabi, Amir Arbabi, Seyedeh Mahsa Kamali, Yu Horie, and Andrei Faraon.  Multiwavelength polarization insensitive lenses based on dielectric metasurfaces with meta-molecules.  Optica, 0(6):1–12, 2016.
  * [17] Amir Arbabi, Yu Horie, Alexander J Ball, Mahmood Bagheri, and Andrei Faraon.  Efficient high NA flat micro-lenses realized using high contrast transmitarrays.  In Spie Opto, volume 9372, page 93720P, 2015.
  * [18] Carl Pfeiffer, Student Member, and Anthony Grbic.  Millimeter-Wave Transmitarrays for Wavefront and Polarization Control.  61(12):4407–4417, 2013.
  * [19] Carl Pfeiffer and Anthony Grbic.  Planar Lens Antennas of Subwavelength Thickness : Collimating Leaky-Waves With Metasurfaces.  63(7):3248–3253, 2015.
  * [20] Li Yi Hsu, Thomas Lepetit, and Boubacar Kant.  Extremely Thin Dielectric Metasurface for Carpet Cloaking.  Progress In Electromagnetics Research, 152(March):33–40, 2015.
  * [21] Steven J. Byrnes, Alan Lenef, Francesco Aieta, and Federico Capasso.  Designing large, high-efficiency, high-numerical-aperture, transmissive meta-lenses for visible light.  Optics Express, 24(5):5110, 2016.
  * [22] Min Zhou, Stig Busk Sørensen, Erik Jørgensen, Peter Meincke, Oleksiy S Kim, and Olav Breinbjerg.  Analysis of Printed Reflectarrays Using Extended Local Periodicity.  pages 1494–1498.
  * [23] Ye Feng Yu, Alexander Y Zhu, Ramón Paniagua-Dom$\backslash$’$\backslash$inguez, Yuan Hsing Fu, Boris Luk’yanchuk, and Arseniy I Kuznetsov.  High-transmission dielectric metasurface with 2$\pi$ phase control at visible wavelengths.  Laser and Photonics Reviews, 9(4):412–418, 2015.
  * [24] Ramat Aviv.  Genetically optimized all-dielectric metasurfaces.  25(3):937–943, 2017.
  * [25] J Cheng, D Ansari-Oghol-Beig, and H Mosallaei.  Wave manipulation with designer dielectric metasurfaces.  Opt Lett, 39(21):6285–6288, 2014.
  * [26] Zhongyi Guo, Lie Zhu, Fei Shen, Hongping Zhou, and Rongke Gao.  Dielectric metasurface based high-efficiency polarization splitters.  RSC Adv., 7(16):9872–9879, 2017.
  * [27] Victor Liu and Shanhui Fan.  S 4: A free electromagnetic solver for layered periodic structures.  Computer Physics Communications, 183(10):2233–2244, 2012.
  * [28] J. Michael Johnson and Yahya Rahmat-Samii.  Genetic algorithms in engineering electromagnetics.  IEEE Antennas and Propagation Magazine, 39(4):7–21, 1997.
  * [29] Guillaume Tremblay, Jean Numa Gillet, Yunlong Sheng, Martin Bernier, and Gilles Paul-Hus.  Optimizing fiber Bragg gratings using a genetic algorithm with fabrication-constraint encoding.  Journal of Lightwave Technology, 23(12):4382–4386, 2005.
  * [30] Dervis Karaboga and Bahriye Basturk.  A powerful and efficient algorithm for numerical function optimization: Artificial bee colony (ABC) algorithm.  Journal of Global Optimization, 39(3):459–471, 2007.
  * [31] P Moscato and Carlos Cotta.  A Gentle Introduction to Memetic Algorithms.  Building, 7(19):1–36, 2003.
  * [32] Asaju Laaro Bolaji, Ahamad Tajudin Khader, Mohammed Azmi Al-Betar, and Mohammed A. Awadallah.  Artificial bee colony algorithm, its variants and applications: A survey.  Journal of Theoretical and Applied Information Technology, 47(2):434–459, 2013.
  * [33] Iztok Fister, Janez Brest, and Viljem Zumer.  Memetic Artificial Bee Colony Algorithm for Large-Scale Global Optimization §.  Proc. IEEE Congress on Evolutionary Computation, pages 1–21, 2012\.
  * [34] Jagdish Chand Bansal, Harish Sharma, K V Arya, and Atulya Nagar.  Memetic search in artificial bee colony algorithmBansal, J. C., Sharma, H., Arya, K. V, & Nagar, A. (2013). Memetic search in artificial bee colony algorithm. Soft Computing, 17(10), 1911–1928. https://doi.org/10.1007/s00500-013-1032-8.  Soft Computing, 17(10):1911–1928, 2013.
  * [35] Sandeep Kumar, Vivek Kumar Sharma, and Rajani Kumari.  Memetic Search in Artificial Bee Colony Algorithm with Fitness based Position Update.  IEEE International Conference on Recent Advances and Innovations in Engineering (ICRAIE-2014),, 2014.
  * [36] Hongfeng Wang, Dingwei Wang, and Shengxiang Yang.  A memetic algorithm with adaptive hill climbing strategy for dynamic optimization problems.  Soft Computing, 13(8-9):763–780, 2009.


#  An illustration of canonical
quantum-classical dynamics
Backreaction, canonical relations and time evolution in the quantum-classical
harmonic oscillator

Mustafa Amin<EMAIL_ADDRESS>Mark A. Walton<EMAIL_ADDRESS>
###### Abstract

Using the example of the harmonic oscillator, we illustrate the use of hybrid
dynamical brackets in analyzing quantum-classical interaction. We only assume
that a hybrid dynamical bracket exists, is bilinear, and reduces to the pure
quantum/classical bracket when acting on pure quantum/classical variables. Any
hybrid bracket obeying these natural requirements will produce the same
dynamics for pure classical or quantum variables, given a hybrid Hamiltonian.
Backreaction is manifested in the evolution of a nonvanishing commutator
between classical variables. The more massive the classical system is, the
less it is affected by backreaction. Interestingly, we show that while pure
variables evolve to violate the pure canonical relations, they always obey the
hybrid canonical relations. The dynamics of hybrid variables, on the other
hand, is shown to require a fully specified and consistent hybrid bracket,
otherwise evolution cannot be defined.

##  Introduction

Quantum-classical systems appear in numerous fields, from chemical and
condensed matter physics to cosmology and black hole physics. In practical
applications, a full quantum treatment is often computationally intractable.
Treating parts of the system as classical simplifies calculations and provide
easier models to conceptualize and approach such problems. We constantly deal
with quantum electrons around classical nuclei, quantum matter on classical
spacetimes or quantum systems interacting with external classical fields. Less
practically, foundational research on the nature of gravity or the measurement
problem in quantum mechanics may also benefit from the study of quantum-
classical systems.

Familiar approaches to modelling quantum-classical interactions usually fall
into two categories. The first treats the classical part of the system as
omnipotent; it acts on, but is not acted back on by, the quantum part.
Examples of this type include textbook quantum mechanics [1], or quantum field
theories on curved spacetime.

The second approach is to use averages of quantum observables $\langle
x_{Q}\rangle$ in classical equations

$\text{Eq.}_{\text{Classical}}\big{(}\langle x_{Q}\rangle\,,y_{C}\big{)}~{},$

and then feed the classical observables $y_{C}$ back into quantum equations

$\text{Eq.}_{\text{Quantum}}\big{(}x_{Q}\,,y_{C}\big{)}~{}.$

While this approach captures a form of quantum backreaction on the classical
system, it suffers from two problems: 1. it averages out quantum fluctuations
so they do not manifest on the classical side, and 2. it introduces
nonlinearity to the quantum Schrödinger equation thus negating the
superposition principle. Even if one ignores the loss of superposition as a
side effect of a simplifying approximation, nonlinear equations can be very
difficult to handle in practical applications.

In this paper, we are concerned with a third approach: directly coupling
quantum and classical variables. This requires defining a consistent framework
for the dynamics of such hybrid quantum-classical variables. Namely, we are
interested in a canonical dynamical structure where the equations of motion
are defined in terms of a hybrid dynamical bracket and canonical relations
between conjugate variables. Many proposals for a hybrid quantum-classical
bracket can be found in the literature [2, 3, 4, 5, 6, 7, 8, 9]. See also [10,
11, 12].

However, impediments to a consistent hybrid dynamics have been found [13, 14,
15, 16, 17]. Specifically, it was shown that hybrid dynamical brackets do not
satisfy the Jacobi identity or the Leibniz rule. These are crucial properties
for dynamical consistency. The Jacobi identity guarantees that the fundamental
canonical relations are preserved with dynamical evolution. The Leibniz rule
is essential for defining dynamical evolution and derivatives, as shall be
demonstrated in Sec. 4.

The present authors, however, have shown in [18] that starting from a full
quantum system and taking the classical limit of a part, one can uncover a
general class of hybrid brackets. Different brackets arise from different
quantization schemes on the classical sector prior to taking the classical
limit. In particular, the bracket proposed by Aleksandrov [2], Gerasimenko
[3], and Boucher and Traschen [4] arises as a special case of the general
bracket when Wigner-Weyl quantization is used. The derivation of the general
bracket makes use of the phase space formulation of quantum mechanics showing
a connection between phase space distributions (like that of Wigner, Husimi,
etc.) and hybrid dynamical brackets.

To address the aforementioned no-go theorems, the authors introduced a hybrid
composition product [18]. The general hybrid bracket is then the commutator of
that product. It follows that the Jacobi identity and Leibniz rule are
automatically satisfied provided that the hybrid composition product is
associative.

Then, for a restricted set of hybrid variables that form an associative
subalgebra with the hybrid composition product, dynamics is consistent. Since
quantum-classical interactions are represented by hybrid terms in the
Hamiltonian, the consistency of quantum-classical dynamics in our scenario
implies that only certain interactions between quantum and classical systems
are allowed.

Here we present a working example of canonical hybrid dynamics that provides a
blueprint for applying the framework to possible problems of interest. The
harmonic oscillator is ubiquitous in physical systems. We use it to illustrate
the features of hybrid dynamics and the cautionary pitfalls of dealing with
systems that are not clearly classical or quantum.

Notably, hybrid evolution of pure variables does not depend on any particular
definition of a hybrid bracket as long as it is bilinear and reduces to a pure
bracket when one of its arguments is pure. Evolution of hybrid variables, on
the other hand, requires specification of the hybrid bracket used. We
emphasize that the methods presented in this paper are applicable to more
general problems.

We introduce basic concepts and general properties of hybrid brackets in Sec.
2. In Sec. 3, we calculate the equations of motion of pure classical and
quantum variables in the quantum-classical harmonic oscillator. Adopting
Anderson’s argument [5], we find the “backreaction” of the quantum on the
classical and show the blurring of the line between quantum and classical
variables due to interaction.

Finally, we describe an example of time evolution of hybrid variables in Sec.
4. In that example, the need for a consistent hybrid bracket is illustrated.
Different brackets, corresponding to different quantization schemes (and, by
extension, different phase space distributions) are shown to produce different
time evolutions for hybrid variables.

##  Hybrid dynamics

There are multiple approaches to the problem of combining quantum and
classical systems [10, 11, 12]. In this paper, we are concerned with the
canonical approach, in which the problem is that of defining a hybrid
dynamical bracket. The hybrid bracket $\\{\\![\cdot\,,\cdot]\\!\\}$ then
plays, for hybrid variables, the role of the commutator
$[\cdot\,,\cdot]/i\hbar$ for quantum variables or the Poisson bracket
$\\{\cdot\,,\cdot\\}$ for classical ones.

Hybrid brackets have been proposed in [2, 3, 4, 5, 6, 7, 8, 9]. The authors
have derived a general class of hybrid brackets in [18] by applying a partial
classical limit to quantum mechanics. Different arguments may produce
different brackets, but a general feature of all of them is their reduction
property. A hybrid bracket should reduce to a pure (quantum or classical)
bracket when one of its arguments is pure:

$\displaystyle\begin{split}\\{\\![\cdot\,,\,A_{Q}]\\!\\}&=\frac{1}{i\hbar}[\cdot\,,\,A_{Q}]~{},\\\
\\{\\![\cdot\,,\,A_{C}]\\!\\}&=\\{\cdot\,,\,A_{C}\\}~{}.\end{split}$ (2.1)

The subscripts $Q$ and $C$ signify quantum and classical variables,
respectively. Another property shared by hybrid brackets is bilinearity:

$\displaystyle\begin{split}\\{\\![A,B+C]\\!\\}=\\{\\![A,B]\\!\\}+\\{\\![A,C]\\!\\}~{},\\\
\\{\\![A+B,C]\\!\\}=\\{\\![A,C]\\!\\}+\\{\\![B,C]\\!\\}~{}.\end{split}$ (2.2)

A general Hamiltonian is given by

$\displaystyle H=H_{Q}+H_{C}+V~{}.$ (2.3)

While the pure $H_{Q}$ and $H_{C}$ encode internal dynamics of the quantum and
classical subsystems, $V$ is a hybrid interaction term that couples quantum
and classical variables. Without knowing the specific definition of the hybrid
bracket, one can find the equation of motion for pure variables using only the
reduction requirement (2.1)

$\displaystyle\frac{d}{dt}A_{\text{Pure}}=\\{\\![A_{\text{Pure}},H]\\!\\}+\frac{\partial}{\partial
t}A_{\text{Pure}}~{}.$ (2.4)

The equation is given in the Heisenberg picture, where time evolution resides
in the dynamical variables instead of the state.

The reduction requirement (2.1) can go a long way in defining hybrid dynamics
as shown by an explicit example in Sec. 3. Specifically, when one is
interested in calculating the time evolution of pure variables, as in Eq.
(2.4), one can rely on the well-defined quantum or classical brackets to
perform calculations.

The dynamics of hybrid variables is more complicated, however. It was shown in
[13, 14, 15, 16, 17] that hybrid brackets cannot define a consistent dynamical
framework for general hybrid variables.

To deal with the problem of evolving hybrid variables or, equivalently, having
both arguments of the dynamical bracket be hybrid, we need a specific
definition of the bracket. The authors have derived a general hybrid bracket
in [18] along with a definite consistency condition for hybrid brackets.
Namely, only a certain subset of hybrid variables is allowed into the theory
for which hybrid dynamics is consistently defined. An immediate consequence of
such restriction is that the consistency of hybrid dynamics dictates what kind
of interaction could exist between quantum and classical systems. The
derivation of the general hybrid bracket in [18], as outlined in Sec. 4,
relies on the methods of phase space quantum mechanics.

##  Backreaction in the harmonic oscillator

One of the main motivations for developing hybrid dynamics is to study the
effect of quantum systems on classical ones. It is to be expected that
initially deterministic classical variables will inherit some uncertainty by
interacting with a quantum system. In [5], Anderson argues that classical
variables will evolve to become “quasiclassical”: variables that exhibit only
secondary fluctuations. That is, the uncertainty in a quasiclassical variable
exists only due to interacting with a quantum system; in the absence of such
interaction, quasiclassical variables evolve in purely classical fashion. This
is in contrast to fully quantum variables, which exhibit primary fluctuations;
those are fundamentally uncertain.

In this section, following Anderson’s pioneering work, we provide an example
of the evolution of (quasi)classical variables. Quantum backreaction on
classical variables is found in terms of a nonvanishing commutator between
classical canonical conjugates. We emphasize that here we only use the
reduction property (Eq. (2.1)) and not a specific definition of the hybrid
bracket.

It should be stressed that the following example only serves as an explicit
illustration of the more general possibilities of hybrid dynamics. The same
logic applies to any Hamiltonian (2.3), even if more involved equations of
motion result. The example presented here demonstrates the basic concepts.

Consider a quantum-classical harmonic oscillator Hamiltonian

$\displaystyle
H=\frac{p_{C}^{2}}{2m_{C}}+\frac{p_{Q}^{2}}{2m_{Q}}+\frac{1}{2}k(x_{C}-x_{Q})^{2}~{}.$
(3.1)

The mass, position and momentum of the classical and quantum particles are
given by $(m_{C},x_{C},p_{C})$ and $(m_{Q},x_{Q},p_{Q})$ respectively, and $k$
is the coupling strength. Using the center-of-mass and relative separation
coordinates

$\displaystyle\begin{split}X&=\frac{m_{C}\,x_{C}+m_{Q}\,x_{Q}}{m_{C}+m_{Q}}~{},\quad
P=p_{C}+p_{Q}\\\ x&=x_{C}-x_{Q}~{},\qquad\qquad
p=\frac{m_{Q}\,p_{C}-m_{C}\,p_{Q}}{m_{C}+m_{Q}}~{},\end{split}$ (3.2)

and the total and reduced masses

$\displaystyle M=m_{C}+m_{Q}~{},\qquad\quad
m=\frac{m_{C}\,m_{Q}}{m_{C}+m_{Q}}~{},$ (3.3)

the Hamiltonian can be written as

$\displaystyle H=\frac{P^{2}}{2M}+\frac{p^{2}}{2m}+\frac{1}{2}kx^{2}~{}.$
(3.4)

Notice that $(X,P,x,p)$ are hybrid variables. However, they are additive
hybrids that do not directly couple quantum and classical variables.

Formally, the new variables obey the canonical relations in terms of the
hybrid bracket. This can be shown using only the general reduction property
(2.1) and bilinearity (2.2) :

$\displaystyle\\{\\![X,P]\\!\\}$
$\displaystyle=\frac{1}{M}\\{\\![m_{C}x_{C}+m_{Q}x_{Q},p_{C}+p_{Q}]\\!\\}$
(3.5)
$\displaystyle=\frac{m_{C}}{M}\\{x_{C},p_{C}\\}+\frac{1}{i\hbar}\frac{m_{Q}}{M}[x_{Q},p_{Q}]=1~{},$
(3.6)

and similarly for the rest of the relations. Using the hybrid fundamental
relations, we can show that $\dot{X}=\\{\\![X,H]\\!\\}$ by explicit
calculation

$\displaystyle\dot{X}$
$\displaystyle=\frac{m_{C}}{M}\dot{x}_{C}+\frac{m_{Q}}{M}\dot{x}_{Q}$ (3.7)
$\displaystyle=\frac{m_{C}}{M}\\{x_{C},H\\}+\frac{1}{i\hbar}\frac{m_{Q}}{M}[x_{Q},H]$
(3.8) $\displaystyle=\\{\\![X,H]\\!\\}~{}.$ (3.9)

The reduction and bilinearity properties were again used to obtain the last
equality.

Now the system can be solved using standard methods. The center-of-mass
position and momentum evolve according to

$\displaystyle P=P(0)~{},\qquad X=X(0)+\frac{P(0)}{M}t~{}.$ (3.10)

The relative separation position and momentum equations of motion are

$\displaystyle x$
$\displaystyle=x(0)\cos(\omega\,t)+\frac{p(0)}{m\omega}\sin(\omega\,t)~{},$
(3.11) $\displaystyle p$
$\displaystyle=p(0)\cos(\omega\,t)-m\omega\,x(0)\sin(\omega\,t)~{},$ (3.12)

where $\omega=\sqrt{k/m}$ . Finally, the equations of motion for the classical
position and momentum are found to be

$\displaystyle\begin{split}x_{C}=&~{}\frac{1}{M}\big{[}m_{C}+m_{Q}\cos(\omega\,t)\big{]}x_{C}(0)\\\
&+\frac{1}{m_{C}M\omega}\big{[}m_{C}\omega\,t+m_{Q}\sin(\omega\,t)\big{]}p_{C}(0)\\\
&+\frac{m_{Q}}{M}\big{[}1-\cos(\omega\,t)\big{]}x_{Q}(0)\\\
&+\frac{1}{M\omega}\big{[}\omega\,t-\sin(\omega\,t)\big{]}p_{Q}(0)\end{split}$
(3.13)

and

$\displaystyle\begin{split}p_{C}=&~{}\frac{1}{M}\big{[}m_{C}+m_{Q}\cos(\omega\,t)\big{]}p_{C}(0)-m\omega\sin(\omega\,t)x_{C}(0)\\\
&+\frac{m_{C}}{M}\big{[}1-\cos(\omega\,t)\big{]}p_{Q}(0)+m\omega\sin(\omega\,t)x_{Q}(0)~{}.\end{split}$
(3.14)

The equations of motion for $x_{Q}$ and $p_{Q}$ are identical and can be found
by swapping subscripts $Q$ and $C$ in (3.13) and (3.14).

We see from the dependence of $x_{C}$ and $p_{C}$ on $x_{Q}(0)$ and $p_{Q}(0)$
that a form of backreaction is present on the classical variables. While the
commutator $[x_{C},p_{C}]$ vanishes initially, it evolves to be non-zero due
to interaction:

$\displaystyle\frac{1}{i\hbar}[x_{C},p_{C}]=\frac{m}{M}\big{(}2-2\cos(\omega\,t)-\omega\,t\sin(\omega\,t)\big{)}~{}.$
(3.15)

If a nonvanishing commutator indicates uncertainty, then the above equation
shows how quantum backreaction is manifested in a quantum-classical harmonic
oscillator. Fig. 1 depicts the evolution of the commutator of $x_{C}$ and
$p_{C}$. Initially, $x_{C}(0)$ and $p_{C}(0)$ were known with complete
certainty, perhaps through measurement. Prediction or retrodiction at times
other than $t=0$ is uncertain. We see that as the classical mass $m_{C}$
becomes large compared to the quantum mass $m_{Q}$, the commutator
$[x_{C},p_{C}]$ approaches zero and classical behaviour dominates $x_{C}$ and
$p_{C}$.

Figure 1: Qualitative evolution of $[x_{C},p_{C}]/i\hbar$ through Eq. (3.15),
with $m_{Q}$ and $k$ set to unity. At $t=0$, classical variables are known
with certainty (vanishing commutator). Before and after this moment of
certainty, the commutator deviates from zero implying uncertain
prediction/retrodiction. As the classical mass $m_{C}$ becomes larger, the
curve flattens approaching zero around the point of certainty at $t=0$.

Further, non-classical behaviour is shown in the change of the value of the
Poisson bracket. While $\\{x_{C},p_{C}\\}$ initially obeys the classical
canonical relation, it evolves to deviate from unity

$\displaystyle\begin{split}\\{x_{C},p_{C}\\}&=\frac{1}{M^{2}}\big{(}m_{C}+m_{Q}\cos(\omega\,t)\big{)}^{2}\\\
&\quad+\frac{m}{m_{C}M}\sin(\omega\,t)\big{(}m_{C}\omega\,t+m_{Q}\sin(\omega\,t)\big{)}~{}.\end{split}$
(3.16)

That deviation from classicality is, again, due to quantum backreaction.

The hybrid nature of $x_{C}$ and $p_{C}$ is now shown elegantly in terms of
the hybrid bracket. While they do not obey either the quantum or the classical
canonical relations (Eqs. (3.15) and (3.16)), they do obey the hybrid
fundamental canonical relation

$\displaystyle\\{\\![x_{C},p_{C}]\\!\\}=1~{}.$ (3.17)

Since $x_{C}$ and $p_{C}$ at $t\neq 0$ are additive hybrids, this can be shown
explicitly using the reduction (2.1) and bilinearity (2.2) properties applied
to Eqs. (3.13) and (3.14) for $x_{C}$ and $p_{C}$. Any hybrid bracket
possessing these properties will produce the same results. A similar
calculation can be done for $x_{Q}$ and $p_{Q}$.

An important note is in order. Thanks to the reduction property (2.1), hybrid
evolution of pure variables is consistent. That is, if the hybrid bracket has
one of its two arguments pure (or a sum of pure variables), it reduces to a
consistent pure bracket. This will be the case if only one hybrid variable is
used in the problem. This covers a wide class of problems. Usually the one
interesting hybrid variable is the interaction Hamiltonian, while all other
variables of interest are either pure classical or pure quantum such as
positions, momenta, fields etc. In the next section we will demonstrate the
use of the full hybrid bracket if one desires to study the evolution of
general hybrid variables.

##  Evolution of hybrid variables

As demonstrated in the previous section, the evolution of pure variables and,
by extension, additive hybrid variables is largely bracket-agnostic. For
nontrivial hybrids, however, a specific definition of the hybrid bracket is
necessary. Consider, for example, a hybrid variable $\eta$ that is a product
of quantum and classical variables

$\displaystyle\eta=\eta_{C}\eta_{Q}~{}.$ (4.1)

On one hand, taking the time derivative of such a variable is unclear, since
the Leibniz rule for derivation now runs into operator ordering ambiguities.
As seen in the previous section, classical variables will evolve into
quasiclassical variables that do not commute with quantum ones. Thus, while
$\eta=\eta_{C}\eta_{Q}=\eta_{Q}\eta_{C}$, one does not know whether
$\dot{\eta}=\dot{\eta}_{C}\eta_{Q}+\eta_{C}\dot{\eta}_{Q}$ or
$\dot{\eta}=\eta_{Q}\dot{\eta}_{C}+\eta_{C}\dot{\eta}_{Q}$, for example. On
the other hand, the reduction property (2.1) is clearly insufficient to
calculate the hybrid bracket when both variables are nontrivial hybrids so
that $\\{\\![\eta,H]\\!\\}$ isn’t defined either. A consistent Leibniz rule
for the hybrid bracket is needed.

The problem intensifies as various no-go theorems assert that the Leibniz rule
and the Jacobi identity are not generally satisfied by hybrid brackets, making
them inconsistent [13, 14, 15, 16, 17]. However, the authors have proposed a
possible reinterpretation of these no-go theorems in light of the work
presented in [18]. A general hybrid bracket is derived there, of the form

$\displaystyle\\{\\![A,B]\\!\\}=\frac{1}{i\hbar}\big{(}A\circledast
B-B\circledast A\big{)}~{},$ (4.2)

where $\circledast$ is the hybrid composition product. Now the consistency of
the bracket is seen in light of this new composition product. If $\circledast$
is associative, the bracket will obey the Jacobi identity, and the Leibniz
rule will take the form

$\displaystyle\\{\\![A,B\circledast C]\\!\\}=\\{\\![A,B]\\!\\}\circledast
C+B\circledast\\{\\![A,C]\\!\\}~{}.$ (4.3)

Now that the problem is cast in definite terms, we propose a possible
circumvention. If we restrict hybrid variables to be only those forming an
associative subalgebra with the hybrid composition product $\circledast$, then
the dynamics of these variables is consistent. This implies that only certain
quantum-classical interactions are allowed, as dictated by the consistency of
the framework.

Since associativity with the $\circledast$-product is the central condition, a
concrete definition of $\circledast$ is needed. In [18], it is shown that a
general hybrid bracket is derived through the application of a partial
classical limit to a full quantum theory. A quantum system is subdivided into
$Q$(uantum) and $C$(lassical) sectors. Then a classical limit is applied on
the $C$ sector using the phase space formulation of quantum mechanics. The
form of $\circledast$ depends on the choice of the quantization scheme (e.g.,
the ordering recipe) on the $C$ sector.

The $\circledast$-product [18]

$\displaystyle A\circledast
B=AB+\frac{i\hbar}{2}\big{[}\\{A,B\\}+\sigma(A,B)\big{]}$ (4.4)

acts only on the classical part of hybrid variables. The symmetric binary
operation $\sigma$ reflects the quantization scheme on the $C$ sector prior to
taking the classical limit. Using (4.4) in (4.2), we obtain the bracket in a
more familiar form

$\displaystyle\begin{split}\\{\\![A,B]\\!\\}=\frac{1}{i\hbar}[A,B]&+\frac{1}{2}\big{[}\\{A,B\\}-\\{B,A\\}\big{]}\\\
&+\frac{1}{2}\big{[}\sigma(A,B)-\sigma(B,A)\big{]}~{}.\end{split}$ (4.5)

For $\sigma=0$, the bracket reduces to that proposed by Aleksandrov [2],
Gerasimenko [3], and Boucher and Traschen [4].

A general formula for $\sigma$ resulting from familiar quantization schemes on
the $C$ sector is given by

$\displaystyle\begin{split}\sigma({A},{B})&=a~{}\frac{\partial A}{\partial
x_{C}}\,\frac{\partial B}{\partial x_{C}}+b~{}\frac{\partial A}{\partial
p_{C}}\,\frac{\partial B}{\partial p_{C}}\\\ &\qquad+c\left(\frac{\partial
A}{\partial x_{C}}\,\frac{\partial B}{\partial p_{C}}+\frac{\partial
A}{\partial p_{C}}\,\frac{\partial B}{\partial x_{C}}\right)~{},\end{split}$
(4.6)

where the constants $(a,b,c)$ reflect the choice of quantization. For example,
$(0,0,0)$ reflects the Weyl ordering associated with the Wigner phase space
distribution, while $(1,1,0)$ reflects a certain parameterization of the anti-
normal ordering associated with the Husimi distribution.

Now a hybrid equation of motion for hybrid variables in the Heisenberg picture
can be defined as

$\displaystyle\frac{dA}{dt}=\\{\\![A,H]\\!\\}+\frac{\partial A}{\partial
t}~{}.$ (4.7)

Once a choice for $\sigma$ (or the quantization scheme prior to the classical
limit) is specified, dynamical evolution is unambiguously defined. Different
$\sigma$’s will produce different evolutions.

As an illustrative example, take

$\displaystyle\eta=p_{C}\,p_{Q}$ (4.8)

to be a hybrid variable of interest. Notice that the time derivative of the
classical part of $\eta$ does not commute with quantum variables. Given the
harmonic oscillator Hamiltonian (3.1), we have

$\displaystyle\dot{\eta}_{C}=\dot{p}_{C}=\\{p_{C},H\\}=-k(x_{C}-x_{Q})~{},$
(4.9)

and thus the ordering of $\frac{d}{dt}\left(\eta_{Q}\eta_{C}\right)$ is
ambiguous. To see the problem, let us naively calculate $\dot{\eta}$ in two
different orderings

$\displaystyle\begin{split}\left(\dot{\eta}\right)_{1}&=\dot{p}_{C}\,p_{Q}+p_{C}\,\dot{p}_{Q}\\\
&=k\left(x_{Q}\,p_{Q}-x_{C}\,p_{Q}-x_{Q}\,p_{C}+x_{C}\,p_{C}\right)~{},\end{split}$
(4.10a)
$\displaystyle\begin{split}\left(\dot{\eta}\right)_{1}&=\dot{p}_{C}\,p_{Q}+p_{C}\,\dot{p}_{Q}\\\
&=k\left(p_{Q}\,x_{Q}-x_{C}\,p_{Q}-x_{Q}\,p_{C}+x_{C}\,p_{C}\right)~{},\end{split}$
(4.10b)

It could be suggested to symmetrize the product of $x_{Q}$ and $p_{Q}$, and
thus get rid of the ambiguity. However, such symmetrization would be an ad hoc
choice of operator ordering on the $Q$ sector. A consistent hybrid bracket is
necessary.

In terms of the hybrid bracket (4.5), time evolution of $\eta$ is given by

$\displaystyle\begin{split}\dot{\eta}&=\\{\\![\eta,H]\\!\\}\\\
&=k\left(\frac{x_{Q}\,p_{Q}+p_{Q}\,x_{Q}}{2}-x_{C}\,p_{Q}-x_{Q}\,p_{C}+x_{C}\,p_{C}\right.\\\
&\qquad\qquad\qquad\qquad\qquad\qquad+\left.\frac{i\hbar}{2}\sigma(x_{C},p_{C})\right)~{}.\end{split}$
(4.11)

Compare (4.11) to (4.10). While a symmetrization of (4.10) gets rid of the
ordering ambiguity (in an ad-hoc manner), it completely misses out on the
$\sigma$ term present in (4.11). To reiterate, $\sigma$ is connected the
operator ordering on the $C$ sector, not the $Q$ sector, before taking the
classical limit.

As promised, different forms of $\sigma$ will have different evolutions. In
the simple example of $\eta=p_{C}\,p_{Q}$, using the definition (4.6) of
$\sigma$, $\dot{\eta}$ becomes

$\displaystyle
k\left(\frac{x_{Q}\,p_{Q}+p_{Q}\,x_{Q}}{2}-x_{C}\,p_{Q}-x_{Q}\,p_{C}+x_{C}\,p_{C}+\frac{i\hbar}{2}c\right)~{}.$
(4.12)

Here, $c$ is a constant determined by the original quantization scheme on the
$C$ sector before the classical limit, as discussed below Eq. (4.6). Finally,
the time dependence of $\eta$ can be found by plugging in the solutions for
$(x_{C},p_{C},x_{Q},p_{Q})$. The classical variables $x_{C}$ and $p_{C}$ have
been explicitly calculated in Eqs. (3.13) and (3.14). The equations for
$(x_{Q},p_{Q})$ can be found similarly.

The argument presented in this section is not specific to the harmonic
oscillator, of course. The procedure is general and applicable to any
Hamiltonian and hybrid variable provided they belong to an associative
subalgebra with the hybrid composition product $\circledast$.

##  Conclusion

The example of the quantum-classical harmonic oscillator presented here
exhibits basic features of hybrid systems. The treatment used is applicable to
general Hamiltonians.

Quantum backreaction in interacting quantum-classical systems can, in
principle, be described through the framework of Hamiltonian hybrid dynamics.
The example of the quantum-classical harmonic oscillator studied here
demonstrates this backreaction: classical variables grow to be uncertain by
interacting with quantum ones. The uncertainty in classical variables is
manifested in the nonvanishing commutator of classical variables.

Along with Anderson [5], we note that despite the fuzziness of classical
variables through interaction, they are still classical in the sense that
their uncertainty is entirely due to the effect on them of truly quantum
variables. That effect is depicted in Fig. 1: a commutator of classical
variables vanishes initially, but evolves to be nonzero at later times. The
figure also illustrates that the more massive the classical system is, the
longer it remains certain.

An important property made plain by the analysis presented here is that the
details of a hybrid bracket are irrelevant for the dynamics of pure variables.
All brackets obeying the reduction and bilinearity properties (Eqs. (2.1) and
(2.2)) are equivalent as far as pure variables are concerned. This covers a
wide class of problems where the only nontrivial hybrid variable is the
interaction Hamiltonian.

Dynamics of hybrid variables requires the use of a fully-defined hybrid
bracket. Without such a bracket, we run into ordering ambiguities and miss
extra effects resulting from the hybrid nature of the system. As shown in
[18], these extra effects arise from the specific quantization scheme (and its
associated phase space distribution) used in deriving the bracket.

## Acknowledgements

This research was funded partly by an NSERC Discovery Grant.

## References

  * [1] John Maddox.  Classical and quantum physics mix.  Nature, 373, 1995.
  * [2] I. V. Aleksandrov.  The Statistical Dynamics of a System Consisting of a Classical and a Quantum Subsystem.  Zeitschrift für Naturforschung A, 36, 1981.
  * [3] V. I. Gerasimenko.  Dynamical equations of quantum-classical systems.  Theor Math Phys, 50, 1982.
  * [4] Wayne Boucher and Jennie Traschen.  Semiclassical physics and quantum fluctuations.  Phys. Rev. D, 37, 1988.
  * [5] Arlen Anderson.  Quantum Backreaction on "Classical" Variables.  Phys. Rev. Lett., 74, 1995.
  * [6] Oleg V. Prezhdo and Vladimir V. Kisil.  Mixing quantum and classical mechanics.  Phys. Rev. A, 56, 1997.
  * [7] Oleg V. Prezhdo.  A quantum-classical bracket that satisfies the Jacobi identity.  J. Chem. Phys., 124, 2006.
  * [8] Hans-Thomas Elze.  Linear dynamics of quantum-classical hybrids.  Phys. Rev. A, 85, 2012.
  * [9] Denys I. Bondar, François Gay-Balmaz, and Cesare Tronci.  Koopman wavefunctions and classical–quantum correlation dynamics.  Proc. R. Soc. A., 475(2229):20180879, September 2019.
  * [10] Carlos Barceló, Raúl Carballo-Rubio, Luis J. Garay, and Ricardo Gómez-Escalante.  Hybrid classical-quantum formulations ask for hybrid notions.  Phys. Rev. A, 86, 2012.
  * [11] Hans-Thomas Elze.  Four questions for quantum-classical hybrid theory.  J. Phys.: Conf. Ser., 361, 2012.
  * [12] Hans-Thomas Elze.  Quantum-classical hybrid dynamics - a summary.  J. Phys.: Conf. Ser., 442, 2013.
  * [13] L. L. Salcedo.  Absence of classical and quantum mixing.  Phys. Rev. A, 54, 1996.
  * [14] J. Caro and L. L. Salcedo.  Impediments to mixing classical and quantum dynamics.  Phys. Rev. A, 60, 1999.
  * [15] Debendranath Sahoo.  Mixing quantum and classical mechanics and uniqueness of Planck’s constant.  J. Phys. A: Math. Gen., 37, 2004.
  * [16] L. L. Salcedo.  Statistical consistency of quantum-classical hybrids.  Phys. Rev. A, 85, 2012.
  * [17] V. Gil and L. L. Salcedo.  Canonical bracket in quantum-classical hybrid systems.  Phys. Rev. A, 95, 2017.
  * [18] Mustafa Amin and Mark A. Walton.  Quantum-Classical Dynamical Brackets.  2020\.  arXiv:2009.09573 [quant-ph].