# A non-local phase field model of Bohm’s quantum potential

Roberto Mauri111Author to whom correspondence should be addressed. Electronic
mail<EMAIL_ADDRESS>Department of Civil and Industrial Engineering,
Laboratory of Multiphase Reactive Flows,
Università di Pisa, I-56126 Pisa, Italy

###### Abstract

Assuming that the free energy of a gas depends non-locally on the logarithm of
its mass density, the body force in the resulting equation of motion consists
of the sum of density gradient terms. Truncating this series after the second
term, Bohm’s quantum potential and the Madelung equation are identically
obtained, showing explicitly that some of the hypotheses that led to the
formulation of quantum mechanics admit a classical interpretation based on
non-locality.

##  1 Introduction

In 1927 Madelung [1] showed that the Schrodinger equation for one-electron
problems can be transformed into hydrodynamical equations. These equations
describe the time evolution of the density and velocity fields of an ideal
(i.e., non-viscous) gas subjected to the action of both classical and quantum
potentials. This latter represents the novelty of the Madelung equation; it is
generally referred to as the Bohm potential [2, 3] and is a central concept of
the de Broglie-Bohm formulation of quantum mechanics (see review in [4]).
Accordingly, the Schrodinger equation can be interpreted as a Hamilton-Jacobi
equation associated to the Madelung hydrodynamic equations, so that Madelung
could conclude [1] that there is thus a chance of erecting the quantum theory
of the atom on this basis. In fact, the Madelung equations are the basis for a
full description of quantum fluids [5].

Since then, many researchers have observed that, as the Bohm quantum potential
derives from the Laplacian term of the Schrodinger equation, quantum effects
must the consequence of some sort of non-local interaction. The non-local
nature of quantum mechanics is further stressed by the fact that, due to the
instantaneous nature of the collapse of the wavefunction, possible actions
performed in a certain space region can induce immediate changes in distant
regions [4, 7]. Nevertheless, a clear explanation of how the Schrodinger
equation derives from the non-locality of the force field is still missing.

In this work, we propose to fill this gap, starting from a non-local
characterization of an ideal, i.e., non-viscous, fluid. Then, considering
particle density as a phase field, we show that its evolution is described by
the Madelung hydrodynamic equations.

##  2 From non-locality to Schrödinger equation

Consider an ideal fluid in isothermal conditions, that is a system composed of
a large number of identical particles of mass m that do not interact with each
other except for instantaneous thermalizing collisions. Its governing
equations, expressing the conservation of mass and momentum, read [8]:

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\mathbf{v}\right)=0,\hskip 40.00006pt$ i.e.,
$\displaystyle\hskip
30.00005pt\frac{d\rho}{dt}=-\rho\left(\nabla\cdot\mathbf{v}\right),$ (1)
$\displaystyle\frac{\partial\left(\rho\mathbf{v}\right)}{\partial
t}+\nabla\cdot\left(\rho\mathbf{v}\mathbf{v}\right)=\rho\mathbf{F},\hskip
10.00002pt$ i.e., $\displaystyle\hskip
30.00005pt\rho\frac{d\mathbf{v}}{dt}=\rho\mathbf{F},$ (2)

with $dA/dt=\partial A/\partial t+\mathbf{v}\cdot\nabla A$ denoting the
material derivative. Here, $\rho$ is the density, $\mathbf{v}$ is the fluid
velocity, while $\mathbf{F}$ is a force per unit mass, so that
$\rho\mathbf{F}$ is a body force, i.e., a force per unit volume. Note that
$\rho$ and $\mathbf{v}$ are mean values, defined within an elementary point
volume, assuming local equilibrium. In addition, $\rho$ is proportional to the
probability distribution to find a particle in a unit volume located in
$\mathbf{r}$ at time $t$. Now, considering that the probability distribution
for a combination of two mutually independent subsystems is equal to the
product of the distribution functions of the two subsystems separately, we
reach the conclusion that $\ln\rho$ must be a scalar and additive integral of
motion, and therefore it must be proportional to the energy (see Landau et al.
[9]. Therefore, we may assume that the free energy per particle mass,
$\widehat{f}\left(\mathbf{r},t\right)$ , of a gas at constant temperature $T$
has the following form:

$\widehat{f}\left(\mathbf{r},t\right)=\frac{kT}{m}\int{u\left(|\mathbf{r}-\mathbf{r}^{\prime}|\right)\ln{\rho}\left(\mathbf{r}^{\prime},t\right)d^{3}\mathbf{r}^{\prime}}+\phi\left(\mathbf{r},t\right),$
(3)

where $\phi\left(\mathbf{r},t\right)$ is the potential energy per unit mass
resulting from the action of an external conservative force, while
$u\left(x\right)$ is an interaction kernel between particles located at a
distance $x=|\mathbf{r}-\mathbf{r}^{\prime}|$, with the normalization
condition, $\int{u\left(x\right)d^{3}\mathbf{x}}=1$ . Dropping for convenience
the time dependence, this expression can be rewritten in the more convenient
form,

$\widehat{f}\left(\mathbf{r}\right)=\widehat{f}_{th}\left(\mathbf{r}\right)+\frac{kT}{m}\int{u\left(|\mathbf{r}-\mathbf{r}^{\prime}|\right)\left[\ln{\rho}\left(\mathbf{r}^{\prime}\right)-\ln{\rho}\left(\mathbf{r}\right)\right]d^{3}\mathbf{r}^{\prime}},$
(4)

where

$\widehat{f}_{th}\left(\mathbf{r}\right)=\frac{kT}{m}\ln\rho\left(\mathbf{r}\right)+\phi\left(\mathbf{r}\right)$
(5)

is the thermodynamic part, corresponding to the free energy of a uniform ideal
gas. Expanding in Taylor series,

$\ln\rho\left(\mathbf{r}+\mathbf{x}\right)=\ln\rho\left(\mathbf{r}\right)+\mathbf{x}\cdot\nabla\ln\rho\left(\mathbf{r}\right)+\frac{1}{2}\nabla\nabla:\ln\rho\left(\mathbf{r}\right)+\cdots$
(6)

and truncating the series after the second term (there can be no
$\nabla\ln\rho$ term, due the isotropy of the fluid), we find,

$\widehat{f}\left(\mathbf{r}\right)=\widehat{f}_{th}\left(\mathbf{r}\right)+\widehat{f}_{nl}^{(2)}\left(\mathbf{r}\right)$
(7)

where the leading-order correction term,

$\widehat{f}_{nl}^{(2)}\left(\mathbf{r}\right)=-\frac{1}{2}\frac{kT}{m}a^{2}\nabla^{2}\ln\rho\left(\mathbf{r}\right)$
(8)

is the non-local part, with

$a^{2}=-\int{x^{2}u\left(x\right)d^{3}\mathbf{x}}$ (9)

denoting the square of a characteristic length, $a$. Note the negative sign in
Eq. (9), revealing that particles are assumed to attract each other.

Now we define a free energy per unit volume, $\rho\widehat{f}$, so that the
total free energy is given by the following functional:

$F\left[\rho\left(\mathbf{r}\right)\right]=\int_{V}{\rho\left(\widehat{f}_{th}-\frac{1}{2}\frac{kT}{m}a^{2}\nabla^{2}\ln\rho\right)d^{3}\mathbf{r}},$
(10)

where $V$ is the total volume, that we assume to be infinite. Integrating by
parts with the assumption that $\rho\rightarrow 0$ exponentially as
$r\rightarrow\infty$, find:

$F\left[\rho\left(\mathbf{r}\right)\right]=\int_{V}{fd^{3}\mathbf{r}}=\int_{V}{\rho\left[\widehat{f}_{th}+\frac{1}{2}\frac{kT}{m}a^{2}\left(\nabla\ln\rho\right)^{2}\right]d^{3}\mathbf{r}},$
(11)

where $f$ denotes the effective free energy per unit volume.

At equilibrium, the free energy functional (11) is minimal, under mass
conservation constraint, $\int{\rho d^{3}\mathbf{r}}=M\text{ constant}$.
Therefore, applying the Euler-Lagrange equation, we see that the Legendre
multiplier defines the following generalized chemical potential,

$\mu=\frac{\delta f}{\delta\rho}=\frac{\partial
f}{\partial\rho}-\nabla\cdot\frac{\partial f}{\partial\nabla
f}=\mu_{th}-\frac{kT}{m}a^{2}\left[\nabla^{2}\ln\rho+\frac{1}{2}\left(\nabla\ln\rho\right)^{2}\right].$
(12)

Here, the first term on the RHS denotes the thermodynamic chemical potential
(i.e., Gibbs’ free energy per unit mass), i.e.,

$\mu_{th}=d\left(\rho\widehat{f}_{th}\right)/d\rho=\widehat{f}_{th}-\widehat{v}d\widehat{f}_{th}/d\widehat{v}=\widehat{f}_{th}+P\widehat{v},$
(13)

where $\widehat{v}=\rho^{-1}$ and $P=-d\widehat{f}_{th}/d\widehat{v}$ are the
specific volume and the pressure, respectively.

The second term on the RHS of Eq. (12) is the non-local contribution to the
generalized chemical potential, i.e.,

$\mu_{nl}^{(2)}=-\frac{kT}{m}a^{2}\left[\nabla^{2}\ln\rho+\frac{1}{2}\left(\nabla\ln\rho\right)^{2}\right]$
(14)

The equations of motion can be determined applying a variational principle
[10] to derive the Euler Eq. (2) for a compressible, inhomogeneous fluid where
$\mathbf{F}$ is a body force, driven by density gradients in the fluid. The
same result can be obtained more heuristically by applying Noether’s theorem
[11] or by the least action principle [12]. At the end, we find Eqs. (1) and
(2), with the following force per unit mass:

$\mathbf{F}=-\nabla\mu=-\nabla\mu_{th}-\nabla\mu_{nl}^{(2)}.$ (15)

Let us consider the two terms on the RHS separately. On one hand we have:

$\mathbf{F}_{th}=-\nabla\mu_{th},\hskip 20.00003pt\text{with}\hskip
20.00003pt\mu_{th}=\frac{kT}{m}\left(\ln\rho+1\right)+\phi,$ (16)

and therefore, considering that $P=\frac{kT}{m}\rho$,

$\rho\mathbf{F}_{th}=-\rho\nabla\mu_{th}=-\nabla P-\rho\nabla\phi.$ (17)

On the other hand, from Eqs. (14) and (15), we define a non-local reversible
body force, which is usually referred to as the Korteweg force:

$\mathbf{F}_{nl}^{(2)}=-\nabla\mu_{nl}^{(2)}=\frac{kT}{m}a^{2}\nabla\left[\nabla^{2}\ln\rho+\frac{1}{2}\left(\nabla\ln\rho\right)^{2}\right].$
(18)

The non-local potential (14) can be expressed in the following equivalent
form:

$\mu_{nl}^{(2)}=-2\dfrac{kT}{m}a^{2}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}.$
(19)

