# Derivation of Dirac Exchange Interaction Potential from Quantum Plasma
Kinetic Theory

Fernando Haas111E-mail<EMAIL_ADDRESS>Orcid: 0000-0001-8480-6877
Instituto de Física, Universidade Federal do Rio Grande do Sul, Av. Bento
Gonçalves 9500, 91501-970 Porto Alegre, RS, Brasil

###### Abstract

The Dirac exchange interaction is derived from recent quantum kinetic theory
for collisionless plasmas. For this purpose, the kinetic equation is written
in the semiclassical and long wavelength approximations. The validity of the
model for real systems is worked out, in terms of temperature and density
parameters. Within the region of applicability, the correlation potential
energy is shown to be always smaller than the exchange contribution. From the
moments of the quantum kinetic equations, macroscopic, hydrodynamic equations
are found, for an electron-ion plasma. The Dirac exchange term is explicitly
derived, in the case of a completely degenerate electron gas. These results
show, within quantum kinetic theory for charged particle systems, a new view
of the Dirac exchange interaction frequently used in density functional theory
parametrization. Finally, a simpler form of the quantum plasma exchange
kinetic theory is also found.

quantum plasmas; exchange interaction; Dirac potential.

##  I Introduction

Exchange effects have been recently included in quantum plasma kinetic theory
Zamanian ; Zamanian2 . For this purpose, the complete antisymmetry of the
N-particle density matrix was taken into account, which is inline with the
Fermi statistics of electrons, in accordance with the Hartree-Fock
approximation. For linear waves in completely degenerate plasmas Ekman , the
results are in exact agreement with older treatments following another methods
Roos -Nozieres .

There has been some discussion about the accordance between the new quantum
plasma kinetic theory and the density functional (DFT) modeling, as well as
about the adequate quantum hydrodynamics taking into account exchange effects
Ekman . Due to the analytical and even computational complexity of phase-space
models, the reduction of variables in averaged, macroscopic theories is an
attractive alternative, for both quantum and classical plasma. Exchange (and
correlation) effects have been considered e.g. in the recent review Manfredi
about quantum plasma hydrodynamics.

The purpose of the present work is to show the complete agreement between
exchange quantum kinetic plasma theory Zamanian ; Zamanian2  and the
traditional Dirac exchange potential often considered in DFT. This goes a step
further than previous treatments where the exchange (and correlation)
potential was inserted ab initio in the Schrödinger equation hh .

The article is organized as follows. In Section II, we consider the
semiclassical limit of the new exchange kinetic theory. This allows
considerable simplification of the central kinetic equation. The validity
conditions and the approximations made are described in detail. Section III
works out the resulting susceptibilities in electron-ion quantum plasmas, for
general electrostatic linear waves. With appropriate averaging methods, in
Section IV the general exchange fluid equations are derived. Section V shows
that in the particular case of a completely degenerate Fermi gas the equation
of state is exactly compatible with the Dirac exchange term. Section VI has
the conclusions, also pointing the possible generalizations.

##  II Exchange kinetic contribution in the semiclassical limit

Assuming the two-particle electron density matrix in the form of an
antisymmetric product of two one-particle density matrices, the Pauli
principle becomes assured. Following this approach, a kinetic equation for the
one-particle Wigner function $f=f({\bf x},{\bf p},t)$ can be deduced Zamanian
; Zamanian2 , neglecting spin polarization, quantum diffraction and
correlation effects, in the long wavelength approximation, reading

$\displaystyle\frac{\partial}{\partial t}\\!\\!\\!$ $\displaystyle f$
$\displaystyle\\!\\!\\!({\bf x},{\bf p},t)+\\!\\!\frac{{\bf p}}{m}\cdot\nabla
f({\bf x},{\bf p},t)+e\nabla\phi({\bf x},t)\cdot\frac{\partial}{\partial{\bf
p}}f({\bf x},{\bf p},t)$ $\displaystyle=$
$\displaystyle\frac{\hbar}{2}\frac{\partial}{\partial
p_{i}}\int\\!\\!d^{3}\rho\,d^{3}q\,e^{-i{\bm{\rho}}\cdot{\bf
q}}\,\frac{\partial V(\bm{\rho})}{\partial\rho_{i}}\,f\left({\bf
x}-\frac{\hbar\bm{\rho}}{2},{\bf p}+\frac{\bf q}{2},t\right)f\left({\bf
x}-\frac{\hbar\bm{\rho}}{2},{\bf p}-\frac{\bf q}{2},t\right)$ $\displaystyle-$
$\displaystyle\\!\\!\frac{i\hbar^{2}}{8}\frac{\partial}{\partial
p_{i}}\frac{\partial}{\partial
p_{j}}\\!\int\\!\\!d^{3}\rho\,d^{3}q\,e^{-i{\bm{\rho}}\cdot{\bf
q}}\,\frac{\partial V(\bm{\rho})}{\partial\rho_{i}}\,\left[f\left({\bf
x}\\!-\\!\frac{\hbar\bm{\rho}}{2},{\bf p}\\!-\\!\frac{\bf
q}{2},t\right)\\!\left(\\!\frac{\overleftarrow{\partial}}{\partial
x_{j}}-\frac{\overrightarrow{\partial}}{\partial
x_{j}}\\!\right)\\!f\left({\bf x}\\!-\\!\frac{\hbar\bm{\rho}}{2},{\bf
p}\\!+\\!\frac{\bf q}{2},t\right)\right].$

Here $\phi=\phi({\bf x},t)$ is the scalar potential and $V=V({\bf
x})=e^{2}/(4\pi\varepsilon_{0}|{\bf x}|)$ is the electron-electron Coulomb
potential. The summation convention holds and the remaining symbols have their
usual meaning while the arrows in some partial derivatives indicate the sense
of operation. By construction Eq. (II) takes into account the bare exchange
effects, represented by the terms on the right hand side, and can be termed
the exchange kinetic equation, which is our starting point. Details are
discussed in Zamanian ; Zamanian2 .

To have the exchange kinetic equation in a semiclassical approximation, we
have to make a formal series expansion in powers of $\hbar$. In the long
wavelength approximation, the calculation (shown in the Appendix) gives the
semiclassical exchange kinetic equation

$\displaystyle\frac{\partial}{\partial t}\\!\\!\\!$ $\displaystyle f$
$\displaystyle\\!\\!\\!({\bf x},{\bf p},t)+\\!\\!\frac{{\bf p}}{m}\cdot\nabla
f({\bf x},{\bf p},t)+e\nabla\phi({\bf x},t)\cdot\frac{\partial}{\partial{\bf
p}}f({\bf x},{\bf p},t)$ (2) $\displaystyle=$
$\displaystyle\frac{e^{2}\hbar^{2}}{2\varepsilon_{0}}\frac{\partial}{\partial
p_{i}}\frac{\partial}{\partial
x_{j}}\int\frac{d^{3}q}{q^{2}}\left(\delta_{ij}-\frac{2q_{i}q_{j}}{q^{2}}\right)f({\bf
x},{\bf p}+{\bf q},t)f({\bf x},{\bf p}-{\bf q},t)\,,$

which is considerably simpler than Eq. (II). The last integral on the right
hand side of Eq. (II) is of higher order and does not contribute in the
semiclassical approximation.

The formal series expansion in powers of $\hbar$ implicitly implies a power
series on a dimensionless quantity $\hbar/(m\,v_{0}\,l_{0})\ll 1$, where
$v_{0}$ is a natural velocity scale and $l_{0}$ is a natural length scale. The
same approximation justifies the neglect of quantum diffraction, also known as
quantum recoil, which was already assumed for the derivation Zamanian ;
Zamanian2  of the basic kinetic equation (II). The semiclassical approximation
keeps the first order exchange effects. It is a welcome avenue, to analyze the
implications of such approximation.

Since exchange effects are more stringent for dense plasmas, we focus on
degenerate plasmas. For a degenerate plasma, a natural velocity scale is the
Fermi velocity $v_{F}$, so that $v_{0}=v_{F}$. On the other hand, a
characteristic length $l_{0}$ could be given by the electron Thomas-Fermi
length $\lambda_{F}=v_{F}/\omega_{p}$ for stationary structures book , where
$\omega_{p}$ is the plasma frequency, or by $\lambda=2\pi/k$ for linear wave
propagation with wavenumber $k$. For the sake of definiteness, setting
$l_{0}=\lambda_{F}$ and taking into account
$v_{F}=\hbar(3\pi^{2}n_{0})^{1/3}/m,\,\,\omega_{p}=[n_{0}e^{2}/(m\varepsilon_{0})]^{1/2}$,
where $n_{0}$ is the equilibrium electrons number density, $m$ is the electron
mass, $e$ is the elementary charge and $\varepsilon_{0}$ is the vacuum
permittivity, the semiclassical approximation traduces into
$a_{0}\,n_{0}^{1/3}\gg 4\pi/(3\pi^{2})^{4/3}=0.14$, where
$a_{0}=4\pi\varepsilon_{0}\hbar^{2}/(me^{2})$ is the Bohr radius. Hence, for
stationary structures in degenerate plasma the semiclassical approximation
applies for $n_{0}\gg 1.7\times 10^{28}\,{\rm m}^{-3}$, a condition on the
number density only.