In fact:

$\displaystyle\nabla^{2}\ln\rho+\frac{1}{2}\left(\nabla\ln\rho\right)^{2}=\frac{1}{\sqrt{\rho}}\left[\sqrt{\rho}\nabla^{2}\ln\rho+\frac{\nabla\rho\cdot\nabla\ln\rho}{2\sqrt{\rho}}\right]=$
$\displaystyle=\frac{1}{\sqrt{\rho}}\nabla\cdot\left(\sqrt{\rho}\nabla\ln\rho\right)=\frac{1}{\sqrt{\rho}}\nabla\cdot\left(\frac{\nabla\rho}{\sqrt{\rho}}\right)=2\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}.$
(20)

When $a$ is the thermal de Broglie wavelength [13],

$a=\frac{\hbar}{\sqrt{4mkT}},$ (21)

then $\mu_{nl}$ reduces to Bohm’s quantum potential [2, 3], i.e.,

$\mu_{nl}^{(2)}=Q=-\frac{\hbar^{2}}{2m^{2}}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}.$
(22)

Finally, summarizing, from Eq. (2), the force balance per unit mass yields the
Madelung equation:

$\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}=\mathbf{F}=-\nabla V-\nabla Q,$ (23)

with

$V=\widehat{f}_{th}=\frac{kT}{m}\ln\rho+\phi,\hskip 10.00002pt\text{and}\hskip
10.00002pt\nabla V=\nabla P/\rho+\nabla\phi.$ (24)

Here, $V$ denotes the thermodynamic potential energy (5), while $Q=\mu_{nl}$
is its non-local counterpart, that is the Bohm potential.

Since the fluid is frictionless, the flow is irrotational (i.e., with zero
curl), and the fluid velocity can be expressed as the gradient of a real
scalar function, namely the velocity potential, $S$, that is,

$\mathbf{v}=\nabla S.$ (25)

Substituting (25) into (23) we see that $S$ satisfies the following Hamilton-
Jacobi equation, describing the motion of a particle of mass $m$ in a
conservative force field, characterized by a classical potential energy, $V$,
and the quantum potential field, $Q$ (see Eqs. (22) and (23)),

$\frac{\partial S}{\partial t}+\frac{1}{2}\left(\nabla
S\right)^{2}+V-\frac{\hbar^{2}}{2m^{2}}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}.$
(26)

This equation, together with the continuity equation (1),

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\nabla S\right)=0,$ (27)

can be obtained using the Madelung transform,

$\psi=\sqrt{\rho}e^{iS},$ (28)

and imposing that the complex function $\psi$ satisfies the Schrödinger
equation,

$i\hbar\frac{\partial\psi}{\partial
t}=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+V\psi,$ (29)

with $\frac{\hbar^{2}}{2m}=2kTa^{2}$. Eq. (23) can also be written as a
Bernoulli equation as follows:

$m\frac{\partial\mathbf{v}}{\partial t}=-\nabla
E=-\nabla\left(\dfrac{1}{2}mv^{2}+mV+mQ\right),$ (30)

where the RHS denotes the gradient of the total energy, $E$, that is the sum
of kinetic energy, potential energy, and quantum potential energy. With some
manipulation, we can write Eq. (30) in the following equivalent form,

$\rho\left(\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}\right)=\rho\mathbf{F}=-\rho\nabla
V-\nabla\cdot\mathbf{P}_{Q},\hskip
10.00002pt\mathbf{P}_{Q}=-\frac{\hbar^{2}}{4m}\rho\nabla\nabla\ln\rho$ (31)

with $\mathbf{P}_{Q}$ denoting the quantum stress tensor. Finally, it must be
noted that the non-local term is not a small correction to its thermodynamic
counterpart, as in Eq. (18) $\nabla^{2}\rho\approx\rho/a^{2}$, so that the
non-local potential $Q$ is not necessarily smaller than the local,
$\mathcal{O}(kT/m)$ term.

Equation (31) (or, equivalently, Eqs. (23) and (31)), together with Eq. (1),
are known as the superfluid hydrodynamic equations [5].

##  3 Higher-order terms

Considering the higher-order terms in the expansion (6), we obtain the
following expansion for the free energy,

$\widehat{f}=\frac{kT}{m}\sum_{n=0}^{\infty}(-1)^{n}\frac{c_{2n}}{(2n)!}a^{2n}\nabla^{2n}\ln\rho,$
(32)

where $c_{0}=1$ and $c_{2}=1$ (this is the definition (9) of $a^{2}$, while,

$c_{2n}=(-1)^{n}\frac{1}{a^{2n}}\int_{V}{x^{2n}u\left(x\right)d^{3}\mathbf{x}};\hskip
5.0ptn=2,3,\cdots$ (33)

are the normalized central moments of the $u(x)$ distribution.

The total free energy then becomes, after integrating by parts:

$\displaystyle
F\left[\rho\left(\mathbf{r}\right)\right]=\int_{V}fd^{3}\mathbf{r}=\frac{kT}{m}\sum_{n=0}^{\infty}(-1)^{n}\frac{c_{2n}}{(2n)!}a^{2n}\int_{V}{\rho\left(\nabla^{2n}\ln\rho\right)d^{3}\mathbf{r}}=$
(34)
$\displaystyle=F_{th}\left(\mathbf{r}\right)-\frac{kT}{m}\sum_{n=1}^{\infty}(-1)^{n}\frac{c_{2n}}{(2n)!}a^{2n}\int_{V}{\nabla\ln\rho\cdot\left(\nabla\nabla^{2n-2}\rho\right)d^{3}\mathbf{r}},$
(35)

where $F_{th}=\frac{kT}{m}\int{\rho\ln\rho d^{3}\mathbf{r}}$. Naturally,
truncating the series (34) after $n=1$, we find Eqs. (10) and (11). Now,
minimizing this functional with the constraint of mass conservation we obtain:

$\mu=\frac{\delta
f}{\delta\rho}=\sum_{m=0}^{\infty}\frac{\left(-1\right)^{m}}{m!}\nabla^{m}\frac{\partial
f}{\partial\nabla^{m}\rho}=\mu_{th}+\sum_{n=1}^{\infty}\mu_{nl}^{(2n)},$ (36)

where $\mu_{th}$ is the thermodynamic chemical potential (13), while the non-
local potentials $\mu_{nl}^{(2n)}$ read:

$\mu_{nl}^{(2n)}=\frac{kT}{m}(-1)^{n}\frac{c_{2n}}{(2n)!}a^{2n}\left[\frac{1}{\rho}\nabla^{2n}\rho+\nabla^{2n}\ln\rho\right].$
(37)

In particular, $c_{2}=1$ and $\mu_{nl}^{(2)}=Q$. Finally, the body force
appearing in the Bernoulli equation (2) is again equal to the expression (15),
with the Korteweg force now being equal to the sum of the gradients of the
non-local potential (37),

$\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}=\mathbf{F}=-\nabla V-\nabla
Q-\sum_{n=2}^{\infty}\nabla\mu_{nl}^{(2n)}.$ (38)

The last term on the RHS of Eq. (38) represents a presumably small correction
to the Madelung hydrodynamic equation (23).

##  4 Comments and conclusions

In this work, Bohm’s quantum potential is derived assuming non-local effects,
namely a non-local dependence of the free energy of the system on the
logarithm of its mass density through an interaction kernel, $u(x)$. This
implies that the body force in the Madelung hydrodynamic equations of quantum
fluids consists of a density gradient expansion (see Eq. (38)), truncated
after the second term. Clearly, that means that considering all the hierarchy
of density gradient terms, the interaction kernel $u(x)$ leads to very small
corrections to the Schrödinger equation. For example, assume that $u(x)$
consists of a hard-core, short-range repulsion term of range $d$ and a long-
range attractive, exponentially decaying potential of range $\ell$, i.e.,

$u\left(x\right)=\left\\{\begin{split}n_{0}\hskip 30.00005pt&(x<d)\\\
-n_{1}\,exp(-\frac{x}{\ell})\hskip 10.00002pt&(x>d)\end{split}\right.$ (39)

in the limit $d/\ell=\mathcal{O}(\epsilon)$ , with $\epsilon\ll 1$, assuming
that $n_{0}d^{3}\approx n_{1}\ell^{3}=\mathcal{O}(\epsilon^{-1})$, the
normalization condition, $\int{u(x)d^{3}\mathbf{x}}=1$ , leads to
$n_{0}=9n_{1}\ell^{3}/d^{3}$, while Eq. (9) yields:
$(a/\ell)^{3}=96\pi(n_{1}\ell^{3})$, showing that
$a/\ell=n_{0}d^{3}=\mathcal{O}(\epsilon^{-1/2})$. In general, Eq. (33) gives:

$c_{2n}=(-1)^{n+1}\frac{\left(2n+2\right)!}{24}\,\left(\frac{\ell}{a}\right)^{2(n-1)},\hskip
10.00002ptn=1,2,\cdots,$ (40)

showing that, as $\ell\ll a$, the moments decrease as $n$ increases. In
particular, we confirm that $c_{2}=1$, while
$c_{4}=-30(\ell/a)^{2}=\mathcal{O}(\epsilon)$,
$c_{6}=1680(\ell/a)^{4}=\mathcal{O}(\epsilon^{2})$ , etc. Consequently,
assuming that the interaction kernel is given by Eq. (39), we see that the
additional terms, correcting the Madelung equation (26), are very small.

Concluding, we have shown that some of the hypotheses that led to the
formulation of quantum mechanics admit a simple classical interpretation,
based on non-locality, so that quantum mechanics could be interpreted as a
sort of mean field theory, that is a coarse-grained approximation of a more
microscopic physical reality. What remains to be seen is the physical nature
of the underlying non-local interaction field (see discussion by L. Smolin
[14]). In fact, although this hypothesis is not used in the mathematical
development of the model, it is essential for its physical understanding,
analogously to the background field hypothesis leading to the formulation of
stochastic mechanics (see Nelson [15]).

## References

  * [1] Madelung, E.: Quantum theory in hydrodynamical form. Zeit. f. Phys., 40, 322 (1927).
  * [2] Bohm D.: A suggested interpretation of the quantum theory in terms of “hidden” variables. I. Phys. Rev. 85, 166 (1952).
  * [3] Bohm, D.: A Suggested Interpretation of the Quantum Theory in Terms of “Hidden” Variables. II. Phys. Rev. 85, 180 (1952).
  * [4] Bohm, D., Hiley, B.J.: Non-locality and locality in the stochastic interpretation of quantum mechanics. Phys. Rep. 172, 93 (1989).
  * [5] Barenghi, C.F.F., Parker, N.G.: A Primer on Quantum Fluids. Springer, Berlin (2016).
  * [6] Holland, P.R.: The de Broglie-Bohm theory of motion and quantum field theory. Phys. Rep. 224, 95 (1993).
  * [7] Bassi, A., Ghirardi G-C.: Dynamical reduction model. Phys. Rep. 379, 257 (2003).
  * [8] De Groot, S.R., Mazur, P.: Non-Equilibrium Thermodynamics. North Holland, Amsterdam (1962).
  * [9] Landau, L.D., Lifshitz, E.M., Pitaevskii, L.P.: Statistical Physics I, Pergamon (1980).
  * [10] Lamorgese, A., Mauri, R., Sagis, L.M.C.: Modeling soft interface dominated systems: A comparison of phase field and Gibbs dividing surface models. Phys. Rep. 675, 1 (2017).
  * [11] D. Anderson, McFadden, G., Wheeler, A.: Diffuse-interface methods in fluid mechanics. Ann. Rev. Fluid Mech. 30, 139 (1998).
  * [12] Lowengrub,J., Truskinovsky, L.: Quasi-incompressible Cahn-Hilliard fluids and topological transitions. P. R. Soc. London A 454, 2617 (1998).
  * [13] Gaveau, B., Schulman, L.S.: Reconciling kinetic and quantum theory. Found. Phys. 50, 55 (2020).
  * [14] Smolin, L., Could quantum mechanics be an approximation to another theory? arXiv:quant-ph/0609109 (2016).
  * [15] Nelson, E.: Quantum Fluctuations. Princeton University Press (1985).


# Demonstration of an Integrated Terahertz Band-Stop Filter Using an Apodized
Bragg Grating

Ali Dehghanian  1,3 Walid Gomaa  2 Mohsen Haghighat1,3 Thomas Darcie1 and Levi
Smith1,3  1Department of Electrical and Computer Engineering, University of
Victoria, Victoria, BC, V8P 5C2 Canada
2Department of Engineering Physics, Military Technical College, Cairo, Egypt
3Centre for Advanced Materials and Related Technology (CAMTEC), University of
Victoria, 3800 Finnerty Rd, Victoria, BC, V8P 5C2, Canada<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

This paper presents the demonstration of an on-chip integrated Terahertz (THz)
Apodized Bragg grating (TABG) which functions as band-stop filter with a
center frequency of 0.8 THz and a bandwidth of 200 GHz. For experimentation,
we integrate the TABG into our THz System-on-Chip to enable wideband (DC - 1.5
THz) device characterization. Using this methodology, we measure the signal
transmission through the TABG and find the experimental results align with
simulation and theory provides a rejection of approximately 20 dB across the
stop-band.

##  1 Introduction

The terahertz (THz) gap spanning from 0.1 to 10 THz has garnered considerable
attention because of its unique capabilities across various fields such as
communication, sensing, imaging, and spectroscopy [1, 2, 3, 4, 5, 6]. However,
the progress in developing terahertz devices and systems has been challenging
due to the absence of integrated system components such as transmitters,
communication channels, receivers, and other passive devices such as filters.
Band-stop filters play a crucial role in many applications because they enable
the suppression of undesired frequencies within a specific spectral range. At
frequencies above 500 GHz, there are limited experimental results for planar
guided-wave band-stop filters [7], but there are a few key examples: in [7],
the authors demonstrate a band-stop filter with a center frequency of 0.6 THz
and a -3dB bandwidth of $\approx$160 GHz using a single conductor Goubau line
loaded with capacitivly-coupled $\lambda/2$ resonators. This filter behaves
well; however, Goubau-lines can be limiting when developing an integrated
system because the single conductor is incompatible with two-terminal elements
without necessitating a mode converter. Next, others have developed THz band-
stop filters using a dielectric Bragg grating with a center frequency of 0.31
THz and a -3dB bandwidth of $\approx$20 GHz[8]. Again, their filter works
well, but, as with most non-transverse electromagnetic (TEM) waveguides, the
finite single-mode bandwidth can lead to dispersion in a wideband system.
Next, in our prior work we loaded a transmisission line with split-ring
resonator elements [9]. Using this method we created a band-stop filter which
possessed a center frequency of 0.51 THz and a -3dB bandwidth of $\approx$100
GHz. These SRR filters are useful, but the design procedure is not systematic,
the filter roll-off rates are sub-optimal, and higher-order resonator modes
must be considered [9]. Lastly, we note that there are several non-planar
guided-wave Bragg gratings which have been demonstrated. In [10], a corrugated
circular waveguide was used as a Bragg grating which obtained a center
frequency was 0.14 THz and the -3dB bandwidth was $\approx$9 GHz. In [11], a
two-wire waveguide was periodically etched to fabricate the grating which
achieved a center frequency of 0.53 THz and the -3dB bandwidth was $\approx$1
GHz.

This work expands upon the aforementioned list of THz band-stop filters by
performing the experimental demonstration of an integrated Terahertz Apodized
Bragg Grating (TABG) proposed in [12]. Previously the TABG was analyzed via
simulations, but no experimental results were presented. Also, we compliment
the theory by adding periodic filter concepts from microwave engineering. The
TABG was inspired by apodized Fiber Bragg Gratings (FBGs) found in optical
communication systems which consist of alternating refractive indices (and
wave impedance) along the propagation direction. FBGs experience detrimental
side-lobes in the reflection spectrum which are proportional the FBG strength
(difference of refractive indices). The magnitude of the side-lobes can be
reduced by gradually tapering the FBG strength with an apodization profile.
The TABG analogously has a periodic modulation of the characteristic impedance
and similarly benefits from gradual tapering of the grating strength.

To characterize the TABG we use our integrated THz System-on-chip (TSoC)
platform which combines the THz transmitter, device-under-test (i.e., TABG),
and THz receiver onto a single wafer. The TSoC platform consists of planar
circuitry which is lithographically defined on an ultra-thin (1 $\mu$m) Si3N4
substrate to ensure signals exhibit low loss and low dispersion at THz
frequencies. The ultra-thin substrate is a key requirement to perform wideband
measurements otherwise radiation loss into the substrate becomes very
problematic. We have used the TSoC platform in previous works to investigate
several other THz components such as split-ring resonators [9], low-pass
filters [13], tapers [14], and power dividers [15]. Overall, we have found the
results from TSoC experiments to closely align with simulation and theory,
thus we expect our experimental methodology to be suitable for investigating
the TABG.

Figure 1: The experimental setup and TABG. (a) The optical path of the
transmitter is modulated by the optical chopper, while the receiver beam
passes through an optical delay line. The transmitter has a DC bias applied
(24V), and the receiver is connected to a lock-in amplifier. The bias lines
for the transmitter and receiver are separated by a gap that blocks DC
signals. (b) illustrates the TABG filter and its distinct sections.

##  2 Design

In this paper, we perform the experimental verification of the TABG shown in
Fig. 1 [12]. Specifically, we focus on a TABG where the stop-band is centered
at $f_{c}$ = 0.8 THz. The filter consists of N = 12 unit cell sections and 4
apodization sections which provide a gradual change in geometry and
characteristic impedance between the feedlines and the grating unit cells.
Figure 2 illustrates an annotated portion of the TABG. The TABG response is
defined by the reflections caused by characteristic impedance discontinuities
along the grating. Thus knowledge of the characteristic impedance of each
section is necessary to predict the filter response. In this work, the
characteristic impedance of each section is obtained by performing full-wave
simulations using ANSYS HFSS at 0.8 THz to ensure accurate values which
include the impact of the thin Si3N4 substrate. Table 1 tabulates the
dimensions and simulated characteristic impedances of each section. Next, the
grating period, $\Lambda$, required to obtain a specified center frequency,
$f_{c}$,is calculated using:

$\Lambda=\frac{c}{2f_{c}\sqrt{\varepsilon_{re}}},$ (1)

where $c$ is the speed of light, $\varepsilon_{re}$ is the effective relative
permittivity of the propagating mode. Again, from simulation (ANSYS HFSS), we
have found $\varepsilon_{re}\approx 1.3$ for the TABG CPS configuration, then
using (1), we find $\Lambda$ = 165 $\mu$m. Next, the fractional bandwidth,
$\Delta f/f_{c}$, is dependent on the characteristic impedance of adjacent
cells and is calculated by [16]:

$\frac{\Delta
f}{f_{c}}=\frac{4}{\pi}\sin^{-1}\left(\frac{Z_{n+1}-Z_{n}}{Z_{n+1}+Z_{n}}\right),$
(2)

using (2), with $Z_{0}$ and $Z_{U}$ (Table 1), we find $\Delta f=0.18$ THz.

Figure 2: Apodization sections and grating unit cell where $\Lambda=165\mu m$. Dimensions and characteristic impedances are found in Table 1. Table 1: TABG dimensions and characteristic impedances for Fig. 2 $n$ | 0 | 1 | 2 | 3 | 4 | U
---|---|---|---|---|---|---
Wn [$\mu$m] | 45 | 35 | 30 | 25 | 20 | 15
Sn [$\mu$m] | 70 | 80 | 85 | 90 | 95 | 100
Zn [$\Omega$] | 234 | 260 | 274 | 290 | 311 | 332

##  3 Theory

In microwave engineering, the Bragg filter of [12] is commonly referred to as
a periodic filter which can be described by a dispersion diagram and Bloch
impedance [17]. To calculate these quantities we construct an ABCD matrix of
the unit cell depicted in Fig. 2 as:

$\begin{split}\begin{bmatrix}A&B\\\
C&D\end{bmatrix}=&\begin{bmatrix}\text{cos}\beta\Lambda/4&jZ_{U}\text{sin}\beta\Lambda/4\\\
jY_{U}\text{sin}\beta\Lambda/4&\text{cos}\beta\Lambda/4\end{bmatrix}\\\
\cdot&\begin{bmatrix}\text{cos}\beta\Lambda/2&jZ_{0}\text{sin}\beta\Lambda/2\\\
jY_{0}\text{sin}\beta\Lambda/2&\text{cos}\beta\Lambda/2\end{bmatrix}\\\
\cdot&\begin{bmatrix}\text{cos}\beta\Lambda/4&jZ_{U}\text{sin}\beta\Lambda/4\\\
jY_{U}\text{sin}\beta\Lambda/4&\text{cos}\beta\Lambda/4\end{bmatrix},\end{split}$
(3)

where $\beta=\omega\sqrt{\varepsilon_{re}}/c$ is the propagation constant.

After constructing the ABCD matrix we can obtain the dispersion diagram and
Bloch impedance. The dispersion diagram is calculated using
$\beta\Lambda=\text{Imag}\\{\mbox{cosh}^{-1}[(A+D)/2]\\}$ and the result is
plotted in Fig. 3. This figure illustrates the stopband associated with the
infinite periodic structure. We find that the filter has a center frequency at
$f_{c}$ = 0.8 THz and a bandwidth of $\Delta f\approx 0.18$ THz which is in
agreement with (2). Next, we calculate the Bloch impedance using
$Z_{B}=B/\sqrt{A^{2}-1}$ then calculate the reflection coefficient between an
infinite periodic filter and a load impedance from:

$\Gamma=\frac{Z_{L}-Z_{B}}{Z_{L}+Z_{B}}.$ (4)

We use (4) to illustrate the impact of apodization. If we negate the
apodization sections, then $Z_{L}=Z_{0}$ which will result in a mismatch
causing a less desirable reflection spectrum. Alternatively, if we add an
apodization section then we can significantly improve the matching and reduce
the reflections. We model apodization section by cascading 4 ABCD matrices to
represent each apodization cell in Fig. 2. The result of this procedure is
plotted in Fig. 4 where we see the introduction of apodization sections
improves the reflection response of the filter and we expect to see a sharp
roll-off at the band-edge of the transmitted signal in our experimental
results.

Figure 3: Dispersion diagram for the TABG. Figure 4: Reflection coefficient at
the load of the TABG with and without an apodization section.

##  4 Simulation

A full-wave frequency domain simulation was performed using ANSYS HFSS to
characterize the scattering parameters of the TABG (Fig. 5). In the
simulation, the material parameters are given by $\varepsilon_{r}$ = 7.6,
$\sigma_{Si_{3}N_{4}}=0$, $\mu_{r}$ = 1, tan $\delta_{e}$ = 0.00526 for the
Si3N4 substrate[18] and $\sigma_{Au}=4.1\times 10^{7}$ S/m for the gold
conductors. All geometric parameters are the same as the fabricated device
(see Fig. 2 and Table 1). The simulated center frequency was found to be
$f_{c}$ = 0.815 THz and the -3 dB bandwidth was found to be 0.22 THz which
illustrates reasonable agreement between theory and simulation. We do not
expect perfect agreement between Fig. 4 ($|\Gamma|$) and Fig. 5 ($|S_{11}|$)
because the reflection coefficient of (4) calculates the reflection between an
infinite periodic structure and a constant load impedance, whereas the
simulation accounts for the finite length grating and the frequency-dependent
characteristic impedance. The insertion loss below the Bragg regime is less
than 2 dB, whereas above the Bragg regime, the insertion loss increases from 5
dB at 0.95 THz up-to 11 dB at 1.5 THz. The increased insertion loss originates
from diffractive grating radiation [19]. This effect is discussed and
illustrated in Appendix A.

Figure 5: Simulated S-Parameters. $f_{c}$ = 0.815 THz and the -3 dB bandwidth
is 0.22 THz.

##  5 Methods

To perform the experimental characterization of the TABG we use a modified THz
Time Domain Spectroscopy (THz-TDS) setup as shown in Fig. 1(a). An optical
pulse train is generated by a 780 nm, 20 mW, 80 femtosecond laser and then is
divided into two beams that are directed towards the transmitter and receiver
photoconductive switches (PCSs) which are thin films of LT-GaAs measuring 70
µm × 40 µm × 1.8 µm. The procedure to fabricate a grid of PCS’s is described
in [20, 13]. For the transmitter, the optical beam passes through a optical
chopper before being focused onto the PCS. The receiver beam passes through a
mechanical delay line before getting focused onto the PCS. The transmitter has
a bias voltage applied (24V) and the receiver is connected to a lock-in
amplifier. Translation of the mechanical delay line and measurement of the
lock-in amplifier current reconstructs a signal proportional to the
transmission characteristics of the TABG.

Figure 6: TSoC with TABG. (a) illustrates a LT-GaAs PCS (Tx) on a gold
transmission line on a 1µm-thin Si3N4 membrane. In (b), the TABG is displayed.
(c) shows the LT-GaAs PCS receiver (Rx), coupling section, and the DC bias
lines.

Figure 6(a) illustrates the transmitter which is placed on top of a
lithographically-defined 200 nm gold CPS TL situated on a 1 µm Si3N4 thin
layer. We use a thin Si3N4 membrane as the substrate to enable the
transmission of sub-picosecond pulses over centimeter-scale distances. If a
thicker substrate were used, radiation from the substrate would cause
significant losses and dispersion [21]. We have successfully used this method
in several other works [9, 13, 14, 15]. After the transmitter PCS, we taper
the CPS TL to achieve a wider cross-section (W = 45 µm and S = 70 µm) which
has lower attenuation [14]. In Fig. 6(b), the TABG structure is displayed
under a microscope. Fig. 6(c) portrays the receiver section, which comprises
of a DC block, DC bias lines for the Tx PCS, and bias lines to connect the Rx
PCS to the lock-in amplifier.

##  6 Results and Discussion

Figure 7: The experimental result of the received THz-bandwidth pulse after
propagation through the TABG. (a) The simulated temporal response (ANSYS
HFSS). (b) The experimental temporal response. (c) The DFT of the experimental
and simulated temporal responses.

Figure 7(a-b) plot the simulated and experimental temporal response of the
received THz-bandwidth pulse after propagating through the TABG. Figure 7(c)
plots the spectral response obtained by applying the Discrete Fourier
Transform (DFT) to the temporal response. First, we note that the spectral
roll-off associated with both the simulation and experiment is an expected
consequence of a finite duration time-domain pulse. Next, we observe good
agreement between experiment and simulation which is observed by comparing the
stop-band center frequency and bandwidth. We do notice a discrepancy in the
stop-band rejection. For the experiment the rejection is $\approx$30 dB,
whereas for simulation it is $\approx$45 dB. The difference originates from a
few different mechanisms. First, our photolithographic fabrication process is
limited to a feature size of $\approx$2 $\mu$m, thus there will not be perfect
agreement between the simulated dimensions and the fabricated dimensions.
Next, the material parameters will differ between simulation and experiment.
While we expect the material parameters to be close, they will not match the
fabricated device and will result in differences. Also, it is possible that
another unexpected mode coupled across the filter experiencing less
attenuation. Lastly, we explain a discrepancy temporal response. The
simulation does not use a DC blocking capacitor, thus it contains a DC
component. The experiment requires a DC block, therefore we observe the
derivative behaviour (high-pass filtering). Note that the cut-off frequency is
below our frequencies of interest (0.2 THz) and it is not observed in the
spectral response.

##  7 Conclusion

This work presented the experimental validation of a Terahertz Apodized Bragg
Grating (TABG) [12]. The TABG was designed to have a center frequency of 0.8
THz and a bandwidth of 0.2 THz. These parameters were confirmed by experiment.
Apodization serves to reduce the side-lobes and increase the roll-off rate of
the filter which was observed in simulation and experiment. The results showed
that the filter had a minimum rejection of approximately 20 dB over the stop-
band.

Novelty points: (1) The first demonstration of a fully integrated
(transmitter, feedlines, TABG, receiver) on-chip Apodized Bragg filter in the
THz gap. (2) The first demonstrated the highest Bragg frequency achieved
within the THz gap using a quasi-TEM feedlines. (3) The first demonstration of
a Bragg grating constructed from guided-wave TLs on an ultra-thin Si3N4
membrane. (4) We verify that the theory for periodic filter is applicable the
grating design. (5) We illustrate the radiative loss mechanism for frequencies
above the Bragg bandgap for a CPS grating (Appendix A).

## Funding

We acknowledge the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC).