The remaining validity conditions are as follows. Neglect of correlations
imply that the average Coulomb energy
$E_{C}=e^{2}n_{0}^{1/3}/(4\pi\varepsilon_{0})$ is much smaller than the
average kinetic energy, which is the Fermi energy $E_{F}=mv_{F}^{2}/2$ in an
order of magnitude estimate. Defining the coupling parameter
$g_{c}=E_{C}/E_{F}$ we then find the necessary condition

$g_{c}=\frac{E_{C}}{E_{F}}=\frac{2}{(3\pi^{2})^{2/3}a_{0}\,n_{0}^{1/3}}\ll
1\,,$ (3)

which at the end is again a function of the number density only. Equation (3)
implies $n_{0}\gg 6.1\times 10^{28}\,{\rm m}^{-3}$, which is very similar to
the semiclassical condition. Moreover, it is often overlooked that the
application of a non-relativistic model is possible only when relativistic
effects are negligible. To avoid a relativistic mass increase, one needs at
least $v_{F}/c\ll 1$, where $c$ is the speed of light. This amounts to
$n_{0}\lambda_{C}^{3}\ll 1/(3\pi^{2})$, where $\lambda_{C}=\hbar/(mc)$ is the
Compton length. Therefore, $n_{0}\ll 5.9\times 10^{35}\,{\rm m}^{-3}$,
excluding very dense plasmas deserving a relativistic treatment. Finally, the
degeneracy condition is $T\ll T_{F}$, where $T$ is the thermodynamic
temperature, $T_{F}=E_{F}/\kappa_{B}$ is the Fermi temperature and
$\kappa_{B}$ is the Boltzmann condition. All in all, the region for which the
present modeling is applicable is shown in the filled area in the density-
temperature diagram in Fig. 1, in a logarithmic scale.

Figure 1: Density-temperature diagram for which the present modeling is
applicable, in a logarithmic scale, where $n_{0}$ is the equilibrium number
density (measured in ${\rm m}^{-3}$) and $T$ is the thermodynamic temperature
in $K$. The left vertical line indicates the minimal number density for which
correlation effects are negligible, also justifying the semiclassical
approximation. The right vertical line indicates the maximal number density to
avoid relativistic effects. Degenerate plasmas are below the straight line
$T=T_{F}$. The modeling is valid in the filled area.

On the other hand, concerning wave propagation and setting $l_{0}=\lambda$,
from $\hbar/(mv_{F}\lambda)\ll 1$ one gets $\lambda
n_{0}^{1/3}>>1/(3\pi^{2})^{1/3}$ which means the wavelength should be much
larger than the mean inter-particle distance. Therefore the semiclassical
approximation is equivalent to a continuous media (fluid) assumption, a long
wavelength condition in this case.

It is interesting to look into the neglect of correlations from the point of
view of density-functional theory. The DFT exchange-correlation potential
$V_{XC}$ is given Hedin ; Brey  by

$V_{XC}=V_{X}+V_{C}=g_{D}\left(\frac{n}{n_{0}}\right)^{1/3}\left[1+\frac{0.034}{a_{0}n^{1/3}}\,\ln\left(1+18.37a_{0}n^{1/3}\right)\right]\,,$
(4)

where

$V_{X}=g_{D}\left(\frac{n}{n_{0}}\right)^{1/3}\,,\quad
g_{D}=0.985\,\frac{(3\pi^{2})^{2/3}}{4\pi}\frac{\hbar^{2}\omega_{p}^{2}}{mv_{F}^{2}}$
(5)

is the Dirac potential form of the exchange potential Dirac  and

$V_{C}=V_{X}\times\left[\frac{0.034}{a_{0}n^{1/3}}\,\ln\left(1+18.37a_{0}n^{1/3}\right)\right]$
(6)

is the correlation potential. Evaluating at the equilibrium number density and
using Eq. (3), one has

$\frac{V_{C}}{V_{X}}=0.16\,g_{c}\,\ln\left(1+\frac{3.84}{g_{c}}\right)\,,$ (7)

shown in Fig. 2. As expected, the relative importance of the correlation
effects (similar to collisionall effects) goes to zero as the coupling
parameter $g_{c}\ll 1$. On the opposite limit, when $g_{c}\gg 1$ one has
$V_{C}/V_{X}=0.63$, which is a considerable value but still such that
$V_{C}<V_{X}$, which is thus always true at least from the point of view of
the effective DFT exchange-correlation potential.

Figure 2: Ratio between effective correlation and exchange potentials as a
function of the coupling parameter $g_{c}$, from Eq. (7). The horizontal
straight line shows the limiting value $V_{C}/V_{X}=0.63$ as $g_{C}\gg 1$.

##  III Susceptibilities in an Electron-Ion Plasma

In an electron-ion plasma we have to consider the ions probability
distribution function $f_{i}=f_{i}({\bf x},{\bf p},t)$. It follows a classical
Vlasov equation

$\frac{\partial}{\partial t}f_{i}({\bf x},{\bf p},t)+\frac{{\bf
p}}{M}\cdot\nabla f_{i}({\bf x},{\bf p},t)-e\nabla\phi({\bf
x},t)\cdot\frac{\partial}{\partial{\bf p}}f_{i}({\bf x},{\bf p},t)=0\,,$ (8)

due to a presumably larger ion mass $M\gg m$. For simplicity ions have been
assumed to be single ionized. To close the system we have Poisson’s equation

$\nabla^{2}\phi=\frac{e}{\varepsilon_{0}}\,\left(\int d^{3}p\,f({\bf x},{\bf
p},t)-\int d^{3}p\,f_{i}({\bf x},{\bf p},t)\right)\,.$ (9)

For small amplitude waves it is assumed plane wave perturbations as

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{0}({\bf p})+\delta f({\bf
p})\,\exp[i({\bf k}\cdot{\bf x}-\omega t)]\,,\quad f_{i}=f_{0(i)}({\bf
p})+\delta f_{i}({\bf p})\,\exp[i({\bf k}\cdot{\bf x}-\omega t)]\,,$
$\displaystyle\phi$ $\displaystyle=$ $\displaystyle\delta\phi\,\exp[i({\bf
k}\cdot{\bf x}-\omega t)]\,,$ (10)

where $\delta$ denotes first order quantities. Linearizing yields the
dispersion relation

$1+\chi_{e}+\chi_{i}=0\,,$ (11)

where the electron-ion susceptibilities $\chi_{e,i}$ are

$\chi_{e}=\frac{e}{\varepsilon_{0}k^{2}}\int d^{3}p\,\frac{\delta
f}{\delta\phi}\,,\quad\chi_{i}=-\,\frac{e}{\varepsilon_{0}k^{2}}\int
d^{3}p\,\frac{\delta f_{i}}{\delta\phi}\,.$ (12)

Linearizing the ions kinetic equation (8) and assuming a cold ions equilibrium
$f_{i(0)}=n_{0}\,\delta({\bf p})$ we get
$\chi_{i}=-\omega_{pi}^{2}/\omega^{2}$, where
$\omega_{pi}=[n_{0}e^{2}/(M\varepsilon_{0})]^{1/2}$. For electrons we
linearize the exchange kinetic equation (2) to get

$\displaystyle\left(-\omega+\frac{{\bf k}\cdot{\bf p}}{m}\right)\,\delta
f({\bf p})+e\,\delta\phi\,{\bf k}\cdot\frac{\partial f_{0}({\bf
p})}{\partial{\bf p}}$ (13)
$\displaystyle=\frac{e^{2}\hbar^{2}}{2\varepsilon_{0}}\,\frac{\partial}{\partial
p_{i}}\int\frac{d^{3}q}{q^{2}}\left(k_{i}-2\,{\bf k}\cdot{\bf
q}\,\frac{q_{i}}{q^{2}}\right)\Bigl{(}f_{0}({\bf p}+{\bf q})\,\delta f({\bf
p}-{\bf q})+f_{0}({\bf p}-{\bf q})\,\delta f({\bf p}+{\bf q})\Bigr{)}\,.$

At this point we remember that under the semiclassical condition the entire
right-hand side of Eq. (13) is itself a correction. Hence it is legitimate to
insert the classical expression

$\delta f({\bf p})=\frac{e\,\delta\phi\,{\bf k}\cdot\partial f_{0}({\bf
p})/\partial{\bf p}}{\omega-{\bf k}\cdot{\bf p}/m}$ (14)

into the exchange term of Eq. (13). In this way the electrons susceptibility
decomposes into classical $\chi_{e}^{C}$ and exchange $\chi_{e}^{X}$
contributions such that $\chi_{e}=\chi_{e}^{C}+\chi_{e}^{X}$, where