## Acknowledgments

This work made use of the 4D LABS core facility at Simon Fraser University
(SFU) supported by the Canada Foundation for Innovation (CFI), British
Columbia Knowledge Development Fund (BCKDF), and Pacific Economic Development
Canada (PacifiCan). We would like to acknowledge CMC Microsystems for the
provision of products and services that facilitated this research.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * [1] H.-J. Song and T. Nagatsuma, “Present and future of terahertz communications,” IEEE transactions on terahertz science and technology 1, 256–263 (2011).
  * [2] P. H. Siegel, “Terahertz technology in biology and medicine,” IEEE transactions on microwave theory and techniques 52, 2438–2447 (2004).
  * [3] E. Pickwell and V. Wallace, “Biomedical applications of terahertz technology,” Journal of Physics D: Applied Physics 39, R301 (2006).
  * [4] Z. Chen, X. Ma, B. Zhang, Y. Zhang, Z. Niu, N. Kuang, W. Chen, L. Li, and S. Li, “A survey on terahertz communications,” China Communications 16, 1–35 (2019).
  * [5] L. Zhu, S.-H. Shin, R. Payapulli, T. Machii, M. Motoyoshi, N. Suematsu, N. M. Ridler, and S. Lucyszyn, “3-d printed rectangular waveguide 123–129 ghz packaging for commercial cmos rfics,” IEEE Microwave and Wireless Technology Letters (2023).
  * [6] M. C. Beard, G. M. Turner, and C. A. Schmuttenmaer, “Terahertz spectroscopy,” (2002).
  * [7] J. Cabello-Sánchez, V. Drakinskiy, J. Stake, and H. Rodilla, “Capacitively-coupled resonators for terahertz planar-goubau-line filters,” IEEE Transactions on Terahertz Science and Technology 13, 58–66 (2023).
  * [8] W. Gao, W. S. L. Lee, C. Fumeaux, and W. Withayachumnankul, “Effective-medium-clad Bragg grating filters,” APL Photonics 6, 076105 (2021).
  * [9] L. Smith, V. Shiran, W. Gomaa, and T. Darcie, “Characterization of a split-ring-resonator-loaded transmission line at terahertz frequencies,” Optics Express 29, 23282–23289 (2021).
  * [10] T. Ma, K. Nallapan, H. Guerboukha, and M. Skorobogatiy, “Analog signal processing in the terahertz communication links using waveguide bragg gratings: example of dispersion compensation,” Opt. Express 25, 11009–11026 (2017).
  * [11] J. Dong, A. Tomasino, G. Balistreri, P. You, A. Vorobiov, E. Charette, B. Le Drogoff, M. Chaker, A. Yurtsever, S. Stivala, M. A. Vincenti, C. De Angelis, D. Kip, J. Azana, and R. Morandotti, “Versatile metal-wire waveguides for broadband terahertz signal processing and multiplexing,” Nature Communications 13, 741 (2022).
  * [12] W. Gomaa and T. Darcie, “Design and simulation of terahertz apodized bragg grating using coplanar stripline transmission line a 1 $\mu$m-thin membrane,” in _Journal of Physics: Conference Series,_ vol. 2304 (IOP Publishing, 2022), p. 012015.
  * [13] W. Gomaa, L. Smith, V. Shiran, and T. Darcie, “Terahertz low-pass filter based on cascaded resonators formed by cps bending on a thin membrane,” Optics Express 28, 31967–31978 (2020).
  * [14] L. Smith, W. Gomma, H. Esmaeilsabzali, and T. Darcie, “Tapered transmission lines for terahertz systems,” Optics Express 29, 17295–17303 (2021).
  * [15] W. Gomaa, R. L. Smith, H. Esmaeilsabzali, and T. E. Darcie, “Terahertz power divider using symmetric cps transmission line on a thin membrane,” IEEE Access 8, 214425–214433 (2020).
  * [16] S. Orfanidis, _Electromagnetic Waves and Antennas_ (Online, 2016).
  * [17] D. Pozar, _Microwave Engineering_ (Wiley, 2011), 4th ed.
  * [18] G. Cataldo, J. A. Beall, H.-M. Cho, B. McAndrew, M. D. Niemack, and E. J. Wollack, “Infrared dielectric properties of low-stress silicon nitride,” Optics Letters 37, 4200–4202 (2012).
  * [19] P. Cheben, R. Halir, J. H. Schmid, H. A. Atwater, and D. R. Smith, “Subwavelength integrated photonics,” Nature 560 (2018).
  * [20] R. D. V. Rios, S. Bikorimana, M. A. Ummy, R. Dorsinville, and S.-W. Seo, “A bow-tie photoconductive antenna using a low-temperature-grown gaas thin-film on a silicon substrate for terahertz wave generation and detection,” Journal of Optics 17, 125802 (2015).
  * [21] R. Smith and T. Darcie, “Demonstration of a low-distortion terahertz system-on-chip using a cps waveguide on a thin membrane substrate,” Optics express 27, 13653–13663 (2019).