$\displaystyle\chi_{e}^{C}$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{\varepsilon_{0}k^{2}}\int d^{3}p\,\frac{{\bf
k}\cdot\partial f_{0}({\bf p})/\partial{\bf p}}{\omega-{\bf k}\cdot{\bf
p}/m}\,,$ (15) $\displaystyle\chi_{e}^{X}$ $\displaystyle=$
$\displaystyle-\,\frac{1}{2}\left(\frac{e^{2}\hbar}{\varepsilon_{0}k}\right)^{2}\int\frac{d^{3}p}{\omega-{\bf
k}\cdot{\bf p}/m}\,\frac{\partial}{\partial
p_{i}}\int\frac{d^{3}q}{q^{2}}\times$ (16) $\displaystyle\times$
$\displaystyle\left(k_{i}-2\,{\bf k}\cdot{\bf
q}\,\frac{q_{i}}{q^{2}}\right)\left(\frac{f_{0}(+)\,{\bf k}\cdot\partial
f_{0}(-)/\partial{\bf p}}{\omega-{\bf k}\cdot({\bf p}-{\bf
q})/m}+\frac{f_{0}(-)\,{\bf k}\cdot\partial f_{0}(+)/\partial{\bf
p}}{\omega-{\bf k}\cdot({\bf p}+{\bf q})/m}\right)\,,$

in terms of the shorthand $f_{0}(\pm)=f_{0}({\bf p}\pm{\bf q})$.

It is instructive to keep the full kinetic equation (II) and perform again the
previous operations, to find the same $\chi_{e}^{C}$ as in Eq. (15) and

$\displaystyle\chi_{e}^{X}$ $\displaystyle=$
$\displaystyle-\,\frac{2\,\hbar\,e^{4}}{\varepsilon_{0}^{2}k^{2}}\int\frac{d^{3}p}{\omega-{\bf
k}\cdot{\bf p}/m}\,\frac{\partial}{\partial
p_{i}}\int\frac{d^{3}q\,(q_{i}+\hbar k_{i}/4)}{|{\bf q}+\hbar{\bf
k}/4|^{2}}\times$ (17) $\displaystyle\times$
$\displaystyle\left(\frac{f_{0}(+)\,{\bf k}\cdot\partial f_{0}(-)/\partial{\bf
p}}{\omega-{\bf k}\cdot({\bf p}-{\bf q})/m}+\frac{f_{0}(-)\,{\bf
k}\cdot\partial f_{0}(+)/\partial{\bf p}}{\omega-{\bf k}\cdot({\bf p}+{\bf
q})/m}\right)$

instead of Eq. (16). However, expanding taking into account the long
wavelength condition $\hbar k\ll mv_{0}$ as discussed in Section II it is easy
to regain Eq. (16) from Eq. (17), as verified from parity properties. In
particular it is immediate to verify that these parity properties show that
the first non-zero contribution is proportional to $\hbar^{2}$.

##  IV Fluid Equations

Macroscopic equations can be derived from the moments of the simplified
exchange kinetic equation (2). In this way it is even possible to derive an
exchange potential $V_{X}$ from first principles. In the case of completely
degenerate electrons it turns out that $V_{X}$ is the Dirac exchange
interaction potential, as will be shown next.

As usual h2021  the number density $n=n({\bf x},t)$, the velocity field ${\bf
u}={\bf u}({\bf x},t)$ and the pressure dyad ${\bf P}={\bf P}({\bf x},t)$ are
defined in terms of the moments

$\displaystyle n$ $\displaystyle=$ $\displaystyle\int d^{3}p\,f\,,$ (18)
$\displaystyle mn{\bf u}$ $\displaystyle=$ $\displaystyle\int d^{3}p\,{\bf
p}\,f\,,$ (19) $\displaystyle{\bf P}$ $\displaystyle=$
$\displaystyle\frac{1}{m}\int d^{3}p\,{\bf p}\otimes{\bf p}\,f-mn{\bf
u}\otimes{\bf u}\,.$ (20)

Higher order moments could be easily implemented as well.

From Eq. (2) and with appropriate boundary conditions in momentum space we
find the continuity equation

$\frac{\partial n}{\partial t}+\nabla\cdot(n{\bf u})=0\,,$ (21)

and the momentum balance equation

$\displaystyle m\,\left(\frac{\partial}{\partial t}+{\bf
u}\cdot\nabla\right)u_{i}$ $\displaystyle=$
$\displaystyle-\,\frac{1}{n}\frac{\partial P_{ij}}{\partial
x_{j}}+e\,\frac{\partial\phi}{\partial x_{i}}$
$\displaystyle-\,\frac{e^{2}\hbar^{2}}{2\,\varepsilon_{0}n}\frac{\partial}{\partial
x_{j}}\int\frac{d^{3}p\,d^{3}q}{q^{2}}\left(\delta_{ij}-\frac{2\,q_{i}q_{j}}{q^{2}}\right)\Bigl{(}f({\bf
x},{\bf p}+{\bf q})f({\bf x},{\bf p}-{\bf q})\Bigr{)}\,,$

expressed in component-wise manner, where the explicit time-dependence was
omitted in $f$, for brevity. The momentum balance equation is valid for any
local $f({\bf x},{\bf p})$.

Equation (IV) can be written as

$m\,\left(\frac{\partial}{\partial t}+{\bf u}\cdot\nabla\right){\bf
u}=-\,\frac{\nabla\cdot{\bf P}}{n}-\,\frac{\nabla\cdot{\bf
P}^{X}}{n}+e\,\nabla\phi\,,$ (23)

in terms of the exchange pressure dyad ${\bf P}^{X}$ defined by

$P_{ij}^{X}=\,\frac{e^{2}\hbar^{2}}{2\,\varepsilon_{0}}\int\frac{d^{3}p\,d^{3}q}{q^{2}}\left(\delta_{ij}-\frac{2\,q_{i}q_{j}}{q^{2}}\right)\Bigl{(}f({\bf
x},{\bf p}+{\bf q})f({\bf x},{\bf p}-{\bf q})\Bigr{)}\,,$ (24)

in component-wise manner. The simplicity of Eq. (23) in comparison with Eq.
(7) of Ref. h2021  comes from starting from Eq. (2) instead of the full
exchange kinetic equation (II), which involves no loss of generality as long
as the validity conditions discussed in Section II remain valid.

For isotropic in momentum distributions, namely if $f=f({\bf x},p),p=|{\bf
p}|$, it follows that $P_{ij}=p\,\delta_{ij}$ and
$P_{ij}^{X}=p_{X}\delta_{ij}$, where $p$ is the scalar pressure $p=(1/3){\rm
Tr}P_{ij}$ (${\rm Tr}$ denotes the trace) and where $p_{X}$ is the exchange
scalar pressure,

$p_{X}=\frac{1}{3}{\rm
Tr}P_{ij}^{X}=\frac{e^{2}\hbar^{2}}{6\,\varepsilon_{0}}\int\frac{d^{3}p\,d^{3}q}{q^{2}}\Bigl{(}f({\bf
x},{\bf p}+{\bf q})f({\bf x},{\bf p}-{\bf q})\Bigr{)}\,.$ (25)

Similarly, for ions one can define

$\displaystyle n_{i}$ $\displaystyle=$ $\displaystyle\int d^{3}p\,f_{i}\,,$
(26) $\displaystyle Mn_{i}{\bf u}_{i}$ $\displaystyle=$ $\displaystyle\int
d^{3}p\,{\bf p}\,f_{i}\,,$ (27) $\displaystyle{\bf P}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{M}\int d^{3}p\,{\bf p}\otimes{\bf p}\,f_{i}-Mn_{i}{\bf
u}_{i}\otimes{\bf u}_{i}\,.$ (28)

The moments of the Vlasov equation (8) yield

$\displaystyle\frac{\partial n_{i}}{\partial t}+\nabla\cdot(n_{i}{\bf u}_{i})$
$\displaystyle=$ $\displaystyle 0\,,$ (29) $\displaystyle
M\,\left(\frac{\partial}{\partial t}+{\bf u}_{i}\cdot\nabla\right){\bf u}_{i}$
$\displaystyle=$ $\displaystyle-\,\frac{\nabla\cdot{\bf
P}_{i}}{n_{i}}-e\,\nabla\phi\,.$ (30)

Finally one has Poisson’s equation

$\nabla^{2}\phi=\frac{e}{\varepsilon_{0}}\,(n-n_{i})\,.$ (31)

To have closure of the system, one needs to express all pressure dyads in
terms of lower order moments. This can be achieved assuming a local
equilibrium distribution function, as shown in the next Section.

##  V Closure for Completely Degenerate Electrons

For a completely degenerate electron gas, one has

$f({\bf x},{\bf p})=A\,\theta(\tilde{p}_{F}-|{\bf p}-m{\bf u}|)\,,\quad
A=\frac{3\,n}{4\pi\tilde{p}_{F}^{3}}\,,$ (32)