## Appendix A - Grating radiation above the stopband

When operating at frequencies above the Bragg bandgap there is an increased
insertion loss which originates from diffractive grating radiation. Figure 8
illustrates this concept where the black arrows illustrate the direction of
the Poynting vector. Note that the color scale is logarithmic and covers three
orders of magnitude. Below the bandgap, Fig. 8(a), the Poynting vector
illustrates a power flow from the left to right with no extra radiative loss.
Within the bandgap, Fig. 8(b), the incident signal from is reflected such that
minimal power is transmitted. Above the bandgap, Fig. 8(c), a portion of the
reflected wave is back–scattered by grating in the tangential direction which
contributes to the additional insertion loss.

Figure 8: Simulated electric field intensity profiles. Black arrows illustrate
the Poynting vector direction. a) Below the stopband, f = 0.6 THz. b) In the
stopband, f = 0.8 THz. c) Above the stopband, f = 1.0 THz.


# Quantum Compiler Optimizations

Jeff Booth<EMAIL_ADDRESS>
###### Abstract

A quantum computer consists of a set of quantum bits upon which operations
called gates are applied to perform computations. In order to perform quantum
algorithms, physicists would like to design arbitrary gates to apply to
quantum bits. However, the physical limitations of the quantum computing
device restrict the set of gates that physicists are able to apply. Thus, they
must compose a sequence of gates from the permitted gate set, which
approximates the gate they wish to apply - a process called _quantum
compiling_.

Austin Fowler proposes a method [2] that finds optimal gate sequences in
exponential time, but which is tractable for common problems. In this paper, I
present several optimizations to this algorithm. While my optimizations do not
improve its overall exponential behavior, they improve its empirical
performance by one to two orders of magnitude.

##  1 Background

In classical computing, we can generally rely on the correctness of hardware
because of the size of the circuit components. For example, if an atom on a
hard disk drive changed its spin orientation, or lost an electron, the hard
drive’s functionality would not be impaired because it takes many thousands of
atoms to represent and store a single bit of data. However, in quantum
computing, data is stored in quantum bits, which are represented by tiny
particles like trapped ions. These qubits are very easy to perturb,
potentially corrupting calculations based on them. Thus, we use redundancy, in
the form of error-correcting codes, to minimize the impact of individual
errors.

The Steane code is one representation of a quantum bit. It uses seven physical
qubits to represent one Steane code qubit, and can tolerate an arbitrary error
in one of the seven qubits. We can perform any desired operation on a Steane
code qubit by applying a combination of $H$ (Hadamard) and $T$ gate operations
[4]. $T$ gates are generally complicated to implement in quantum computing
hardware, so we seek to use a minimal number. For practical purposes, in
addition to $H$, we can use the Pauli X operator $X$, the Pauli Z operator
$Z$, the single qubit phase gate $S$, and its inverse $S^{\dagger}$
[Fowler2010]. The gates $H$, $X$, $Z$, $S$, and $S^{\dagger}$ generate a group
under multiplication, called the Clifford group. Thus, any sequence of gates
we choose will alternate between a member of the Clifford group and a $T$
gate. A $T^{\dagger}$ gate is also used in this implementation, bringing the
total number of non-identity gates in Fowler’s gate set to 25.

A _single-qubit quantum compiler_ finds sequences of gates which yield
matrices that are “close” to a gate we would like to apply to a quantum bit.
Each gate has a corresponding matrix that represents the operation it would
perform on a quantum bit. How close one gate is to another is given by the
“Fowler” distance:

$dist(U,U_{l})=\sqrt{\frac{2-\left|tr\left(U\cdot
U_{l}^{\dagger}\right)\right|}{2}}$ (1)

The longer the gate sequence is, the more closely it can approximate a desired
target gate that is not in the universal instruction set. However, a longer
gate sequence takes more time to compute on a real quantum computer,
increasing the probability of a computation error. An optimal quantum compiler
will find gate sequences which:

  1. 1.

have a minimal Fowler distance from the target gate.

  2. 2.

have a minimal length.

##  2 Fowler’s Algorithm

Austin Fowler presents an algorithm that iterates over sequences in order from
smallest to largest [Fowler2010]. For each sequence, it multiplies the
sequence gates’ matrices together to generate a $2\times 2$ unitary matrix
representing the complete operation that sequence would perform. The simple
brute-force iteration runs in time exponential in sequence length, since all
sequences of length $n$ are produced by appending all elements of the
universal instruction set to all sequences of length $n-1$.

To reduce the run time, Fowler’s algorithm intelligently skips redundant
sequences. The algorithm creates a list of unique sequences for all sequences
of length $N$. Then, for each sequence $S$ of length $N+1$, it searches for
sub-sequences of length $N$. If a sub-sequence $Y$ is not in the list, then it
is not unique. That means it performs the same operation as a sequence $V$
which is in the unique sequence list. Since $Y$ and $V$ are the same, then if
you were to replace $Y$ with $V$ in your sequence $S$, you would get a
sequence $W$ that does the same thing $S$ does.

Since the algorithm iterates over all sequences of length $N$, it will
encounter $W$ anyway (or it has already encountered $W$). Therefore, it should
skip sequence $S$. In fact, it should increment the sub-sequence $Y$ until it
is a unique sequence $U$. Fowler’s algorithm contains a tree lookup structure
which, for any given sequence, records the next unique sequence $U$. The
algorithm can determine what sequence to skip to by simply accessing this
tree. It still requires time exponential in sequence length, but interesting
results can now be obtained in mere days using consumer computer hardware.

As I will demonstrate later, Fowler’s tree data structure requires memory that
scales exponentially with the sequence length. Thus, the algorithm consists of
two stages:

  1. 1.

During the first stage, it builds the structure until it stores all unique
sequences up to length $W$, where $W=15$ for most of the experiments.

  2. 2.

After the structure is built, it enters the second stage, where it generates
sequences – and uses the structure to skip them – but it doesn’t add them to
the structure.

This dramatic change in behavior between stages explains some interesting
features in the following graphs. Also, it means I can only infer behavior on
longer sequences from behavior in the second stage, which explains my focus on
data produced during that stage.

##  3 Experimental Goal

In order to empirically measure the impact of my optimizations, I need a
consistent experimental goal to test on every version of the algorithm. For
this research, I chose to approximate the $\frac{\pi}{6}$ gate
$\exp(i\frac{\pi}{12}\sigma_{z})$ to $10^{-7}$ accuracy. Since this
approximation is currently very time-consuming, I can use it to empirically
evaluate the impact of my enhancements.

##  4 Existing Performance

In order to better understand the performance characteristics of the Fowler
algorithm, I modified its C source code to obtain performance-related
statistics. In this section, I present the data I gathered, along with some
explanations for unusual data features and speculations on how the statistics
should change for a meaningful performance improvement. Most of these
benchmarks ran on an Amazon Elastic Compute Cloud Medium computer, which
contains a 2-2.4 GHz processor and 3.75 GB of memory.

###  4.1 Code Profiling

I ran a profiler (gprof) to determine where the performance bottlenecks are.
Initially, I thought that memory accesses would dominate the program’s
runtime, because of the size of the data structures involved. However, the
program spends 92.49% of its time inside mathematical functions, meaning that
calculation is the dominant operation. In the first stage of the algorithm,
the program spends 98.5% of its time checking for unique matrices. In the
second stage, it spends most of its time multiplying gate matrices together to
calculate the matrix for a given sequence.

###  4.2 Calculation Time vs. Fowler Distance

The figure below indicates how much time will be required to obtain a given
Fowler distance using Fowler’s original source code. For the purposes of this
paper, the Fowler source code is “unoptimized”, as it does not contain my
optimizations and is a baseline for comparison. This graph is perhaps the most
important graph of them all, since we often want gates with a certain specific
precision.

Since there aren’t very many unique distances, there are not enough data
points to establish a clear trend. A power function appears to fit the data
somewhat closely, though. This power function predicts that the unoptimized
version of the program will take about 110 years to approximate the gate to a
distance closer than $10^{-7}$! This massive exponential expansion explains
why Fowler’s original paper had no gates with a precision better than
$10^{-4}$, since it would take at least a day for the $\frac{\pi}{6}$ gate to
compile to even that precision!

###  4.3 Time vs. Sequence Length

This metric is related to the above metric because longer sequences tend to
have better precision. However, the relationship between time and sequence
length is much clearer, as can be witnessed by the much smoother curve. While
this graph may not have as much practical significance, it is much easier to
relate this graph to the underlying implementation of the algorithm.

From sequence length 0 to 2, the line has a steep slope. This feature probably
exists because the processor cache has not warmed up yet. Between 2 and 15,
every sequence generated by the algorithm is checked against a list of unique
sequences, to see if it’s unique. This check only occurs up to a certain
sequence length: 15 in this case. After that, the algorithm speeds up very
rapidly until it reaches about a sequence length of 30. Then, the graph
becomes a clean exponential curve.

To improve performance, I will effectively need to shift this curve down,
producing longer sequences in less time.

###  4.4 Unique Sequences Per Sequence Length

This metric provides insight into the algorithm’s storage requirements. It is
clear that Fowler’s optimizations have not altered the fundamental exponential
nature of the problem. For sequences longer than about 3, the number of unique
sequences grows exponentially with the sequence length. Since I am more
worried about time rather than space, I will not mind if this curve shifts up.
However, I do need to make sure that my optimizations do not consume too much
memory.

##  5 Ways to Improve Performance

To optimize the performance, I need to:

  1. 1.

Speed up calculations such as matrix multiplication.

  2. 2.

Reduce the number of calculations required for a given gate sequence length.

There are quite a few possible approaches to approaches 1 and 2. Some of these
approaches were taken this quarter, yet others will be left for future work.

###  5.1 “Meet in the Middle” Bidirectional Search

A traditional “uni-directional” search seeks a path from a start state to a
goal state by starting from the start state and exploring all possible paths.
A bidirectional search starts searching from the goal state as well. Thus, the
search paths will “meet in the middle”: each search only has to take
$\frac{N}{2}$ steps to meet the other search. Thus, instead of taking
$O\left(a^{N}\right)$ time, the algorithm only takes
$O\left(a^{\frac{N}{2}}\right)$ time. One will need some data structure to
store the paths, but inserting into this data structure does not require
exponential time. Thus, for a given amount of time, the algorithm could
compute gate sequences that are twice as long. This approach is the most
promising, and it was implemented in software.

###  5.2 Optimized Unique Matrix Lookup

The algorithm checks to see if a matrix is unique by calculating the distance
between it and all other matrices. Since 98.5% of the application’s run time
is spent in this function, optimizing it could yield significant improvements
in performance in the first stage. However, in the second stage, no more
unique matrix checks are performed; therefore, no time will be spent in this
function. Unless the first stage lasts a long time, it may not be worth the
implementation trouble. This optimization was easy to implement since the C++
standard template library provides a red-black binary search tree.