where $\theta$ is the Heaviside step function of the indicated argument and
where $\tilde{p}_{F}=\hbar\,(3\pi^{2}n)^{1/3}$ is the local Fermi momentum
depending on the local number density $n$. Equation (32) represents a local,
velocity displaced zero temperature Fermi-Dirac distribution.

In the reference frame of the electrons fluid, the quasi-equilibrium (32) is
isotropic in momentum space. It is immediate to evaluate the scalar pressure

$p=\frac{2}{5}\,n_{0}E_{F}\left(\frac{n}{n_{0}}\right)^{5/3}\,,\quad
E_{F}=\frac{p_{F}^{2}}{2\,m}\,,\quad p_{F}=\hbar\,(3\pi^{2}n_{0})^{1/3}\,,$
(33)

which is the equation of state for the completely degenerate electron gas
where $E_{F},p_{F}$ are the Fermi energy and momentum, and $n_{0}$ is the
equilibrium number density.

The derivation of the equation of state for the exchange scalar pressure (25)
involves the evaluation of the volume intersection between two spheres of
identical radius $\tilde{p}_{F}$, centered at ${\bf p}=\pm{\bf q}$, which is

$\int d^{3}p\,\theta(\tilde{p}_{F}-|{\bf p}+{\bf
q}|)\,\theta(\tilde{p}_{F}-|{\bf p}-{\bf
q}|)=\frac{4\pi}{3}\,(\tilde{p}_{F}+\frac{q}{2})\,(\tilde{p}_{F}-q)^{2}\,\theta(\tilde{p}_{F}-q)\,,$
(34)

as found from elementary calculus Kern . The remaining integration in Eq. (25)
gives

$p_{X}=\frac{3}{16}\,\frac{e^{2}}{\varepsilon_{0}}\,\frac{n^{4/3}}{(3\pi^{2})^{2/3}}\,,$
(35)

a barotropic equation of state. Remarkably, the contributions from the
exchange pressure from Eq. (35) and from the Dirac exchange potential are
entirely equivalent, namely

$\nabla V_{X}=\frac{\nabla p_{X}}{n}\,,$ (36)

as shown from elementary algebra. This is the main result of this work.
Finally, for ultra-cold ions obviously one has ${\bf P}_{i}=0$.

##  VI Conclusions

A complete derivation of the Dirac exchange potential has been made, from
newly introduced quantum plasma models taking into account the antisymmetry of
the two-particle Wigner distribution function. The derivation from the
exchange kinetic theory is quite different from the original derivation by
Dirac, which is based on the Thomas-Fermi atomic model Dirac .

The simpler kinetic equation (2) is more amenable to nonlinear and numerical
analysis, in comparison with Eq. (II), considered in the literature. Even more
simplicity appears in the hydrodynamic modeling, which however needs an
equation of state for closure. A local completely degenerate Fermi equilibrium
shows a pressure term entirely equivalent to Dirac’s exchange potential.
Although evident from the start, it worth to mention that a classical, locally
maxwellian equilibrium would gives something else h2021  than the Dirac
expression.

The approximations made, restrict the results to collisionless, semiclassical
and completely degenerate plasmas. The domain of validity include real, dense
fully degenerate and cold plasma systems, as identified in Fig. 1.
Nevertheless, the procedure can be in principle extended to cover a finite-
temperature electron gas, or a systematic derivation of correlation
interaction potentials. The later would need a first principle kinetic
equation including correlation besides exchange effects. However, at least in
terms of the usual DFT parametrization of the exchange-correlation potential,
the correlation contribution is shown to be a small correction in comparison
with the exchange contribution, in the case of ideal plasmas. Finally,
relativistic effects would be a further ingredient, in even more general
theories.

##  Appendix A Derivation of the semiclassical exchange kinetic equation

Denoting the right-hand side of Eq. (II) by $I$ and expanding in a formal
power series of $\hbar$, it becomes

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\,\frac{e^{2}\hbar}{8\pi\varepsilon_{0}}\frac{\partial}{\partial
p_{i}}\int
d^{3}\rho\,d^{3}q\,\frac{\rho_{i}}{\rho^{3}}\,e^{-i\bm{\rho}\cdot{\bf
q}}\,f(+)f(-)$ (37) $\displaystyle+$
$\displaystyle\frac{e^{2}\hbar^{2}}{16\pi\varepsilon_{0}}\frac{\partial}{\partial
p_{i}}\frac{\partial}{\partial x_{j}}\int
d^{3}\rho\,d^{3}q\,\frac{\rho_{i}\rho_{j}}{\rho^{3}}\,e^{-i\bm{\rho}\cdot{\bf
q}}\,f(+)f(-)$ $\displaystyle+$
$\displaystyle\frac{e^{2}\hbar^{3}}{64\pi\varepsilon_{0}}\frac{\partial}{\partial
p_{i}}\frac{\partial}{\partial p_{j}}\int
d^{3}\rho\,d^{3}q\,\frac{\rho_{i}\rho_{k}}{\rho^{3}}\,e^{-i\bm{\rho}\cdot{\bf
q}}\times\Bigl{[}f(+)\,\frac{\partial^{2}f(-)}{\partial x_{j}\partial
x_{k}}-f(-)\,\frac{\partial^{2}f(+)}{\partial x_{j}\partial x_{k}}$
$\displaystyle-$ $\displaystyle\frac{\partial f(+)}{\partial
x_{j}}\frac{\partial f(-)}{\partial x_{k}}+\frac{\partial f(-)}{\partial
x_{j}}\frac{\partial f(+)}{\partial x_{k}}\Bigr{]}+\,{\cal O}(\hbar^{4})\,,$

summation convention implied, where for simplicity the time-dependence was
omitted since it is not relevant for this discussion and where now the
shorthand

$f(\pm)=f\left({\bf x},{\bf p}\pm\frac{\bf q}{2}\right)\,.$ (38)

is used.

It happens that the terms proportional to $\hbar$ and $\hbar^{3}$ vanishes in
Eq. (37) due to parity properties, as can be verified by means of the
simultaneous change of variables ${\bf q}\rightarrow-{\bf
q},\bm{\rho}\rightarrow-\bm{\rho}$. For the surviving term proportional to
$\hbar^{2}$ we use

$\int d^{3}\rho\,\frac{\rho_{i}\rho_{j}}{\rho^{3}}\,e^{-i\bm{\rho}\cdot{\bf
q}}=i\,\frac{\partial}{\partial q_{i}}\,\int
d^{3}\rho\,\frac{\rho_{j}}{\rho^{3}}\,e^{-i\bm{\rho}\cdot{\bf
q}}=\frac{4\pi}{q^{2}}\left(\delta_{ij}-\frac{2q_{i}q_{j}}{q^{2}}\right)\,,$
(39)

as can be verified using a partial integration and the Fourier transform of
the Coulomb potential,

$\int\frac{d^{3}\rho}{\rho}\,e^{-i\bm{\rho}\cdot{\bf
q}}=\frac{4\pi}{q^{2}}\,.$ (40)

Equation (37) becomes

$I=\frac{e^{2}\hbar^{2}}{4\,\varepsilon_{0}}\frac{\partial}{\partial
p_{i}}\frac{\partial}{\partial
x_{j}}\int\frac{d^{3}q}{q^{2}}\,\left(\delta_{ij}-\frac{2q_{i}q_{j}}{q^{2}}\right)f(+)f(-)+\,{\cal
O}(\hbar^{4})\,.$ (41)

A few more simple calculations finally give Eq. (2). It is worth mentioning
that expanding to higher orders the terms proportional to $\hbar^{\nu}$ where
$\nu$ is an odd integer vanishes due to parity properties. Hence, actually one
has a power series on $\hbar^{2}$, as expected from quantum perturbation
theory in general.

###### Acknowledgements.

The author acknowledges the support by Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq) and the help of Jackson Galvão for one figure.
Data Availability Statement: the data that support the findings of this study
are available from the corresponding author upon reasonable request.

## References

  * (1) J. Zamanian, M. Marklund and G. Brodin, Phys. Rev. E 88, 063105 (2013).
  * (2) J. Zamanian, M. Marklund and G. Brodin, Eur. Phys. J. D 69, 25 (2015).
  * (3) R. Ekman, J. Zamanian and G. Brodin, Phys. Rev. E 92, 013104 (2015).
  * (4) O. von Roos and J. S. Zmuidzinas, Phys. Rev. 121, 941 (1961).
  * (5) P. Burt and H. D. Wahlquist, Phys. Rev. 125, 1785 (1962).
  * (6) H. Kanazawa, S. Misawa, and K. Fujita, Progr. Theoret. Phys. 23, 426 (1960).
  * (7) P. Nozieres and D. Pines, Phys. Rev. 111, 442 (1958).
  * (8) G. Manfredi, P. A. Hervieux and J. Hurst, Rev. Mod. Plasma Phys. 5, 7 (2021).
  * (9) F. Haas, Plasma Phys. Control. Fusion 61, 044001 (2019).
  * (10) F. Haas, Quantum Plasmas: an Hydrodynamic Approach. Springer, New York, 2011.
  * (11) L. Hedin and B. I. Lundqvist, J. Phys. C 4, 2064 (1971).
  * (12) L. Brey, J. Dempsey, N. F. Johnson and B. I. Halperin, Phys. Rev. B 42, 1240 (1990).
  * (13) P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930).
  * (14) F. Haas, Contrib. Plasma Phys. 62, e202100046 (2021).
  * (15) W. F. Kern and J. R. Bland, Solid Mensuration with Proofs, 2nd ed. Wiley, New York, 1948.


# Fourier Quantum Process Tomography

Francesco Di Colandrea<EMAIL_ADDRESS>Nexus for Quantum
Technologies, University of Ottawa, K1N 5N6, Ottawa, ON, Canada  Nazanin
Dehghan  Nexus for Quantum Technologies, University of Ottawa, K1N 5N6,
Ottawa, ON, Canada  Alessio D’Errico  Nexus for Quantum Technologies,
University of Ottawa, K1N 5N6, Ottawa, ON, Canada  Ebrahim Karimi  Nexus for
Quantum Technologies, University of Ottawa, K1N 5N6, Ottawa, ON, Canada

###### Abstract

The characterization of a quantum device is a crucial step in the development
of quantum experiments. This is accomplished via Quantum Process Tomography,
which combines the outcomes of different projective measurements to deliver a
possible reconstruction of the underlying process. The tomography is typically
performed by processing an overcomplete set of measurements and extracting the
process matrix from maximum-likelihood estimation. Here, we introduce a new
technique, referred to as Fourier Quantum Process Tomography, which requires a
reduced number of measurements, and benchmark its performance against the
standard maximum-likelihood approach. Fourier Quantum Process Tomography is
based on measuring probability distributions in two conjugate spaces for
different state preparations and projections. Exploiting the concept of phase
retrieval, our scheme achieves a complete and robust characterization of the
setup by processing a near-minimal set of measurements. We experimentally test
the technique on different space-dependent polarization transformations,
reporting average fidelities higher than 90% and significant computational
advantage.

August 28, 2024

##  I Introduction

The functionalities of a _black-box_ quantum device can be assessed via
Quantum Process Tomography (QPT) techniques. These techniques prescribe a set
of experimental measurements to identify the unknown parameters of the
underlying process matrix [1]. QPT is routinely performed across various
quantum architectures, ranging from nuclear magnetic resonances [2] to cold
atoms [3], trapped ions [4, 5], superconducting circuits [6, 7] and photonic
setups [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18].

In principle, one could extract the analytical relations between the operator
parameters and the outcomes of suitable projective measurements [19]. However,
this proves to be often incompatible with realistic experimental noise,
typically yielding nonphysical reconstructions. This inconvenience can be
overcome by formulating the process tomography as an optimization problem, as
first proposed for the tomography of quantum states [20].

In this framework, the most elementary scenario is the characterization of an
SU(2) gate $\hat{U}$ acting on a two-level quantum system (qubit).
Polarization of photons provides a natural way of encoding qubits, with
$\hat{U}$ implemented via one or multiple birefringent waveplates.
Accordingly, the characterization of devices acting on light polarization can
be accomplished via QPT [21].

Here, we address the more challenging scenario of characterizing optical SU(2)
gates that are dependent on some $d$-dimensional degree of freedom, hereafter
referred to as lattice. We introduce a new technique, named Fourier Quantum
Process Tomography (FQPT), that allows retrieving all the parameters of the
unknown transformation by processing only three sets of projective
measurements collected in 2 conjugate planes. This method applies to
SU$(2\times d)$ transformations, which can be decomposed in a $2\times 2$
block-diagonal form.

FQPT is validated experimentally on complex polarization transformations
realized via liquid-crystal metasurfaces (LCMSs) patterned with high spatial
frequencies [22], and its performance is compared with a standard maximum-
likelihood (ML) approach. In this experiment, the measurements can be
conveniently chosen to be performed in two conjugate planes, namely the near
and far field, wherein the light distributions are directly connected via a
Fourier transform. If the near field is associated with an intermediate plane
or two intermediate planes are selected, then the Fourier transform is
replaced by a paraxial Fresnel propagator [23]. FQPT can also be implemented
in other platforms. For instance, integrated photonic technologies [24] can
support additional chips specifically implementing the Quantum Fourier
Transform (QFT) algorithm [25, 26, 27]. At the same time, SU(2) operations
could be implemented either in the polarization [28, 29, 30] or path encoding.
In the latter case, the waveguide array would simulate a composite lattice.
Similar schemes have also been reported in several non-photonic platforms [31,
32, 33].

##  II Theory

A qubit rotation of an angle 2$E$ around the axis
${\bm{n}=(n_{1},n_{2},n_{3})}$, with $0\leq E<\pi$ and
$\absolutevalue{\bm{n}}=1$, is described by an SU(2) operator

$\hat{U}=e^{-iE\bm{n}\cdot\boldsymbol{\sigma}}=\cos(E)\sigma_{0}-i\sin(E)(\bm{n}\cdot\bm{\sigma}),$
(1)

where $\sigma_{0}$ is the 2$\times$2 identity matrix and
$\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$ is the vector of the three
Pauli matrices.

The characterization of an optical SU(2) gate is typically performed by
processing an overcomplete set of 16 projective measurements of the form

$I_{ab}=\absolutevalue{\matrixelement{b}{\hat{U}}{a}}^{2},$ (2)

where $\ket{a}$ and $\ket{b}$ are extracted from the three sets of states
forming the Mutually Unbiased Bases (MUB) of SU(2) [1, 8, 10]. The process
tomography of the gate is then accomplished via an ML approach, i.e., by
minimizing a cost function expressing the distance between the experimental
outcomes $I_{ab}^{\text{exp}}$ and the corresponding theoretical predictions
$I_{ab}^{\text{th}}$ [10, 21]:

$\mathcal{L}=\sum_{ij}(I_{ab}^{\text{exp}}-I_{ab}^{\text{th}})^{2}.$ (3)

This approach may become unfavourably time-consuming and less accurate in the
case of transformations acting on high-dimensional Hilbert spaces. Here, we
consider the case of unitaries, which depend on some additional parameter,
e.g., a spatial variable or a lattice position. More precisely, we assume the
parameters $E,\,n_{1},\,n_{2},\,n_{3}$ to be functions of $\mathbf{r}$, where
$\mathbf{r}$ can be any set of either discrete or continuous variables. In
this current study, we assume position $\mathbf{r}$ as a continuous variable.
However, the discrete case can be obtained by replacing integrals with
summations, i.e., $\int f(x)\,dx\rightarrow\sum_{x}\,f(x)$.

The unknown unitary process acts on quantum states whose Hilbert space is the
tensor product of a qubit space $\mathcal{H}_{i}$, associated with an internal
degree of freedom, and a high-dimensional space $\mathcal{H}_{r}$, associated
with the lattice.

We consider an input state uniformly distributed along $\mathbf{r}$, with an
internal state prepared as one of the MUB states $\ket{a}$. $\ket{a}$ is
assumed to be the positive eigenvector of $\sigma_{1},\sigma_{2}$ or
$\sigma_{3}$. The unknown unitary $U(\mathbf{r})$ acts on this state, and then
a projection onto the same internal state is performed. The probability
distribution along $\mathbf{r}$ is
${I_{aa}(\mathbf{r})=\absolutevalue{\matrixelement{a}{U(\mathbf{r})}{a}}}^{2}$.
Another measurement is performed to retrieve the probability distribution
$\tilde{I}_{aa}(\mathbf{k})$ in the reciprocal space of $\mathbf{r}$. This is
enabled by a QFT of the final state:

$\tilde{I}_{aa}(\mathbf{k})=\absolutevalue{\int\text{d}^{2}r\,\bra{a}U(\mathbf{r})\ket{a}e^{i\mathbf{k}\cdot\mathbf{r}}}^{2}.$
(4)

Figure 1: Fourier Quantum Process Tomography. a) A space-dependent
polarization transformation $U(x,y)$ is implemented via LCMSs. The process
tomography is performed by preparing and projecting onto MUB states. The
resulting intensity distributions are measured in the far field,
$\tilde{I}_{aa}(k_{x},k_{y})$, and image plane, $I_{aa}(x,y)$, of the optical
operator. b) Plots of the optic-axis pattern of two plates used to test the
technique. These are patterned along orthogonal directions to create a complex
2D modulation. c) Reconstruction of the process $U(x,y)=T_{x}(x)T_{y}(y)$
using FQPT and the traditional ML technique. The arrows represent the local
eigenvectors $\bm{n}(x,y)$ and their color is associated with the local
rotation angle $E(x,y)$ (see also Fig. 3d) for details).