##  6 Bidirectional Search

Searching for the correct gate is like searching through nodes in a tree: for
a given sequence of gates, the computer must choose which gate to add to the
sequence to come closer to the target gate. In the diagrams below, the arrows
represent a choice of gate, and the boxes represent matrices. When an arrow is
drawn from some box A to a box B, box B is the matrix resulting from
multiplying A by some gate matrix.

In the example shown in these figures, the existing code must go through five
levels of searching in order to reach the target gate. At each new level, the
algorithm considers adding all of the available gates to _each_ sequence
generated by the previous level. Thus, each step multiplies the number of
matrices to consider by 25. So, for a sequence of length N, there will be
$25^{N}$ operations. The “meet in the middle” figure reveals that starting the
search from the start and the goal results in the computer exploring fewer
levels. Each side would only have to explore half as many levels since the
searches meet in the middle. Instead of $25^{N}$ operations, the computer can
ideally perform $2\cdot 25^{N/2}$ operations using the MITM (meet in the
middle) algorithm.

###  6.1 The Search Index

The critical component of the MITM algorithm is the structure that allows the
paths to connect. This structure effectively creates the red arrow in the MITM
figure above, matching up left matrices with right matrices. It must be
designed carefully to ensure optimal performance of the algorithm. For a given
left matrix, it should find a minimal number of right matrices which are close
to the left matrix. Thus, the data structure needs a way to parameterize all
of the matrices stored in it, using parameters that are related to the Fowler
distance between two matrices.

The simplest approach is to choose some reference matrix M, and store the
right matrices in a tree map, using their distances from M as keys. Then, to
find right matrices that are “close” to a left matrix L, the algorithm simply
measures the distance from L to M, and performs a range query for all right
matrices that have about the same distance to M. This trick works because the
Fowler distance measure obeys the triangle inequality: if two matrices L and R
are within some distance d of each other, then the difference in their
distances to some other matrix M will not be greater than d. In the figure
below, this fact is true for all matrices inside the circle.

For my implementation, I use the target gate as the reference matrix, and I
choose d to be $10^{-10}$ less than the smallest distance found so far. Since
the left matrix must check its distance from the target gate anyway, we can
re-use the distance calculation without having to cache it. Note that it is
possible for two matrices to be far away from each other while still having
the same distance to M. Thus, the range query may return false positives,
which are shown between the red lines in the figure. The triangle inequality
property simply guarantees that the range query will not leave out potential
candidates.

###  6.2 Building the Structure

For each sequence S the algorithm generates, a corresponding matrix M is
generated. M represents the transformation that S would perform on a quantum
bit. The algorithm usually assumes that S is a prefix of the solution, meaning
that other gates will be added to the end of S to reach the target gate G.
However it’s also possible to consider S as a suffix, in which gates are added
onto the beginning of S. In this case, S would work backwards from G,
attempting to come close to the identity matrix, rather than the other way
around. If the computer knows M, it can work backwards by multiplying the
inverse of M with G to get a matrix M2. Then, prefix sequences can see if S is
their suffix by comparing their matrices to M2. If a prefix matrix is close to
M2, then it would be close to G if it were multiplied by M.

Therefore, the middle structure simply needs to store as many matrices N as
possible, with pointers to their corresponding sequences. It stores a list of
binary search trees by sequence length, so that all short sequences can be
examined before long sequences.

The middle structure only has so much room to store entries, though. Since the
number of unique sequences scales exponentially with the sequence length, the
structure can store entries up to some length L before running out of memory.
Thus, the MITM algorithm does not always cut the number of search levels in
half; instead, it subtracts L from the number of search levels required to
find a solution. This approach replaces the $O(25^{N})$ cost of exploring
sequences of length N with a $O(25^{N-L})$ cost, since a well-optimized middle
structure should not have an exponential lookup time.

###  6.3 Performing the Search

Whenever the algorithm finds a new unique sequence P, it checks the middle
structure to see if one of the suffixes S can connect it to the target gate G.
Since suffixes are searched by ascending length, the first result should be of
optimal length. The search function is given a distance parameter that
indicates the maximum tolerable Fowler distance for the match; all matrices
that are farther away are skipped. If a result is found, the search function
also returns the distance D from P’s matrix to S’s matrix, so that the
distance threshold can be reduced to $D-\epsilon$ (some small value). That
way, future searches will only return more precise matches.

One problem that I noted after obtaining my results is that the real sequence
may not be of optimal length. The Clifford group contains elements that are
composed of multiple real gates, but each Clifford group element is considered
to be one gate in this algorithm. Since every sequence alternates between
Clifford group elements and T gates, the number of real gates in the sequence
of length n returned by the algorithm is about $n/2+3(n/2)$. However, the
resulting sequence will still have an optimal real length: the Clifford group
elements are ordered such that the ones comprised of multiple real gates are
visited later by the algorithm, meaning they are added to the structure at a
later time. Thus, if the structure uses a stable sort, these longer sequences
will be considered later. I am not entirely certain that my structure does so,
however, which would be a good topic for future research.

Another potential problem is that a very good suffix may be skipped because a
“sufficient” suffix was encountered first. For speed, the MITM algorithm
returns the first suffix that is within the desired distance threshold.
Technically, if this event occurs, the improved suffix would be discovered at
the next search level, so this problem should not impact correctness. However,
that means the best result might not be returned as early as possible. One
sufficient correction would be to continue the search; it won’t impact
performance because new sequences are rarely found. This fix could be
implemented in future work.

###  6.4 Results

As the graph below shows, the “meet in the middle” (MITM) optimization
improved performance by an order of magnitude. Instead of taking about one
hour to calculate a gate sequence that is within $10^{-3}$ of the target gate,
it takes about ten minutes. The Unoptimized and MITM Width 15 lines both used
a “width” of 15, meaning that the middle structure and Fowler’s data
structures stored sequences of length 15. The actual improvement appears to
depend on the width of the middle structure: when sequences of length 30 are
stored in it, the time is cut by two orders of magnitude instead of one.
_Note: ”unoptimized” refers to Fowler’s existing algorithm without the MITM
optimization, not to a simple brute-force enumeration._

Fowler’s unoptimized algorithm also improves performance when the width is
increased, because his data structures can cache more data. Thus, it makes
sense that increasing the width to 30 from 15 results in a larger improvement
than just turning on the MITM optimization.

The memory requirements are much clearer as well: the number of unique
sequences increases exponentially with the sequence length. I omitted data for
sequences of length less than five because they adversely affect the
exponential curve fit.

Finally, I noticed that the number of sequences per unit of time was much
larger in the optimized versions than in the unoptimized versions, confirming
my hypothesis. It clearly makes sense to keep expanding the middle structure
if possible: beyond sequences of length 30, the MITM implementation with width
15 slows down relative to the implementation with width 30. However, in the
long run, the MITM optimization does not change the base of the exponential
that governs the algorithm run time: notice that all of the lines are roughly
parallel towards the right side of the graph.

I managed to approximate the $\frac{\pi}{6}$ gate to $6.8\times 10^{-5}$
precision in about 3 hours and 5 minutes. The result is 72 gates long:

$\begin{array}[]{c}HTHT(HS)THTHTHT(HS)THT(HS)T(HS)\\\
T(HS)T(HS)THTHT(HS)THTHT(HXZ)\\\ THTHTHTHTHT(HS)THTHT(HS)THT(HS)\\\
THTHTHTHT(HS)THT(HXS)T^{\dagger}\end{array}$

##  7 Change of Basis

Since the Fowler distance is phase independent, we can adjust gates to remove
their global phase. Thus, it is possible to represent a quantum gate in
$SU\left(2\right)$ by using just four real numbers. In the equation below,
$\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$ are the Pauli basis matrices.
Since they are multiplied by $i$, the basis is called the _modified Pauli
basis_.

$\displaystyle A$ $\displaystyle=$ $\displaystyle a_{0}\cdot
I+a_{1}\cdot\sigma_{x}+a_{2}\cdot\sigma_{y}+a_{3}\cdot\sigma_{z}$ (2)
$\displaystyle=$ $\displaystyle a_{0}\cdot\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)+a_{1}\cdot\left(\begin{array}[]{cc}0&i\\\
i&0\end{array}\right)+a_{2}\cdot\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)+a_{3}\cdot\left(\begin{array}[]{cc}i&0\\\
0&-i\end{array}\right)$ (11) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}a_{0}+a_{3}\cdot i&a_{2}+a_{1}\cdot
i\\\ -a_{2}+a_{1}\cdot i&a_{0}-a_{3}\cdot i\end{array}\right)$ (14)

One advantage of this new basis is that the trace distance between two gates
$A$ and $B$ is just double the dot product of their vectors $a$ and $b$:

$\displaystyle tr\left(A\cdot B^{\dagger}\right)$ $\displaystyle=$
$\displaystyle tr\left(\left(\begin{array}[]{cc}a_{0}+a_{3}\cdot
i&a_{2}+a_{1}\cdot i\\\ -a_{2}+a_{1}\cdot i&a_{0}-a_{3}\cdot
i\end{array}\right)\cdot\left(\begin{array}[]{cc}b_{0}+b_{3}\cdot
i&b_{2}+b_{1}\cdot i\\\ -b_{2}+b_{1}\cdot i&b_{0}-b_{3}\cdot
i\end{array}\right)^{\dagger}\right)$ (19) $\displaystyle=$ $\displaystyle
tr\left(\left(\begin{array}[]{cc}a_{0}+a_{3}\cdot i&a_{2}+a_{1}\cdot i\\\
-a_{2}+a_{1}\cdot i&a_{0}-a_{3}\cdot
i\end{array}\right)\cdot\left(\begin{array}[]{cc}b_{0}-b_{3}\cdot
i&-b_{2}-b_{1}\cdot i\\\ b_{2}-b_{1}\cdot i&b_{0}+b_{3}\cdot
i\end{array}\right)\right)$ (24) $\displaystyle=$ $\displaystyle
tr\left(\left(\begin{array}[]{cc}\begin{array}[]{c}\left(a_{0}+a_{3}i\right)\cdot\left(b_{0}-b_{3}i\right)+\\\
\left(a_{2}+a_{1}i\right)\cdot\left(b_{2}-b_{1}i\right)\end{array}&\ldots\\\
\ldots&\begin{array}[]{c}\left(-a_{2}+a_{1}i\right)\cdot\left(-b_{2}-b_{1}i\right)+\\\
\left(a_{0}-a_{3}i\right)\cdot\left(b_{0}+b_{3}i\right)\end{array}\end{array}\right)\right)$
(31) $\displaystyle=$
$\displaystyle\begin{array}[]{c}\left(a_{0}+a_{3}i\right)\cdot\left(b_{0}-b_{3}i\right)+\left(a_{2}+a_{1}i\right)\cdot\left(b_{2}-b_{1}i\right)+\\\
\left(-a_{2}+a_{1}i\right)\cdot\left(-b_{2}-b_{1}i\right)+\left(a_{0}-a_{3}i\right)\cdot\left(b_{0}+b_{3}i\right)\end{array}$
(34) $\displaystyle=$
$\displaystyle\begin{array}[]{c}\left(a_{0}b_{0}+a_{3}b_{3}-a_{0}b_{3}i+b_{0}a_{3}i\right)+\left(a_{2}b_{2}+a_{1}b_{1}-a_{2}b_{1}i+a_{1}b_{2}i\right)\\\
\left(a_{2}b_{2}+a_{1}b_{1}+a_{2}b_{1}i-a_{1}b_{2}i\right)+\left(a_{0}b_{0}+a_{3}b_{3}+a_{0}b_{3}i-a_{3}b_{0}i\right)\end{array}$
(37) $\displaystyle=$ $\displaystyle
2\left(a_{0}b_{0}+a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3}\right)=2\cdot a\cdot b$
(38)

This result is important because multiplication is an expensive operation in
computer calculation, relative to addition. Traditionally, calculating the
trace distance between two $2\times 2$ matrices $A$ and $B$ requires one to
obtain the diagonal elements of the product $AB$, which requires 4 complex
number multiplications. Since every complex number multiplication requires 4
real-number multiplications, 16 real multiplications must be performed in
total. If the matrices are in the modified Pauli basis, on the other hand,
only four real multiplications are required.

The other advantage is that multiplying two gates requires only 16 real
multiplications. A traditional $2\times 2$ matrix multiplication, on the other
hand, requires 8 complex number multiplications, or 32 real multiplications.

The final advantage is storage size: this new basis can be stored in half the
space that a full $2\times 2$ matrix would require.

The advantages of this basis are outlined in this table:

Task |  Regular Matrices |  Pauli Basis |  Improvement
---|---|---|---
Find trace distance |  16 real multiplies |  4 real multiplies |  4x speedup
Multiply matrices |  32 multiplies |  16 multiplies |  2x speedup
Store a matrix |  8 real numbers |  4 real numbers |  1/2 storage

##  8 Future Work

###  8.1 Using multidimensional spatial indices for the bidirectional search
middle structure

The bidirectional search index only uses one parameter to index the right
matrices. For the reference matrices $M$ which I chose, many matrices had
similar Fowler distances to $M$. Thus, while the algorithm was able to avoid
iterating over some right matrices, it still had to iterate over many matrices
that were not close to a given left matrix. In fact, only .0003% of the
matrices returned by the index were actual matches.

The modified Pauli basis offers an excellent way to parameterize the right
matrices in a spatial index:

  1. 1.

Its compact representation requires less space than a full matrix would. In
fact, since one can derive one component from any other three components, only
three components are strictly required. Space is not the only advantage;
certain spatial indices, such as k-d trees, perform better with low-
dimensional data. Hardware implementations of the algorithm also benefit from
simpler calculation circuitry.

  2. 2.

Since the trace distance is just the dot product of a left matrix vector $a$
with a right matrix vector $b$, all right matrices $b$ that are close to some
left matrix satisfy this equation:

$\displaystyle-D\leq a\cdot b\leq D$ (39)

where $D$ is some constant related to the maximum trace distance between the
two gates. Geometrically, this means all of the close right matrices are
between two parallel hyperplanes. The process of finding points between the
hyperplanes should be straightforward to optimize. Many spatial indices group
points into bounding volumes like boxes or spheres; checking to see if these
volumes are between the parallel hyperplanes is a simple process.

Libraries such as FLANN [3] provide a wide variety of spatial indices to use.

###  8.2 Map-Reduce Parallelism

The Fowler algorithm can be broken down into a cycle for each sequence length.
Each cycle is essentially a map-reduce job. During the map phase, we assign
one gate to each computer, and that computer will consider all sequences of
length $n$ which start with that gate. Once all computers have finished the
cycle, the reduce phase will merge the data structures for unique matrices, as
well as the discovered gate sequences.

There are several advantages to map-reduce parallelism: Unique sequence data
structures can be shared with all the units between cycles. Thus, all units
can benefit from each unit’s work in each subsequent calculation cycle. If you
keep track of the data structure contents after the final stage, you can
restart the algorithm from this final stage. No specialized hardware (such as
a FPGA) is required. Anyone with access to Amazon’s Elastic MapReduce service,
or a Hadoop cluster, can use a map-reduce algorithm.

Map-reduce parallelism will probably divide the algorithm’s run-time for a
given sequence length by the number of computers involved. Thus, if there are
25 computers (for 25 gates), then the algorithm ought to run up to 25 times
faster. However, since all of the computers must merge their data after each
cycle, the faster computers must wait for the slower ones. Due to the
complexity of the map-reduce setup, this method was not implemented this
quarter. However, Amazon provides a map-reduce framework that should be
straightforward to use and scale, should someone decide to adapt the program.

##  9 Related Work

A variation of the MITM algorithm was independently invented by researchers at
the Institute for Quantum Computing at the University of Waterloo [1]. This
group also seeks to find quantum circuits of optimal length implementing a
given quantum gate. There are a few key differences between their research and
the work presented here:

  1. 1.

Their work applies the algorithm to multiple-qubit gates, and does not combine
it with Fowler’s algorithm.

  2. 2.

They focus on finding _exact_ matches, rather than approximate ones.

Their future work may benefit from the approximate matching technique
discussed in this paper, as well as the brief discussion of using spatial
indices and a change of basis to accelerate matching. My research will benefit
from their more rigorous treatment of the algorithm, as well as its extension
to multiple qubits.

##  10 Summary

I considered a variety of optimizations to Fowler’s quantum compiler
algorithm. Then, I implemented the “meet in the middle” algorithm in software,
as well as a change of basis technique, and I presented the results here.
While the algorithm certainly provides a dramatic performance boost, it also
requires a lot of memory to maintain the middle index structure I introduced.
Future work involves using map-reduce parallelism and better spatial indices
to improve performance.

##  11 Acknowledgements

I performed most of this research independently, but received significant
guidance and assistance from the following individuals and organizations.
Without their involvement, this research project would not have happened!

  1. 1.

Paul Pham – the UW graduate student who suggested the research topic for this
project, and who provided essential quantum computing context and advice. I
had weekly meetings with him, and I worked with him on his pulse sequence
board two years ago. He is working on his own quantum compiler based on the
Solovay-Kitaev Theorem.

  2. 2.

Austin Fowler – a Research Fellow in Quantum Computer Science at the
University of Melbourne. He wrote the original paper describing the sequence-
skipping optimization, upon which my research is based. He also supplied the C
source code to his algorithm, so that I could test my optimizations.

  3. 3.

Aram Harrow – my faculty advisor, who came up with smart suggestions for error
accumulation analysis and calculation optimization. He also indirectly
proposed the MITM algorithm at the beginning of this research project.

## References

  * [1] Matthew Amy, Dmitry Maslov, Michele Mosca, and Martin Roetteler.  A meet-in-the-middle algorithm for fast synthesis of depth-optimal quantum circuits.
  * [2] Austin G. Fowler.  Constructing arbitrary steane code single logical qubit fault-tolerant gates.  Quantum Info. Comput., 11(9-10):867–873, September 2011.
  * [3] Marius Muja and David G. Lowe.  Fast approximate nearest neighbors with automatic algorithm configuration.  In International Conference on Computer Vision Theory and Application VISSAPP’09), pages 331–340. INSTICC Press, 2009.
  * [4] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, Cambridge, U.K., 2000.