One can extract the wave function phase in each plane from the two probability
distributions by applying phase-retrieval techniques. We specifically employed
the Gerchberg–Saxton (GS) algorithm [34]. In particular, we can retrieve the
amplitude ${A_{a}(\mathbf{r})=\sqrt{I_{aa}(\mathbf{r})}}$ and the phase
$\alpha_{a}(\mathbf{r})=\arg(\matrixelement{a}{U(\mathbf{r})}{a})$. However,
the latter is determined up to an unknown constant $\xi_{a}$. The amplitude
and phase are related to the process parameters according to

$\displaystyle\matrixelement{a}{U(\mathbf{r})}{a}$ $\displaystyle=$
$\displaystyle A_{a}(\mathbf{r})e^{i\alpha_{a}(\mathbf{r})+i\xi_{a}}$
$\displaystyle=$
$\displaystyle\cos{E(\mathbf{r})}-i\,n_{a}(\mathbf{r})\sin{E(\mathbf{r})},$

where ${a=1,2,3}$, depending on which Pauli matrix is considered
($\sigma_{a}\ket{a}=\ket{a}$). Thus, we obtain

$\displaystyle E(\mathbf{r})$
$\displaystyle=\arccos{\left(A_{a}(\mathbf{r})\cos{(\alpha_{a}(\mathbf{r})+\xi_{a})}\right)};$
(6a) $\displaystyle n_{a}(\mathbf{r})$
$\displaystyle=-A_{a}(\mathbf{r})\dfrac{\sin{(\alpha_{a}(\mathbf{r})}+\xi_{a})}{\sin{E(\mathbf{r})}}.$
(6b)

Equations (6a)-(6b) show that the extracted parameters depend on the global
phase shifts $\xi_{a}$, which cannot be estimated from the phase retrieval
method. Indeed, any phase that differs from $\alpha_{a}(\mathbf{r})$ by a
constant global shift yields the same measured amplitude in the direct and
reciprocal space. Considering the ambiguity due to the global phase shift, we
list all the possible energy modulations compatible with the measurements. In
practice, we select $N$ values of the global phase shift $\xi_{a,j}=2\pi j/N$,
with ${j=0,1,...,N-1}$. We specifically set ${N=64}$. Only one of the
candidates can best describe the process under investigation. To find it, we
perform an additional measurement in the reciprocal space, obtained by
evolving any input state without projection on the internal degree of freedom,
e.g.,

$\tilde{I}_{0}(\mathbf{k})=\absolutevalue{\int\text{d}^{2}r\,U(\mathbf{r})\ket{b}e^{i\mathbf{k}\cdot\mathbf{r}}}^{2},$
(7)

where $\ket{b}$ can be chosen arbitrarily. Crucially, this last measurement
also provides the normalization factor for all the data. In this way, by
numerically simulating the far field obtained from each of the $N$ possible
SU(2) evolutions, we can isolate the realization associated with the physical
setup by sifting the one that minimizes the distance with the measurement of
Eq. (7).

Figure 2: One-dimensional complex polarization transformations. a) Far-field
intensity distribution $\tilde{I}_{HH}(k_{x})$ for a 1D periodic SU(2)
transformation. A power spectrum $P(m)$ is extracted after discretizing the
reciprocal space. b)-c) Comparison between the experimental total far-field
distribution for a $\ket{L}$ input state and the reconstruction obtained from
FQPT b) and ML c). The insets show the reconstructed near-field amplitude
$\sqrt{I_{HH}(x)}$, compared with the experimental measurement. d)
Reconstructions of the process parameters across one period performed by FQPT
(red curves) and ML (blue curves). The expected periodicity is captured by
both approaches. Figure 3: Two-dimensional complex polarization
transformations. a) Far-field intensity distribution
$\tilde{I}_{HH}(k_{x},k_{y})$ for a 2D periodic SU(2) transformation. A power
spectrum $P(m_{x},m_{y})$ is extracted. b) Comparison between the total far-
field distribution for a $\ket{L}$ input state and the reconstruction obtained
from FQPT and ML. c) Reconstructed near-field amplitude $\sqrt{I_{HH}(x,y)}$,
compared with the experimental measurement. d) Reconstructions of the process
parameters across one period performed by FQPT and ML. The expected
periodicity is captured by both approaches.

We finally remark that, in principle, the method also works with only 2 sets
of measurements. For instance, one could process $I_{11}(\mathbf{r})$ and
$I_{22}(\mathbf{r})$, together with $\tilde{I}_{11}(\mathbf{k})$ and
$\tilde{I}_{22}(\mathbf{k})$, to retrieve $n_{1}(\mathbf{r})$ and
$n_{2}(\mathbf{r})$. In this case, the third component could be directly
computed from the normalization condition:
$n_{3}=\pm\sqrt{1-n_{1}^{2}-n_{2}^{2}}$. Combined with the final measurement
in the reciprocal space, these add up to a minimal set of 5 total
measurements, in perfect agreement with the argument provided in Ref. [18].
However, we observed that experimental noise may yield nonphysical results,
i.e., $n_{3}$ features a non-zero imaginary part at some location. This
suggested integrating the minimal setting with an additional set of
measurements [20, 21, 18].

When accounting for the noise, the proposed technique thus requires only 7
measurements instead of the conventional 16: 3 measurements in the lattice
space, 3 in the reciprocal space, and a last one, still in the reciprocal
space, to fix normalization and remove all the ambiguities on the parameters
of the unitary.

##  III Experimental Results

In photonic setups, a qubit can be encoded into photon polarization, which is
typically manipulated via optical waveplates. In the circular polarization
basis, where $\ket{L}=(1,0)^{T}$ and $\ket{R}=(0,1)^{T}$ are left and right
circular polarization states, respectively, a waveplate $L_{\delta,\theta}$
having the birefringence $\delta$ and the optic axis oriented at $\theta$ with
respect to the horizontal direction can be expressed in the matrix form

$L_{\delta,\theta}=\begin{pmatrix}\cos{(\delta/2)}&i\sin{(\delta/2)}e^{-2i\theta}\\\
i\sin{(\delta/2)}e^{2i\theta}&\cos{(\delta/2)}\end{pmatrix}.$ (8)

A single waveplate thus implements a rotation of an angle $-\delta/2$ around
the equatorial axis ${\bm{n}=(\cos{2\theta},\sin{2\theta},0)}$. Nevertheless,
one can cascade multiple waveplates to implement more general operations [35,
36].

We apply FQPT to periodic polarization transformations induced by complex
LCMSs [22]. These can be modeled as optical waveplates having patterned optic-
axis modulation $\theta=\theta(x,y)$ and fixed, but tunable, birefringence
[37]. In particular, we test our method with LCMSs featuring high spatial-
frequency modulations along the $x$ and $y$ directions. This scenario is
experimentally more challenging than the one addressed in Ref. [18], where
simple combinations of polarization gratings were considered. We benchmark our
technique against the standard ML approach processing a whole set of 16
polarimetric measurements (all taken in the near field), taking into account
both the timing and the accuracy of the final reconstruction. For the
minimization, we employed the _NMinimize_ routine from Wolfram Mathematica
[38].

The experiments are realized with the setup sketched in Fig. 1a). A Ti:Sa
laser (wavelength ${\lambda=810}$ nm) is coupled to a single-mode fiber. The
output Gaussian mode is magnified with a telescope lens system,
$f_{1}=125\text{ mm}$ and $f_{2}=200$ mm (not shown in the figure). The beam
waist is measured to be $2.6\pm 0.1\text{ mm}$. In this way, the overall beam
size is larger than the largest periodicity on the plates, that is
$\Lambda=2.5\text{ mm}$. A combination of a linear polarizer ($\text{P}_{1}$),
a half-wave plate ($\text{HWP}_{1}$) and a quarter-wave plate
($\text{QWP}_{1}$) is needed to prepare any input polarization state. The beam
then propagates through one or multiple LCMSs, engineering a space-dependent
SU(2) optical operator. Another set
$\text{QWP}_{2}$-$\text{HWP}_{2}$-$\text{P}_{2}$ is adjusted to project onto
any state. The last element ($\text{P}_{2}$) is removed when performing the
measurement of Eq. (7). Each polarimetric measurement is collected on a CCD
camera, placed either after a 4$f$ system (${f=150\text{ mm}}$) or in the
focal plane of a lens (${f^{\prime}=250\text{ mm}}$), depending on if the
measurement is realized in the near field or in the far field, respectively.
Recall that the far-field light distribution is proportional to the transverse
momentum distribution, i.e., to the Fourier transform of the input field.

The first experiment is realized with a single LCMS displaying a complex
periodic modulation along the $x$ axis (see Fig. 1b)). The optical retardation
is set at $\delta=\pi$. The periodicity of the process simplifies the
analysis, as only a discrete spectrum of Fourier components is expected to
appear in the far field for all the measurements. Furthermore, in this one-
dimensional (1D) realization, the intensity modulations can also be integrated
along the $y$ direction to mitigate experimental imperfections.

Figure 2a) shows the far-field intensity distributions recorded for the
$\matrixelement{H}{U}{H}$ configuration, with
${\ket{H}=(\ket{L}+\ket{R})/\sqrt{2}}$. As discussed above, the far-field
distribution can be discretized and a normalized power spectrum $P(m)$ is
extracted. Figures 2b)-c) illustrate the comparison between the experimentally
measured total far-field distribution for a $\ket{L}$ input state (see Eq.
(7)) and the reconstruction performed via FQPT and ML, respectively. The
agreement with the experimental observation is quantified in terms of the
similarity estimator
${s=(\sum_{m}\sqrt{P_{\text{exp}}(m)P_{\text{rec}}(m)})^{2}}$, where
$P_{\text{exp}}(m)$ and $P_{\text{rec}}(m)$ are the (normalized) experimental
and reconstructed far-field distributions. We obtain
${s_{\text{FQPT}}=97.2\%}$ and ${s_{\text{ML}}=97.0\%}$. Another metric that
can be considered is the absolute distance between the two distributions,
computed as
${\Delta=(\sum_{m}\absolutevalue{P_{\text{exp}}(m)-P_{\text{rec}}(m)})^{2}}$.
We report ${\Delta_{\text{FQPT}}=0.088}$ and ${\Delta_{\text{ML}}=0.083}$. In
both figures, the insets show the same comparison for the 1D near-field
amplitude ${\sqrt{I_{HH}(x)}}$. Reconstructions of individual parameters are
plotted in Fig. 2d). The agreement between the predictions of the two methods
is quantified in terms of the fidelity
${F=\absolutevalue{\Tr(U_{\text{ML}}^{{\dagger}}U_{\text{FQPT}})}/n}$, where
${n=2}$ is the dimension of the internal degree of freedom [39]. An excellent
average fidelity is obtained, $\bar{F}=95.7\%$, where $\bar{F}$ denotes the
average fidelity computed over all the pixels. These results prove that both
methods provide reliable reconstructions. It must be noted that FQPT achieves
satisfactory reconstructions by only processing a near-optimal set of
measurements. Moreover, a brute-force minimization approach tends to jump
between the parameters associated with processes $U$ and ${e^{i\pi}U=-U}$, as
these both generate the same experimental outcomes [40, 41, 18, 42]. For this
reason, a continuity constraint must be enforced between consecutive pixels.
Conversely, our method is assumption-free as it does not rely on any a-priori
hypothesis on the process parameters. The technique is also extremely
efficient on the computational level. The total times required for a complete
reconstruction are ${t_{\text{FQPT}}\approx 1\text{ min}}$ and
${t_{\text{ML}}\approx 30\text{ min}}$.

Recall that the GS algorithm is executed to retrieve the unknown phases. This
algorithm is based on an iterative strategy that presents intrinsic
limitations [43]. The convergence speed is sensitive to the initial guess of
the phase, and the convergence to a global minimum is not guaranteed.
Moreover, the presence of noise in the input data can severely affect the
accuracy of phase retrieval. To overcome these limitations, for each set of
polarimetric measurements, the algorithm is run ${N_{T}^{\text{1D}}=100}$
times for ${N_{I}^{\text{1D}}=1000}$ iterations, randomly initializing the
phase guess at each trial. In doing so, the best phase reconstruction can be
selected (up to global shifts) as the one minimizing the total distance
between the reconstructed and measured near-field amplitudes.

The second experiment is realized by cascading the previous LCMS with a second
one, patterned along the $y$ direction (see Fig. 1b)). The sequence thus
implements a two-dimensional (2D) periodic SU(2) transformation. The
birefringence setting is ${\delta_{1}=\pi/2}$ and ${\delta_{2}=\pi}$. Figure
3a) shows the experimental far-field distribution
${\tilde{I}_{HH}(k_{x},k_{y})}$, from which a 2D power spectrum
$P(m_{x},m_{y})$ is extracted. The total far-field reconstructions for a
$\ket{L}$ input state are plotted in Fig. 3b), where the comparison with the
experimental result is also provided. We obtain ${s_{\text{FQPT}}=88.5\%}$ and
${s_{\text{ML}}=96.7\%}$, with a total error of ${\Delta_{\text{FQPT}}=0.080}$
and ${\Delta_{\text{ML}}=0.013}$. The slight deterioration in the final
prediction can be ascribed to the reduced performance of the GS algorithm in
two spatial dimensions [44]. Figure 3c) shows the near-field amplitude
${\sqrt{I_{HH}(x,y)}}$, as reconstructed from FQPT and ML, compared with the
experimental measurement. In this case, the GS algorithm runs
${N_{T}^{\text{2D}}=50}$ times for ${N_{I}^{\text{2D}}=500}$ iterations. In
this realization, we compress the experimental images by integrating light
intensity over ${11\text{ pixel}\times 11\text{ pixel}}$ regions equally
distributed on the camera. This allows for both minimizing the errors due to
local intensity fluctuations in the image area and keeping the computation
time within the same range as the 1D experiment. The reconstructions of
individual parameters are shown in Fig. 3d). An alternative visualization of
the reconstructed process in terms of the local eigenvector and rotation angle
is provided in Fig. 1c). A good agreement is observed between the two
predictions, with average fidelity $\bar{F}=91.4\%$.

##  IV Discussion and conclusions

In this work, a new technique for fast and accurate Quantum Process Tomography
is demonstrated. This is accomplished via a non-interferometric scheme
requiring no a-priori information on the unknown operator. In the case of
complex polarization transformations, our method achieves performances very
close to the standard tomography based on an overcomplete set of measurements.
It offers experimental advantage, only requiring a near-minimal set of
measurements, and computational speed-up, outperforming the standard approach
by at least one order of magnitude. It appears naturally suitable for all
experimental setups that allow for easy access to conjugate domains. For this
reason, this method can also be implemented on other physical platforms, such
as quantum circuits for atoms [45, 46] and electron beams [47, 48].

Nevertheless, we believe that the performance of the method can be further
improved by adopting optimized strategies for phase retrieval. For instance,
Convolutional Neural Networks represent a promising solution [49, 50, 51]. We
also expect that FQPT can detect non-unitary evolutions [52] if equipped with
some minor modifications, such as including multiple intermediate planes in
the final analysis. At the same time, it would be interesting to adapt the
present method to retrieve multi-photon gates and complex operations in high-
dimensional Hilbert spaces.

In principle, this procedure can also be applied to processes acting on a
$m\times d$ dimensional space, having the irreducible form
$\bigoplus_{i=1}^{d}\mathcal{U}^{(i)}$, with
${\mathcal{U}^{(i)}\in\text{SU}(m)}$. In this case, one must consider the
decomposition of $\mathcal{U}^{(i)}$ in the generators of $\text{SU}(m)$, the
generalized Gell-Mann matrices. Although the analytical relations between the
process parameters and the measurement outcomes become significantly more
complicated, the number of measurements required by FQPT will still be
optimal. This will be investigated in successive works.

## Acknowledgements

This work was supported by the Canada Research Chair (CRC) Program, NRC-
uOttawa Joint Centre for Extreme Quantum Photonics (JCEP) via the Quantum
Sensors Challenge Program at the National Research Council of Canada, and
Quantum Enhanced Sensing and Imaging (QuEnSI) Alliance Consortia Quantum
grant.

## References

  * Chuang and Nielsen [1997] I. L. Chuang and M. A. Nielsen, J. Mod. Opt. 44, 2455 (1997).
  * Childs  _et al._ [2001] A. M. Childs, I. L. Chuang, and D. W. Leung, Phys. Rev. A 64, 012314 (2001).
  * Myrskog  _et al._ [2005] S. H. Myrskog, J. K. Fox, M. W. Mitchell, and A. M. Steinberg, Phys. Rev. A 72, 013615 (2005).
  * Roos  _et al._ [2004] C. F. Roos, G. P. T. Lancaster, M. Riebe, H. Häffner, W. Hänsel, S. Gulde, C. Becher, J. Eschner, F. Schmidt-Kaler, and R. Blatt, Phys. Rev. Lett. 92, 220402 (2004).
  * Riebe  _et al._ [2006] M. Riebe, K. Kim, P. Schindler, T. Monz, P. O. Schmidt, T. K. Körber, W. Hänsel, H. Häffner, C. F. Roos, and R. Blatt, Phys. Rev. Lett. 97, 220407 (2006).
  * Yamamoto  _et al._ [2010] T. Yamamoto, M. Neeley, E. Lucero, R. C. Bialczak, J. Kelly, M. Lenander, M. Mariantoni, A. D. O’Connell, D. Sank, H. Wang, M. Weides, J. Wenner, Y. Yin, A. N. Cleland, and J. M. Martinis, Phys. Rev. B 82, 184515 (2010).
  * Bialczak  _et al._ [2010] R. C. Bialczak, M. Ansmann, M. Hofheinz, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, H. Wang, J. Wenner, M. Steffen, A. N. Cleland, and J. M. Martinis, Nat. Phys. 6, 409 (2010).
  * Mitchell  _et al._ [2003] M. W. Mitchell, C. W. Ellenor, S. Schneider, and A. M. Steinberg, Phys. Rev. Lett. 91, 120402 (2003).
  * Altepeter  _et al._ [2003] J. B. Altepeter, D. Branning, E. Jeffrey, T. C. Wei, P. G. Kwiat, R. T. Thew, J. L. O’Brien, M. A. Nielsen, and A. G. White, Phys. Rev. Lett. 90, 193601 (2003).
  * O’Brien  _et al._ [2004] J. L. O’Brien, G. J. Pryde, A. Gilchrist, D. F. V. James, N. K. Langford, T. C. Ralph, and A. G. White, Phys. Rev. Lett. 93, 080502 (2004).
  * Lobino  _et al._ [2008] M. Lobino, D. Korystov, C. Kupchak, E. Figueroa, B. C. Sanders, and A. I. Lvovsky, Science 322, 563 (2008).
  * Bongioanni  _et al._ [2010] I. Bongioanni, L. Sansoni, F. Sciarrino, G. Vallone, and P. Mataloni, Phys. Rev. A 82, 042307 (2010).
  * Rahimi-Keshari  _et al._ [2013] S. Rahimi-Keshari, M. A. Broome, R. Fickler, A. Fedrizzi, T. C. Ralph, and A. G. White, Opt. Express 21, 13450 (2013).
  * Ndagano  _et al._ [2017] B. Ndagano, B. Perez-Garcia, F. S. Roux, M. McLaren, C. Rosales-Guzman, Y. Zhang, O. Mouane, R. I. Hernandez-Aranda, T. Konrad, and A. Forbes, Nat. Phys. 13, 397 (2017).
  * Antón  _et al._ [2017] C. Antón, P. Hilaire, C. A. Kessler, J. Demory, C. Gómez, A. Lemaître, I. Sagnes, N. D. Lanzillotti-Kimura, O. Krebs, N. Somaschi, P. Senellart, and L. Lanco, Optica 4, 1326 (2017).
  * Jacob  _et al._ [2018] K. V. Jacob, A. E. Mirasola, S. Adhikari, and J. P. Dowling, Phys. Rev. A 98, 052327 (2018).
  * Bouchard  _et al._ [2019] F. Bouchard, F. Hufnagel, D. Koutný, A. Abbas, A. Sit, K. Heshami, R. Fickler, and E. Karimi, Quantum 3, 138 (2019).
  * Di Colandrea  _et al._ [2023a] F. Di Colandrea, L. Amato, R. Schiattarella, A. Dauphin, and F. Cardano, Opt. Express 31, 31698 (2023a).
  * Le Roy-Brehonnet and Le Jeune [1997] F. Le Roy-Brehonnet and B. Le Jeune, Prog. Quantum Electron. 21, 109 (1997).
  * James  _et al._ [2001] D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, Phys. Rev. A 64, 052312 (2001).
  * Aiello  _et al._ [2006] A. Aiello, G. Puentes, D. Voigt, and J. P. Woerdman, Opt. Lett. 31, 817 (2006).
  * Di Colandrea  _et al._ [2023b] F. Di Colandrea, A. Babazadeh, A. Dauphin, P. Massignan, L. Marrucci, and F. Cardano, Optica 10, 324 (2023b).
  * Goodman [2005] J. W. Goodman,  _Introduction to Fourier optics_ (Roberts and Company publishers, 2005).
  * Wang  _et al._ [2020] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, Nat. Photonics 14, 273 (2020).
  * Politi  _et al._ [2009] A. Politi, J. C. Matthews, and J. L. O’Brien, Science 325, 1221 (2009).
  * Crespi  _et al._ [2016] A. Crespi, R. Osellame, R. Ramponi, M. Bentivegna, F. Flamini, N. Spagnolo, N. Viggianiello, L. Innocenti, P. Mataloni, and F. Sciarrino, Nat. Commun. 7, 10469 (2016).
  * Flamini  _et al._ [2018] F. Flamini, N. Viggianiello, T. Giordani, M. Bentivegna, N. Spagnolo, A. Crespi, G. Corrielli, R. Osellame, M. A. Martin-Delgado, and F. Sciarrino, J. Opt. 20, 074001 (2018).
  * Sansoni  _et al._ [2012] L. Sansoni, F. Sciarrino, G. Vallone, P. Mataloni, A. Crespi, R. Ramponi, and R. Osellame, Phys. Rev. Lett. 108, 010502 (2012).
  * Heilmann  _et al._ [2014] R. Heilmann, M. Gräfe, S. Nolte, and A. Szameit, Sci. Rep. 4, 4118 (2014).
  * Pitsios  _et al._ [2017] I. Pitsios, F. Samara, G. Corrielli, A. Crespi, and R. Osellame, Sci. Rep. 7, 11342 (2017).
  * Weinstein  _et al._ [2001] Y. S. Weinstein, M. Pravia, E. Fortunato, S. Lloyd, and D. G. Cory, Phys. Rev. Lett. 86, 1889 (2001).
  * Fowler  _et al._ [2004] A. G. Fowler, S. J. Devitt, and L. C. Hollenberg, Quantum Info. Comput. 4, 237–251 (2004).
  * Mariantoni  _et al._ [2011] M. Mariantoni, H. Wang, T. Yamamoto, M. Neeley, R. C. Bialczak, Y. Chen, M. Lenander, E. Lucero, A. D. O’Connell, D. Sank, _et al._ , Science 334, 61 (2011).
  * Gerhberg and Saxton [1972] R. Gerhberg and W. Saxton, Optik 35, 237 (1972).
  * Simon and Mukunda [1990] R. Simon and N. Mukunda, Phys. Lett. A 143, 165 (1990).
  * Sit  _et al._ [2017] A. Sit, L. Giner, E. Karimi, and J. S. Lundeen, J. Opt. 19, 094003 (2017).
  * Rubano  _et al._ [2019] A. Rubano, F. Cardano, B. Piccirillo, and L. Marrucci, J. Opt. Soc. Am. B 36, D70 (2019).
  * Wolfram Research [2021] Wolfram Research,  _Numerical Nonlinear Global Optimization_ (2021) https://reference.wolfram.com/language/tutorial/ConstrainedOptimizationGlobalNumerical.html.
  * Wang  _et al._ [2009] X. Wang, Z. Sun, and Z. D. Wang, Phys. Rev. A 79, 012105 (2009).
  * Fläschner  _et al._ [2016] N. Fläschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. Lühmann, K. Sengstock, and C. Weitenberg, Science 352, 1091 (2016).
  * Tarnowski  _et al._ [2019] M. Tarnowski, F. N. Ünal, N. Fläschner, B. S. Rem, A. Eckardt, K. Sengstock, and C. Weitenberg, Nat. Commun. 10, 1728 (2019).
  * Yi  _et al._ [2023] C.-R. Yi, J. Yu, H. Yuan, R.-H. Jiao, Y.-M. Yang, X. Jiang, J.-Y. Zhang, S. Chen, and J.-W. Pan, Phys. Rev. Res. 5, L032016 (2023).
  * Fienup [1982] J. R. Fienup, Appl. Opt. 21, 2758 (1982).
  * Zhao and Chi [2020] T. Zhao and Y. Chi, Entropy 22 (2020).
  * Guidoni and Verkerk [1999] L. Guidoni and P. Verkerk, J. opt., B Quantum semiclass. opt. 1, R23 (1999).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Science 357, 995 (2017).
  * Harris  _et al._ [2015] J. Harris, V. Grillo, E. Mafakheri, G. C. Gazzadi, S. Frabboni, R. W. Boyd, and E. Karimi, Nat. Phys. 11, 629 (2015).
  * Grillo  _et al._ [2015] V. Grillo, G. C. Gazzadi, E. Mafakheri, S. Frabboni, E. Karimi, and R. W. Boyd, Phys. Rev. Lett. 114, 034801 (2015).
  * White  _et al._ [2021] J. White, S. Wang, W. Eschen, and J. Rothhardt, Opt. Express 29, 9283 (2021).
  * Li  _et al._ [2021] B. Li, O. K. Ersoy, C. Ma, Z. Pan, W. Wen, Z. Song, and W. Gao, J. Mod. Opt. 68, 1108 (2021).
  * Ding  _et al._ [2022] G. Ding, W. Xiong, P. Wang, Z. Huang, Y. He, J. Liu, Y. Li, D. Fan, and S. Chen, J. Opt. 24, 025701 (2022).
  * Wang  _et al._ [2023] C. Wang, Z. Fu, W. Mao, J. Qie, A. D. Stone, and L. Yang, Adv. Opt. Photon. 15, 442 (2023).