# Quantum Simulation of the Radical Pair Dynamics of the Avian Compass

Yiteng Zhang  Independent Researcher  Zixuan Hu  Department of Chemistry,
Department of Physics, and Purdue Quantum Science and Engineering Institute
Yuchen Wang  Department of Chemistry, Department of Physics, and Purdue
Quantum Science and Engineering Institute  Sabre Kais<EMAIL_ADDRESS>Department of Chemistry, Department of Physics, and Purdue Quantum Science and
Engineering Institute

###### Abstract

The simulation of open quantum dynamics on quantum circuits has attracted wide
interests recently with a variety of quantum algorithms developed and
demonstrated. Among these, one particular design of a unitary-dilation-based
quantum algorithm is capable of simulating general and complex physical
systems. In this paper, we apply this quantum algorithm to simulating the
dynamics of the radical pair mechanism in the avian compass. This application
is demonstrated on the IBM QASM quantum simulator. This work is the first
application of any quantum algorithm to simulating the radical pair mechanism
in the avian compass, which not only demonstrates the generality of the
quantum algorithm, but also opens new opportunities for studying the avian
compass with quantum computing devices.

##  I Introduction

An open quantum system is a quantum system that interacts with an external
environment or bath. The interaction between the system and the bath is often
too complex to be simulated exactly, and thus requires approximations to
average out the effects of the bath – this results in the non-unitary dynamics
of open quantum systems. Simulating the dynamics of quantum systems has been a
main focus of quantum computing research [1, 2, 3, 4, 5, 6], yet relatively
few quantum algorithms have been developed for simulating the dynamics of open
quantum systems [7, 8, 9, 10, 11, 12, 13, 14]. To this end, we have developed
and demonstrated a general quantum algorithm for open quantum dynamics [15,
16, 17] that is capable of simulating general and complex physical systems.
The quantum algorithm leverages the Sz.-Nagy unitary dilation approach to
convert non-unitary time evolution operators into corresponding unitary
operators, which can then be implemented on a quantum circuit. This quantum
algorithm has been applied to a variety of physical systems, including the
amplitude damping channel described by the Kraus representation [15], the
Jaynes-Cummings model described by the Kraus representation [18], the Fenna-
Matthews-Olson (FMO) complex described by the Lindblad master equation [16],
and the spin-boson model described by the generalized quantum master equation
(GQME) [17].

In this work, we apply the general quantum algorithm to simulating the radical
pair mechanism in the avian compass and further demonstrate its generality.
The radical pair mechanism (RPM) is a theory proposed to explain the magneto-
reception and navigation abilities of certain bird species [19, 20, 21]. Many
animals possess extraordinary abilities to sense the direction by perceiving
the geomagnetic field. This is probably the result of natural selection over a
very long time of evolution, since the ability to sense the direction is
crucial for certain animals to find their habitats, such as migratory birds
that change habitats from season to season.

In brief, the RPM involves two spatially separated electrons, which are
correlated with each other and affected by the external weak magnetic field
and internal nuclear spin couplings. The basic scheme of the RPM includes
three stages: 1) the photons with certain energies activate a certain type of
molecules located in the bird’s eyes, enabling an electron transfer reaction
and generating a radical pair in the singlet state; 2) the state of the
radical pair converts between the singlet state and the triplet state under
the influence of the external magnetic field (the geomagnetic field) and the
internal magnetic field (the hyperfine coupling effects); 3) the radical pairs
in different states will generate different chemical products which can induce
a detectable signal for birds to recognize the direction [22]. The RPM is a
promising hypothesis that can explain the three unusual properties of the
avian compass: 1. the inclination compass: the functional mode of the avian
magnetic compass is based on the inclination of the field lines instead of
their polarity [23, 24, 25, 26]; 2. the light dependence: light with an energy
above a certain threshold is needed for the RPM to work [27, 28, 29, 30, 31,
32, 33]; 3. the narrow range of responsive magnetic field intensities: both
higher and lower magnetic fields will disable birds’ ability of navigation
[34].

To understand the RP mechanism, E. M. Gauger $et.\ al.$ proposed a way to
model the dynamics of the RPM system with a Lindblad master equation by adding
two “shelving states” for the singlet yield and triple yield [35]. In the
following, we treat the same Lindblad formulation of the RPM dynamics with our
general quantum algorithm for open quantum dynamics and simulate the RPM
dynamics on the IBM QASM quantum simulator. To our best knowledge, this is the
first ever demonstration of any quantum algorithm applied to simulating the
RPM dynamics. This work not only shows the generality of the quantum
algorithm, but also opens new potential ways to study the avian compass with
quantum computing devices.

##  II Methodology

###  II.1 The General Quantum Algorithm for Open Quantum Dynamics

The general quantum algorithm for open quantum dynamics has different versions
that can evolve the Kraus representation [15], the Lindblad master equation
[16], and the generalized quantum master equation [17]. In this work we use
the version for the Lindblad master equation [16].

We first review how the generalized quantum algorithm converts the non-unitary
matrix that encodes the open quantum dynamics into a unitary evaluation based
on the Sz.-Nagy unitary dilation procedure [15, 16]. We assume the initial
density matrix that describes the physical system is composed of a set of
unique pure quantum states $|\phi_{i}\rangle$ that are weighted by their
corresponding probabilities $p_{i}$:

$\rho=\sum_{i}p_{i}|\phi_{i}\rangle\langle\phi_{i}|$

we want to simulate the time evolution of $\rho(t)$ given the initial $\rho$
and the Kraus operators $M_{k}$’s. This task can be achieved by preparing each
input state $|\phi_{i}\rangle$ in a vector form $v_{i}$ in a given basis and
then building a quantum circuit that generates the quantum state:

$|\phi_{ik}(t)\rangle=M_{k}v_{i}\xrightarrow{\text{unitary
dilation}}U_{M_{k}}(v_{i}^{T},0,\cdots,0)^{T}.$ (1)

The $U_{M_{k}}$ is generated via the 1-dilation of $M_{k}$:

$U_{M_{k}}=\begin{pmatrix}M_{k}&D_{M_{k}^{\dagger}}\\\
D_{M_{k}}&-M_{k}^{\dagger}\end{pmatrix},$ (2)

where
$D_{M_{k}}=\sqrt{I-M_{k}^{\dagger}M_{k}},\;D_{M_{k}^{\dagger}}=\sqrt{I-M_{k}M_{k}^{\dagger}}$
[36]. After obtaining each $|\phi_{ik}(t)\rangle$, we can calculate the
population of each basis state in the current basis from the diagonal vector:

$\mathrm{diag}(\rho(t))=\sum_{ik}p_{i}\cdot\mathrm{diag}(|\phi_{ik}(t)\rangle\langle\phi_{ik}(t)|),$
(3)

where $\mathrm{diag}(|\phi_{ik}(t)\rangle\langle\phi_{ik}(t)|)$ can be
efficiently obtained by applying projection measurements on the first half
subspace of $U_{M_{k}}(v_{i}^{T},0,\cdots,0)^{T}$.

For the dynamics of an open quantum system, the time evolution of the density
matrix can be represented as:

$\rho(s+\delta s)=\sum_{k}M_{sk}\rho(s)M_{sk}^{\dagger}$ (4)

where $\rho(s)$ is the density matrix at time step $s$, and $\delta s$ is
considered as the discrete time step, during which the Kraus operators
$M_{sk}$ are assumed to be constant. The formula eq.(4) can be used
iteratively until reaching the time of interest. Explicitly, the dynamics of
the density matrix is described as:

$\displaystyle\rho(1)$ $\displaystyle=\rho(1\delta
t)=\sum_{k}M_{0k}\rho(0)M_{0k}^{\dagger}$ (5) $\displaystyle\rho(2)$
$\displaystyle=\rho(2\delta t)=\sum_{k}M_{1k}\rho(1)M_{1k}^{\dagger}$ (6)
$\displaystyle=\sum_{k}\sum_{j}M_{1k}M_{0j}\rho(0)M_{0j}^{\dagger}M_{1k}^{\dagger}$
$\displaystyle\rho(3)$ $\displaystyle=\rho(3\delta
t)=\sum_{k}M_{2k}\rho(2)M_{2k}^{\dagger}$ (7)
$\displaystyle=\sum_{k}\sum_{j}\sum_{i}M_{2k}M_{1j}M_{0i}\rho(0)M_{0i}^{\dagger}M_{1j}^{\dagger}M_{2k}^{\dagger}$
$\displaystyle\dots$

Here without losing any generality, the Kraus operators $M_{sk}$ are indexed
by the time step $s$, which allows each time step to have a different set of
Kraus operators. However, as will be discussed in Section II.3, the RPM
dynamical model used in this work is a Markovian process described by the
Lindblad master equation, therefore all the time steps have the same set of
Kraus operators $M_{k}$ where the time step index $s$ has been removed.

###  II.2 The Radical Pair Mechanism Theory and Dynamics

Figure 1: The basic scheme of the radical pair mechanism. After absorbing
light, donor (D) and acceptor (A) molecules form radical pairs in its singlet
state. Then, under the influence of the magnetic fields, the states of the
radical pair interconverse between the singlet states and triplet states.
Finally, the singlet and triplet radical pairs end up with different products.

The basic scheme of the RPM is shown in Fig. 1. To simplify the fields, we
assume that only the electron near the donor interacts with the nucleus, and
the electron away from the donor is not affected by the anisotropic hyperfine
coupling [35]. Therefore, the Hamiltonian of the system is

$H=\gamma[\hat{I}\cdot\mathbf{A}\cdot\hat{S}_{1}+\mathbf{B}\cdot(\hat{S}_{1}+\hat{S}_{2})]$
(8)

where $\mathbf{A}$ is the anisotropic hyperfine tensor coupling the nucleus
and one of two spatially separated electrons, and
$\mathbf{A}=diag\\{A_{x},A_{y},A_{z}\\}$ with $A_{x}=A_{y}=A_{z}/2$;
$\mathbf{B}=B_{0}(\cos\varphi\sin\theta,\sin\varphi\sin\theta,\cos\theta)$,
and $B_{0}$ is the magnitude of the geomagnetic field; $\varphi$ is the angle
between the x-axis of the radical pair and the external magnetic field;
$\theta$ is the angle between the z-axis of the radical pair and the external
magnetic field; $\gamma=\frac{1}{2}\mu_{0}g$, and $\mu_{0}$ is the Bohr
magneton and $g=2$ is the electron-spin $g$-factor.

To model the dynamics of the system with a quantum master equation
formulation, two “shelving states” were added to the 8-dimensional Hilbert
space of the three spins (two electron spins and one nuclear spin) [35]. We
employ operators as shown in Eq.(9) to represent the spin-selective relaxation
into the singlet shelf $|S\rangle$ from the electron singlet state, or the
triplet shelf $|T\rangle$ from the electron triplet state. The final
populations of $|S\rangle$ and $|T\rangle$ give the singlet and triplet
yields.

With the usual definition of singlet $|s\rangle$ and triplet states
$|t_{i}\rangle$ in the electronic subspace, while $\mid\uparrow\rangle$ and
$\mid\downarrow\rangle$ describing the states of the nuclear spin, we define
the following decay operators:

$\displaystyle P_{1}=P_{S,\uparrow}=|S\rangle\langle
s,\uparrow|,P_{2}=P_{T_{0},\uparrow}=|T\rangle\langle t_{0},\uparrow|$ (9)
$\displaystyle P_{3}=P_{T_{+},\uparrow}=|T\rangle\langle
t_{+},\uparrow|,P_{4}=P_{T_{-},\uparrow}=|T\rangle\langle t_{-},\uparrow|$
$\displaystyle P_{5}=P_{S,\downarrow}=|S\rangle\langle
s,\downarrow|,P_{6}=P_{T_{0},\downarrow}=|T\rangle\langle t_{0},\downarrow|$
$\displaystyle P_{7}=P_{T_{+},\downarrow}=|T\rangle\langle
t_{+},\downarrow|,P_{8}=P_{T_{-},\downarrow}=|T\rangle\langle
t_{-},\downarrow|$

Parameter Details
---
Symbol |  Description |  Values
$A_{x}$ |  Anisotropic hyperfine tensor |  $1\times 10^{-4}T$
$B_{0}$ |  Magnitude of the geomagnetic field |  $5\times 10^{-5}T$
$\gamma$ |  Half of the product of the Bohr magneton and electron-spin $g$-factor |  $9.27\times 10^{-24}J/T$
$\hbar$ |  Reduced Planck constant |  $1.05457\times 10^{-32}J\cdot s$
$\varphi$ |  Angle between $x$-axis of the radical pair and the magnetic field |  $0$
$k_{d}$ |  Decay Rate of the singlet and triplet states |  $1\times 10^{4}s^{-1}$
Table 1:  Parameter values used in the calculation.

This gives a standard Lindblad master equation:

$\dot{\rho}=-\frac{i}{\hbar}[H,\rho]+k_{d}\sum_{i=1}^{8}[P_{i}\rho
P_{i}^{\dagger}-\frac{1}{2}(P_{i}^{\dagger}P_{i}\rho+\rho
P_{i}^{\dagger}P_{i})],$ (10)

where $k_{d}$ is the decay rate of the singlet and triplet states. Note the
decay rate $k_{d}$ is independent of the radical pair states, so we have
assigned the same decay rate $k_{d}$ to all eight projectors.

###  II.3 Calculation and Simulation

Now to apply the general quantum algorithm to the RPM dynamics, we first
consider the non-unitary part on the right side of Eq. (10) which can be
rewritten as:

$\frac{\delta\rho(t)}{\delta
t}=\mathcal{L}(\rho)=k_{d}\sum_{i=1}^{8}[P_{i}\rho
P_{i}^{\dagger}-\frac{1}{2}\\{P_{i}^{\dagger}P_{i},\rho\\}]$ (11)

Given a very small $\delta t$, Eq. (11) becomes:

$\displaystyle\rho(t+\delta t)-\rho(t)=$ $\displaystyle k_{d}\delta
t\sum_{i=1}^{8}[P_{i}\rho(t)P_{i}^{\dagger}-\frac{1}{2}\\{P_{i}^{\dagger}P_{i},\rho(t)\\}]$
(12) $\displaystyle+\mathcal{O}(\delta t^{2}).$

Figure 2: The comparison of the singlet(in blue) and triplet(in red) yields
between the results obtained from exact calculation and the quantum simulation
of the RPM dynamics. The exact curves are generated from the cubic
interpolation of the exact calculation of the yields at each data points. The
dots represent the results simulated by the general quantum algorithm as
implemented on the IBM QASM simulator. The parameters used are shown in Table
1. The yields are calculated around $7.5\times 10^{-4}s$ after the system has
already reached the steady-state. The y-axis is the final singlet/triplet
yields – i.e. the populations of singlet/triplet shelf state; and the x-axis
is the angle between z-axis of the radical pair and the magnetic field. Figure
3: The comparison of the dynamics of the singlet(in blue) and triplet(in red)
yields between the exact and quantum simulated results. The exact curves are
generated from the cubic interpolation of the exact calculation of the yields
at each data points. The dots represent the results simulated by the general
quantum algorithm as implemented on the IBM QASM simulator. After about
$2\times 10^{-4}s$, both yields reach steady-state. The y-axis is the final
singlet/triplet yields – i.e. the populations of singlet/triplet shelf state;
and the x-axis is the time.

Now assuming $M_{0}=\sqrt{\mathbf{I}-\frac{1}{2}k_{d}\delta
t\sum_{k=1}^{8}P_{k}^{\dagger}P_{k}}$ and $M_{k}=\sqrt{k_{d}\delta t}P_{k}$
for $k>0$, and ignoring the second order of $\delta t$ as $\delta t\rightarrow
0$, Eq. ((12)) can be rewritten as:

$\displaystyle\rho(t+\delta
t)=M_{0}\rho(t)M_{0}^{\dagger}+\sum_{k=1}^{8}M_{k}\rho(t)M_{k}^{\dagger}$ (13)

Eq.(13) is in the same form of Eq.(4), thus knowing the initial state
$\rho(0)$ we can evolve the density matrix to a certain time with the
iterative procedure described in Eq.(5) to Eq.(7).

In the meantime,

$\displaystyle\sum_{k=0}^{8}M_{k}^{\dagger}M_{k}$
$\displaystyle=M_{0}^{\dagger}M_{0}+\sum_{k=1}^{8}M_{k}^{\dagger}M_{k}$ (14)
$\displaystyle=\mathbf{I}-k_{d}\delta
t\sum_{k=1}^{8}P_{k}^{\dagger}P_{k}++\mathcal{O}(\delta t^{2})+k_{d}\delta
t\sum_{k=1}^{8}P_{k}^{\dagger}P_{k}$
$\displaystyle=\mathbf{I}+\mathcal{O}(\delta t^{2})$

When $\delta t\rightarrow 0$, according to the above equation, we have:

$M_{0}=\sqrt{\mathbf{I}-k_{d}\delta t\sum_{k=1}^{8}P_{k}^{\dagger}P_{k}}$ (15)

With Eq. (15), we can formulate the Kraus operators to satisfy the condition
$\sum_{k}M_{k}^{\dagger}M_{k}=\mathbf{I}$. Thus, we have defined all the 9
Kraus operators $M_{k}$ required to describe the RPM dynamics.

There is one additional term containing the Hamiltonian,
$-\frac{i}{\hbar}[H,\rho]$, in Eq. (10) as compared with Eq. (11). This
“oscillating part” of the dynamics is unitary and thus can be easily realized
by multiplying each Kraus operator by a unitary matrix obtained through the
diagonalization of the Hamiltonian [16].

With the parameters in Table 1, we simulated the RPM dynamics by iteratively
applying Eq. (13) on the IBM QASM quantum simulator and then used the output
results to calculate the singlet and triplet yields. Also, we assume the
initial state of the two electron spins is
$\frac{1}{\sqrt{2}}(\mid\uparrow\downarrow\rangle-\mid\downarrow\uparrow\rangle)$,
and the initial state of the nuclear spin is
$\frac{1}{\sqrt{2}}(\mid\uparrow\rangle-\mid\downarrow\rangle)$. In our
simulation, the time interval $\delta t$ is set as $0.5/k_{d}=5\times
10^{-5}s$. We then apply the procedure shown in Eq.(5), Eq.(6) and Eq.(7). As
mentioned above, the populations of $\mid S\rangle$ and $\mid T\rangle$ is the
singlet and triplet yields respectively, after the system reaches the steady
state. The populations are calculated by the procedure explained in Eq. (3),
where the diagonal elements of the density matrix are obtained by projection
measurements into the computational subspace.

The results are compared with those obtained from classical methods in Fig. 2,
where the quantum algorithm results are highly consistent with the classical
method results. Fig. 3 shows the dynamical evolution of the singlet and
triplet yields (when $\theta=\frac{\pi}{2}$) as simulated by the general
quantum algorithm on the IBM quantum simulator. After $2\times 10^{-4}s$, the
yields almost reach steady-state, which is consistent with the chosen decay
rate of $k=1\times 10^{4}s^{-1}$.

###  II.4 Complexity analysis

One factor that contributes to the complexity of the quantum algorithm is the
system’s size. For a density matrix of the size $n\times n$, the cost to
realize the unitary dilation of a most general $n\times n$ Kraus operator
$M_{k}$ is $\mathcal{O}(n^{2})$ [15, 16]. However, in our calculation of the
dynamics of the RPM, the Kraus operators each represents a single elementary
physical process and thus the $M_{k}$ matrices are often sparse with few non-
zero elements. This means the practical complexity scaling of implementing
each $M_{k}$ matrix on a quantum circuit can be greatly reduced to
$\mathcal{O}(\log^{2}n)$. Taking into account the total $K$ number of $M_{k}$
matrices to be simulated on the quantum circuit, the total complexity scaling
is $\mathcal{O}(K\log^{2}n)$ for our given system. It is worth noting that the
$K$ is determined on a case-by-case basis by the dynamical model and different
$M_{k}$ matrices can be evolved in parallel, therefore the scaling in $K$ is a
“soft” scaling that does not contribute to either the depth or the width of
each individual quantum circuit [15, 16]. Another contributing factor to the
quantum algorithm’s complexity is the number of time steps. In the most
general case, as can be seen from Eq.(5) to Eq.(7), taking $s$ steps requires
$K^{s}$ matrices to be evolved, which is an exponential scaling in the number
of time steps. However, fortunately the actual number of matrix terms to be
simulated can be greatly reduced once again due to the sparsity of the $M_{k}$
matrices. As mentioned above, the Kraus operators represent elementary
physical processes and thus the $M_{k}$ matrices are often sparse with very
few non-zero elements: this means that most matrix product terms in e.g.
Eq.(7) are zero matrices or matrices with negligible norms. The actual number
of matrix products we need to evaluate is determined on a case-by-case basis.
In the current simulation, in theory the total number of terms in $n$-th
iteration will be $9^{n}$ with nine Kraus operators
$\\{M_{k}|k=0,1,2,...,8\\}$. However, since the product of each pair of the
the decay operators $\\{P_{k}|k=1,2,3,...,8\\}$ is $0$, the product of each
pair of the Kraus operator $\\{M_{k}|k=1,2,3,...,8\\}$ is 0. Therefore, there
will be only 8 more terms when adding one more iteration. Thus, there will be
$8\times n+1$ terms in $n$-th iteration, reducing the terms significantly.
More simulation details are in the Supplementary Information.

##  III Conclusion

Based on the calculations and results, we can conclude that the general
quantum algorithm based on the Sz.-Nagy dilation can accurately simulate the
RPM dynamics described by the Lindblad master equation. As discussed above,
the RPM is an important theory that can explain the magneto-reception process
of the avian compass. The RPM also acts as an ideal model to help explaining
the isotope effects in xenon anaesthesia and lithium treatment of
hyperactivity, magnetic field effects on the circadian clock, as well as
hypomagnetic field effects on neurogenesis and microtubule assembly [37, 38].
Our simulation of the RPM model with the quantum algorithm not only
demonstrates the generality of the algorithm but also helps bridging the gap
between applying tools of quantum-information science to the investigation of
new areas of quantum biology.

###### Acknowledgements.

We acknowledge the financial support of the National Science Foundation under
award number 2124511, CCI Phase I: NSF Center for Quantum Dynamics on Modular
Quantum Devices (CQD-MQD). We also acknowledge the financial support from the
U.S. Department of Energy under Award No. DE-SC0019215. We acknowledge the use
of IBM Quantum services for this work. The views expressed are those of the
authors, and do not reflect the official policy or position of IBM or the IBM
Quantum team.

## References

  * Kais [2014] S. Kais, Introduction to quantum information and computation for chemistry, Quantum Information and Computation for Chemistry , 1 (2014).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Quantum chemistry in the age of quantum computing, Chemical reviews 119, 10856 (2019).
  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Quantum algorithms for quantum chemistry and quantum materials science, Chemical Reviews 120, 12685 (2020).
  * Head-Marsden  _et al._ [2020] K. Head-Marsden, J. Flick, C. J. Ciccarino, and P. Narang, Quantum information and algorithms for correlated quantum matter, Chemical Reviews 121, 3061 (2020).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Childs  _et al._ [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Toward the first quantum simulation with quantum speedup, Proceedings of the National Academy of Sciences 115, 9456 (2018).
  * Wang  _et al._ [2011] H. Wang, S. Ashhab, and F. Nori, Quantum algorithm for simulating the dynamics of an open quantum system, Physical Review A 83, 062317 (2011).
  * Wang  _et al._ [2013] D.-S. Wang, D. W. Berry, M. C. De Oliveira, and B. C. Sanders, Solovay-kitaev decomposition strategy for single-qubit channels, Physical review letters 111, 130504 (2013).
  * Wei  _et al._ [2016] S.-J. Wei, D. Ruan, and G.-L. Long, Duality quantum algorithm efficiently simulates open quantum systems, Scientific Reports 6, 1 (2016).
  * Di Candia  _et al._ [2015] R. Di Candia, J. S. Pedernales, A. Del Campo, E. Solano, and J. Casanova, Quantum simulation of dissipative processes without reservoir engineering, Scientific reports 5, 1 (2015).
  * Sweke  _et al._ [2015] R. Sweke, I. Sinayskiy, D. Bernard, and F. Petruccione, Universal simulation of markovian open quantum systems, Physical Review A 91, 062308 (2015).
  * García-Pérez  _et al._ [2020] G. García-Pérez, M. A. Rossi, and S. Maniscalco, Ibm q experience as a versatile experimental testbed for simulating open quantum systems, npj Quantum Information 6, 1 (2020).
  * Childs and Li [2016] A. M. Childs and T. Li, Efficient simulation of sparse markovian quantum dynamics, arXiv preprint arXiv:1611.05543 (2016).
  * Kliesch  _et al._ [2011] M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert, Dissipative quantum church-turing theorem, Physical review letters 107, 120501 (2011).
  * Hu  _et al._ [2020] Z. Hu, R. Xia, and S. Kais, A quantum algorithm for evolving open quantum dynamics on quantum computing devices, Scientific Reports 10, 3301 (2020).
  * Hu  _et al._ [2022] Z. Hu, K. Head-Marsden, D. A. Mazziotti, P. Narang, and S. Kais, A general quantum algorithm for open quantum dynamics demonstrated with the fenna-matthews-olson complex, Quantum 6, 726 (2022).
  * Wang  _et al._ [2022] Y. Wang, E. Mulvihill, Z. Hu, N. Lyu, S. Shivpuje, Y. Liu, M. B. Soley, E. Geva, V. S. Batista, and S. Kais, Simulation of open quantum system dynamics based on the generalized quantum master equation on quantum computing devices, arXiv preprint arXiv:2209.04956 (2022).
  * Head-Marsden  _et al._ [2021] K. Head-Marsden, S. Krastanov, D. A. Mazziotti, and P. Narang, Capturing non-markovian dynamics on near-term quantum computers, Physical Review Research 3, 013182 (2021).
  * Pauls  _et al._ [2013] J. A. Pauls, Y. Zhang, G. P. Berman, and S. Kais, Quantum coherence and entanglement in the avian compass, Physical Review E 87, 062704 (2013).
  * Zhang  _et al._ [2014a] Y. Zhang, G. P. Berman, and S. Kais, Sensitivity and entanglement in the avian chemical compass, Physical Review E 90, 042707 (2014a).
  * Zhang  _et al._ [2014b] Y. Zhang, G. P. Berman, and S. Kais, The radical pair mechanism and the avian chemical compass: Quantum coherence and entanglement, International Journal of Quantum Chemistry 115, 1327 (2014b).
  * Dodson  _et al._ [2013] C. A. Dodson, P. J. Hore, and M. I. Wallace, A radical sense of direction: Signalling and mechachanism in cryptochrome magnetoreception, Trends in Biochemical Sciences 38 (2013).
  * Wiltschko and Wiltschko [2005] W. Wiltschko and R. Wiltschko, Magnetic orientation and magnetoreception in birds and other animals, Journal of Comparative Physiology A 191, 675 (2005).
  * Wiltschko and Wiltschko [1972] W. Wiltschko and R. Wiltschko, Magnetic compass of european robins, Science 176, 62 (1972).
  * Wiltschko and Wiltschko [2010] R. Wiltschko and W. Wiltschko, Avian magnetic compass: its functional properties and physical basis, Current Zoology 56, 265 (2010).
  * Davila  _et al._ [2003] A. F. Davila, G. Fleissner, M. Winklhofer, and N. Petersen, A new model for a magnetorecepter in homing pigeons based on interacting clusters of super-paramagnetic magnetite, Physics and Chemistry of the Earth 28, 647 (2003).
  * Wiltschko  _et al._ [1993] W. Wiltschko, U. Munro, H. Ford, and R. Wiltschko, Red light disrupts magnetic orientation of migratory birds, Nature 364, 525 (1993).
  * Wiltschko and Wiltschko [1995] W. Wiltschko and R. Wiltschko, Migratory orientation of european robins is affected by the wavelength of light as well as by a magnetic pulse, Journal of Comparative Physiology A 177, 363 (1995).
  * Wiltschko and Wiltschko [1998] W. Wiltschko and R. Wiltschko, Pigeon homing: Effect of various wavelengths of light during displacement, Naturwissenschaften 85, 164 (1998).
  * Wiltschko and Wiltschko [1999] W. Wiltschko and R. Wiltschko, The effect of yellow and blue light on magnetic compass orientation in european robins, $erithacus$ $rubecula$, Journal of Comparative Physiology A 184, 295 (1999).
  * Wiltschko  _et al._ [2005] R. Wiltschko, T. Ritz, K. Stapput, P. Thalau, and W. Wiltschko, Two different types of light-dependent responses to magnetic fields in birds, Current Biology 15, 1518 (2005).
  * Stapput  _et al._ [2008] K. Stapput, P. Thalau, R. Wiltschko, and W. Wiltschko, Orientation of birds in total darkness, Current Biology 18, 602 (2008).
  * Wiltschko  _et al._ [2011] R. Wiltschko, S. Denzau, D. Gehring, P. Thalau, and W. Wiltschko, Magnetic orientation of migratory robins, $erithacus$ $rubecula$, under long-wavelength light, Journal of Experimental Biology 214, 3096 (2011).
  * Wiltschko and Wiltschko [1978] R. Wiltschko and W. Wiltschko, Evidence for the use of magnetic outward-journey information in homing pigeons, Naturwissenschaften 65, 112 (1978).
  * Gauger  _et al._ [2011] E. M. Gauger, E. Rieper, J. J. L. Morton, S. C. Benjamin, and V. Vedral, Sustained quantum coherence and entanglement in the avian compass, Phys. Rev. Lett. 106, 040503 (2011).
  * Levy and Shalit [2014] E. Levy and O. M. Shalit, Dilation theory in finite dimensions: the possible, the impossible and the unknown, Rocky Mountain Journal of Mathematics 44, 203 (2014).
  * Kominis [2015] I. K. Kominis, The radical-pair mechanism as a paradigm for the emerging science of quantum biology, Modern Physics Letters B 29, 1530013 (2015).
  * Zadeh-Haghighi and Simon [2022] H. Zadeh-Haghighi and C. Simon, Magnetic field effects in biology from the perspective of the radical pair mechanism, Journal of the Royal Society Interface 19, 20220325 (2022).

##  Appendix A Quantum Simulation Details

Here we give an example of the quantum circuit of the $M_{1}=\sqrt{k_{d}\delta
t}P_{1}$ at the first time step, $5\times 10^{-5}s$ to evolve the states.
After multiplied by the unitary matrix accounting for the coherent part as in
Eq. (10), $M_{1}$ becomes a $10\times 10$ sparse matrix shown in Eq. (16) and
the non-zero values are shown in Eq. (LABEL:m1_values).

$\begin{bmatrix}0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
m_{1}&m_{2}&m_{3}&m_{4}&m_{5}&m_{6}&m_{7}&m_{8}&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
\end{bmatrix}$ (16)

$\displaystyle m_{1}=-5.22\times 10^{-2}-4.05\times 10^{-1}i$ (17)
$\displaystyle m_{2}=-6.53\times 10^{-4}+1.23\times 10^{-2}i$ $\displaystyle
m_{3}=-2.93\times 10^{-1}+3.80\times 10^{-3}i$ $\displaystyle
m_{4}=-1.70\times 10^{-2}+1.38\times 10^{-4}i$ $\displaystyle m_{5}=2.93\times
10^{-1}+3.78\times 10^{-3}i$ $\displaystyle m_{6}=-7.39\times
10^{-1}-2.62\times 10^{-3}i$ $\displaystyle m_{7}=-5.22\times
10^{-3}+4.05\times 10^{-1}i$ $\displaystyle m_{8}=5.63\times
10^{-4}+5.37\times 10^{-3}i$

After we apply the unitary dilation described in Eq. (2) on the $M_{1}$, we
obtain a $20\times 20$ unitary matrix $U_{M_{1}}$. To simulate this operation
on quantum simulator, we have to use 5 qubits to cover the 20 dimensions.
Leveraging the Qiskit’s transpile function (qiskit.compliler.transpile), we
decomposed the resulted unitary operator $U_{M_{1}}$ to 2097 gates where a
portion of the circuit is shown below (Fig. 4) using the basis gates: ‘u3’,
‘cx’ and ‘rz’ as shown in Eq.(18) , on the backend of qasm_simulator. The
details of the decomposition of all the quantum circuits used are available on
reasonable request.

$\displaystyle\begin{aligned}
\hfil\displaystyle\begin{split}U_{3}(\theta,\phi,\lambda)=\begin{pmatrix}\cos\left(\frac{\theta}{2}\right)&-e^{i\lambda}\sin\left(\frac{\theta}{2}\right)\\\
e^{i\phi}\sin\left(\frac{\theta}{2}\right)&e^{i(\phi+\lambda)}\cos\left(\frac{\theta}{2}\right)\end{pmatrix}\end{split}\\\
\hfil\displaystyle\begin{split}CX&=I\otimes|0\rangle\langle
0|+X\otimes|1\rangle\langle 1|=\begin{pmatrix}1&0&0&0\\\ 0&0&0&1\\\ 0&0&1&0\\\
0&1&0&0\end{pmatrix}\end{split}\\\
\hfil\displaystyle\begin{split}Rz(\theta)=\begin{pmatrix}e^{-i\frac{\theta}{2}}&0\\\
0&e^{i\frac{\theta}{2}}\end{pmatrix}\end{split}\end{aligned}$ (18)

Figure 4: An example of a portion of the quantum gate sequence of a Kraus
operator. This is only a small portion of the circuit. The full circuit has
2097 gates and the details will be available from the authors on reasonable
request.


# ON FAST CHARGED PARTICLES SCATTERING ON A FLAT RELATIVISTIC BEAM OF CHARGED
PARTICLES IN APPROXIMATION OF CONTINUOUS POTENTIAL

N.F. Shul’ga∗, National Science Center ‘‘Kharkiv Institute of Physics and
Technology’’, Kharkiv, Ukraine  V.N. Karazin Kharkiv National University,
Kharkiv, Ukraine  V.D. Koriukina ∗∗, National Science Center ‘‘Kharkiv
Institute of Physics and Technology’’, Kharkiv, Ukraine  V.N. Karazin Kharkiv
National University, Kharkiv, Ukraine

###### Abstract

The differential scattering cross section for the problem of fast charged
particles motion near a flat relativistic beam of charged particles was
obtained. The problem is considered in the eikonal approximation in the
representation of the beam by a continuous potential.

## Introduction

††∗shulga<EMAIL_ADDRESS>(corresponding author)

Starting with the works of M.L. Ter-Mikaelyan (see the book [1] and references
therein), the effects associated with the scattering of fast charged particles
during interaction with various targets have been actively investigated. The
collision of relativistic beams of charged particles is of particular
interest, since this process often occurs in modern research on particle
accelerators. The scattering process study is not only interesting itself, but
also helps to examine other processes (for example, radiation) that are caused
by scattering. This is possible due to the fact that the radiation cross
section is proportional to the scattering cross section.

Using the quantum electrodynamics basis [2, 3], we introduced the approach [5,
6] for considering problems of particle scattering on various targets from a
single point of view both in the Born and in the eikonal approximations. This
approach enables us to obtain differential scattering cross sections on
targets of complex configuration in a relatively simple way, creating an
alternative to the method of simulating particle trajectories in matter. This
work shows how to extend this approach to consider the scattering problems on
moving relativistic targets without additional transformations of coordinate
system. The method demonstrated here involves the the idea of a continuous
potential suggested by E. Lindhard [4]. This investigation is interesting from
a theoretical point of view and is also useful for the theoretical description
of modern beam collision experiments held at particle accelerators, in
particular at the SuperKEKB accelerator (Japan).

##  1 The Scattering Problem Formulation

Let us consider the problem of fast charged particles scattering on a beam of
charged relativistic particles in the eikonal approximation of quantum
electrodynamics. For simplification, we will study the interaction of one
incident particle with the beam. We choose beam properties close to the
properties available at modern experimental facilities, namely at the
SuperKEKB accelerator (Japan).

Figure 1: The motion of a fast charged particle incident parallel to the flat
beam

Let the particle incident on the beam parallel to the beam velocity, along the
$z$-axis. The potential of a relativistic beam is determined by the four-
potential

$\displaystyle
A_{\mu}^{(N)}=\left\\{1,0,0,-\frac{v_{b}}{c}\right\\}A_{0}^{(N)},$ (1)

where the beam velocity has the following magnitude and direction:
$\vec{v}_{b}=-v_{b}\vec{e}_{z}$. The total scalar potential $A_{0}^{(N)}$ of
the target, which consists of $N$ particles of the same type (electrons or
positrons), is the sum of the scalar potentials of individual target particles

$\displaystyle A_{0}^{(N)}=\sum_{n=1}^{N}A_{0,n}.$ (2)

The scalar potential of each individual target particle, in its turn, is the
relativistic Coulomb potential

$\displaystyle
A_{0,n}=\frac{e}{\sqrt{(z+v_{b}t)^{2}+(\vec{\rho}-\vec{\rho}_{n})^{2}/\gamma_{b}^{2}}},$
(3)

where $\vec{\rho}=(x,y)$ are the coordinates transverse to the initial
direction of particle movement, $\vec{\rho}_{n}=(x_{n},y_{n})$ is the position
of each individual particle of the beam relative to the center of the beam. We
note that we neglect the deflection of the target due to the interaction with
the incident particle, so these coordinates coincide both in the coordinate
system associated with the beam and in the laboratory coordinate system
throughout the scattering process.

The particle distribution function in the target for simplification will be
considered as follows: along the $x$-axis, the beam particles are distributed
according to the Gaussian law with a random deviation $x_{n}$ and variance
$\langle u^{2}_{x}\rangle$; along the $y$ axis, the beam has a zero thickness,
which is equivalent to a zero variance $\langle u^{2}_{y}\rangle$:

$\displaystyle\vec{\rho}_{n}:\begin{cases}f(x_{n})=\frac{1}{\sqrt{2\pi\langle
u^{2}_{x}\rangle}}\>exp\left[\frac{-x^{2}_{n}}{2\langle
u^{2}_{x}\rangle}\right],\\\ y_{n}=0,\>\langle u^{2}_{y}\rangle=0.\end{cases}$
(4)

##  2 Wave Function of Scattered Particle

From the Klein-Gordon and Dirac equations, we find in the eikonal
approximation the wave function for a particle that scatters on a relativistic
beam. The Klein-Gordon equation for a particle of mass $m$, four-momentum
$p_{\mu}$, moving in a four-potential $A_{\mu}^{(N)}$ is

$\displaystyle\left[\left(p_{\mu}-eA_{\mu}^{(N)}\right)^{2}-m^{2}\right]\Phi=0.$
(5)

The wave function $\Phi$ in the eikonal approximation is

$\displaystyle\Phi=f\>e^{iS}.$ (6)

The amplitude $f$ and the phase $S$ of this wave function can be expanded into
a series of the reciprocal powers of the momentum $f_{k}\sim
p^{-k},\chi_{k}^{(N)}\sim p^{-k}$ respectively:

$\displaystyle f=f_{0}+f_{1}(\vec{r},t)+...\ ,$ (7) $\displaystyle
S=\vec{p}\vec{r}-Et+\chi_{0}^{(N)}(\vec{r},t)+\chi_{1}^{(N)}(\vec{r},t)+...\
.$ (8)

Let us determine $S$ through the relativistic Hamilton-Jacobi equation

$\displaystyle\left(\partial_{t}S+eA_{0}^{(N)}\right)^{2}-\left(\nabla
S-e\vec{A}^{(N)}\right)^{2}-m^{2}=0.$ (9)

It leads to the equation

$\displaystyle
c\partial_{z}\chi_{0}^{(N)}+\partial_{t}\chi_{0}^{(N)}+2eA_{0}^{(N)}=0,$ (10)

from which the $\chi_{0}$-function for scattering on a target of $N$ particles
has the following form

$\displaystyle\chi_{0}^{(N)}=-\frac{1}{c}\int_{-\infty}^{\xi}d\xi
eA_{0}^{(N)}(\vec{\rho},\xi),$ (11)

where $\xi=z+v_{b}t$.

So from the Dirac equation, we obtaine the solution, which satisfies equations
(5) and (9):

$\displaystyle f_{0}=u,$ (12)

where $u$ is a bispinor.

Hence, the wave function for the given case in the zero eikonal approximation
is

$\displaystyle\Phi=ue^{i\left[\vec{p}\vec{r}-Et+\chi_{0}^{(N)}(\vec{\rho},z+v_{b}t)\right]}.$
(13)

##  3 Scattering Cross Section in the Eikonal Approximation

Knowing the general form of the wave function for the scattering of a fast
charged particle on a relativistic beam of charged particles, we will find the
general form of the differential scattering cross section for this problem. In
accordance with quantum electrodynamics, the differential scattering cross
section for a particle, which from the initial state with energy
$\varepsilon^{\prime}$ and momentum $\vec{p}\>^{\prime}$ transits to the final
state with energy $\varepsilon$ and momentum $\vec{p}$, is defined as [2]

$\displaystyle
d\sigma=\frac{1}{j}\frac{w_{fi}}{t}\frac{V\>d^{3}p}{(2\pi)^{3}},$ (14)

where $w_{fi}$ is the probability of a particle transition from the initial
state to the final state, $t$ is the observation time of the scattering
process, $V$ is the normalizing volume, $j$ is the flow density of interacting
particles, which is expressed as

$\displaystyle j=\frac{2v}{V},$ (15)

where $v$ is the interacting particles velocity, assuming that the particles
move with the same velocity magnitudes. The probability of the particle
transition from the initial state to the final state is given as follows

$\displaystyle w_{fi}=\frac{\left|T_{fi}\right|^{2}}{(2\varepsilon
V)(2\varepsilon^{\prime}V)},$ (16) $\displaystyle T_{fi}=\int
dtd^{3}r\>e^{i(\varepsilon^{\prime}t-\vec{p}\>^{\prime}\vec{r})}\bar{u}^{\prime}\>(\gamma_{0}+\gamma_{3}v)eA_{0}^{(N)}(\vec{\rho},z+vt)\Phi(\vec{r},t),$
(17)

where $\bar{u}^{\prime}$ is the bispinor of the initial state, $\gamma_{0}$
and $\gamma_{3}$ are the Dirac matrices, $\Phi$ is the wave function of the
final state in the zero eikonal approximation (13).

Considering that for ultrarelativistic particles $\varepsilon\rightarrow pc$,
and also making the substitution $z+vt\rightarrow z$, we obtain

$\displaystyle T_{fi}\approx 2\bar{u}^{\prime}\gamma_{0}u\int
dtd^{3}r\>e^{i(\varepsilon^{\prime}-\varepsilon)t+i\vec{q}\vec{r}}eA_{0}^{(N)}(\vec{\rho},z+vt)\>e^{i\chi_{0}(\vec{\rho},z+vt)}=$
(18) $\displaystyle=2\bar{u}^{\prime}\gamma_{0}u\int
dtdzd^{2}\rho\>e^{i(\varepsilon^{\prime}-\varepsilon-
q_{z}v)t+i\vec{q}\vec{r}}eA_{0}^{(N)}(\vec{\rho},z)\>e^{i\chi_{0}(\vec{\rho},z)}$

where $\vec{q}=\vec{p}-\vec{p}\>^{\prime}$ is momentum transfer.

##  4 Discussion of the Violation of the Energy Conservation Law

Integrating (18) over $t$ leads to the emergence of the Dirac delta function
$\delta(\varepsilon^{\prime}-\varepsilon-q_{z}v)$, which in the stationary
target case is $\delta(\varepsilon^{\prime}-\varepsilon)$ and corresponds to
the energy conservation law. The "new" energy conservation law has no physical
sense. Let us consider this expression in detail.

According to the energy-momentum conservation laws, the momentum transferred
along the motion direction is expressed with some accuracy as

$\displaystyle
q_{z}\approx\frac{\varepsilon^{2}-\varepsilon^{{}^{\prime}2}+\vec{q}^{2}_{\perp}}{2p},$
(19)

where $\vec{q}_{\perp}$ is momentum transfer in perpendicular to the movement
direction. Then

$\displaystyle\varepsilon^{\prime}-\varepsilon-q_{z}v\approx
2(\varepsilon-\varepsilon^{\prime}).$ (20)

As a result, after substituting (20), the formula (18) becomes

$\displaystyle T_{fi}\approx 2\bar{u}^{\prime}\gamma_{0}u\int
dtd^{3}r\>e^{2i(\varepsilon^{\prime}-\varepsilon)t+i\vec{q}\vec{r}}eA_{0}^{(N)}(\vec{\rho},z)\>e^{i\chi_{0}(\vec{\rho},z)},$
(21)

which corresponds to $T_{fi}$ in the zero eikonal approximation. In this from,
the energy conservation law is fulfilled.

Beyond the zero eikonal approximation, another way to "restore" the energy
conservation law can be demonstrated. In the first eikonal approximation

$\displaystyle T_{fi}\approx 2\bar{u}^{\prime}\gamma_{0}u\int
dtd^{3}r\>e^{i(\varepsilon^{\prime}-\varepsilon-
q_{z}v)t+i\vec{q}\vec{r}}eA_{0}^{(N)}(\vec{\rho},z)\>e^{i\chi_{0}(\vec{\rho},z)+i\chi_{1}(\vec{\rho},z)},$
(22)

where $\chi_{1}(\vec{\rho},z)$ from the Hamilton-Jacobi equation (9) is

$\displaystyle\chi_{1}=-\frac{1}{4p}\int^{z+vt}_{-\infty}dz\
\left(\partial_{\vec{\rho}}\chi_{0}\right)^{2}.$ (23)

At a certain choice of the coordinate-time system: $z=vt$, then the
integration with given limits can be replaced by integration over the region
where the particle interacts with the beam, that is, from $z+vt=0$ to
$z+vt=2L_{z}$, where $L_{z}$ is the longitudinal size of the beam. Then, using
the method of integration by parts, we get

$\displaystyle\chi_{1}=-\frac{z\left(\partial_{\vec{\rho}}\chi_{0}\right)^{2}\Big{|}^{z=2L_{z}}_{z=0}}{4p}+\frac{1}{4p}\int^{z+vt}_{-\infty}d\left(\partial_{\vec{\rho}}\chi_{0}\right)^{2}\
z.$ (24)

Considering that in accordance with the quasi-classical approximation
$\partial_{\vec{\rho}}\chi_{0}=-\vec{q}_{\perp}$, the first term in (24)
becomes $\frac{\vec{q}_{\perp}^{2}L_{z}}{2p}$ and is equivalent to the term
$q_{z}vt$. It follows that the term $q_{z}vt$ vanishes due to the change in
the definition of the $\chi_{1}(\vec{\rho},z)$-function and the energy
conservation law is fulfilled.

Integrating (21) over $t$ and $z$, we obtaine the formula for the differential
scattering cross section, which coincides with the corresponding formula for
scattering on a stationary target

$\displaystyle\frac{d^{2}\sigma}{dq_{x}dq_{y}}=|a(\vec{q}_{\perp})|^{2},$ (25)

where

$\displaystyle a(\vec{q}_{\perp})=\int_{-\infty}^{\infty}d^{2}\rho
e^{\frac{i}{\hbar}\vec{q}\vec{\rho}}\left\\{1-exp\left[\frac{i}{\hbar}{\chi}_{0}^{(N)}\left(\vec{\rho}\right)\right]\right\\}.$
(26)

Using the method described in [6] with accuracy up to linear in
$\bar{\chi}^{(1)}_{0}$ terms ($\bar{\chi}^{(1)}_{0}$ is ${\chi}^{(1)}_{0}$
averaged over the particles positions in the target and reduced to one
particle of the target so that ${\chi}^{(N)}_{0}=N{\chi}^{(1)}_{0}$) we obtain
the formula for the differential scattering cross section

$\displaystyle\langle\frac{d^{2}\sigma}{dq_{\perp}^{2}}\rangle\approx\frac{1}{4\pi^{2}}\left|\int
d^{2}\rho
e^{\frac{i}{\hbar}\left[\vec{q}_{\perp}\vec{\rho}+N\bar{\chi}^{(1)}_{0}\left(\vec{\rho}\right)\right]}\right|^{2}.$
(27)

The coincidence of differential cross sections for relativistic and stationary
targets is justified, since the given cross section is invariant. We note that
neglecting quadratic in $\bar{\chi}^{(1)}_{0}$ terms simplifies the problem
consideration and corresponds to the continuous potential approximation.

##  5 The Differential Cross Section for a Particle Scattering on a Flat Beam

Let us obtain the necessary functions for the scattering cross section
calculation, namely $\bar{\chi}_{0}^{(1)}$ and its derivatives. Averaging
${\chi}_{0}^{(1)}$ using the particle distribution function in the target (4)
gives

$\displaystyle\bar{\chi}_{0}^{(1)}=\frac{\alpha}{2}\left\\{\frac{(y+ix)^{2}}{2\langle{u_{x}^{2}}\rangle}\>_{2}F_{2}\left[1,1;\frac{3}{2},2;\frac{(y+ix)^{2}}{2\langle{u_{x}^{2}}\rangle}\right]-\frac{\pi}{2}Erfi\left[\frac{y+ix}{\sqrt{2\langle{u_{x}^{2}}\rangle}}\right]+\right.$
(28)
$\displaystyle\left.+\frac{(y-ix)^{2}}{2\langle{u_{x}^{2}}\rangle}\>_{2}F_{2}\left[1,1;\frac{3}{2},2;\frac{(y-ix)^{2}}{2\langle{u_{x}^{2}}\rangle}\right]-\frac{\pi}{2}Erfi\left[\frac{y-ix}{\sqrt{2\langle{u_{x}^{2}}\rangle}}\right]\right\\}+const$

where $\>{}_{2}F_{2}$ is generalized hypergeometric function, $Erfi$ is
imaginary error function, $const$ is determined by the cutoff which we need to
make, since the long-range potential of the beam leads to the divergence of
the $\bar{\chi}_{0}^{(1)}$-function .

Figure 2: The $\bar{\chi}_{0}^{(1)}$-function for the flat beam

(a) $\partial_{x}\bar{\chi}_{0}^{(1)}$

(b) $\partial_{y}\bar{\chi}_{0}^{(1)}$

Figure 3: Derivatives of the $\bar{\chi}_{0}^{(1)}$-function for the flat beam

(a) $q_{x}=-0.5\ q_{x\ max},\ q_{y}=-0.18\ q_{y\ max}$

(b) $q_{x}=-0.9\ q_{x\ max},\ q_{y}=-0.12\ q_{y\ max}$

Figure 4: Graphic solution of (32) for specific values of momentum transfer
$\left(q_{x},q_{y}\right)$

(a) $q_{y}=0.001\ q_{y\ max}$

(b) $q_{y}=0.3\ q_{y\ max}$

(c) $q_{y}=0.56\ q_{y\ max}$

(d) $q_{y}=0.89\ q_{y\ max}$

Figure 5: The differential cross-section for scattering on the flat beam for
$N=100$ particles in the beam and various $q_{y}$ values

The first derivatives for $x$ and $y$ are determined respectively

$\displaystyle\partial_{x}\bar{\chi}_{0}^{(1)}=\alpha\>\sqrt{\frac{2\pi}{\langle
u^{2}_{x}\rangle}}\>Im\left\\{w\left[\frac{-x+i|y|}{\sqrt{2\langle
u^{2}_{x}\rangle}}\right]\right\\},$ (29)

$\displaystyle\partial_{y}\bar{\chi}_{0}^{(1)}=-\alpha\sqrt{\frac{{\pi}}{{2\langle
u^{2}_{x}\rangle}}}sgn(y)Re\left\\{w\left[\frac{-x+i|y|}{\sqrt{2\langle
u^{2}_{x}\rangle}}\right]\right\\},$ (30)

where $w[z]$ is Faddeeva function ( $w[z]=e^{-z^{2}}Erfc\left[-iz\right]$,
where $Erfc$ is complementary error function).

According to the stationary phase method, (27) can be approximately written as

$\displaystyle\langle\frac{d^{2}\sigma}{dq_{x}dq_{y}}\rangle=\left.\frac{1}{N^{2}\left|\partial_{x}^{2}\bar{\chi}_{0}^{(1)}\partial_{y}^{2}\bar{\chi}_{0}^{(1)}-\left(\partial_{xy}\bar{\chi}_{0}^{(1)}\right)^{2}\right|}\right|_{x=\tilde{x},y=\tilde{y}},$
(31)

where the points $\tilde{x},\tilde{y}$, which correspond to the phase minimum
condition, are given by the system of equations

$\displaystyle\begin{cases}q_{x}+N\partial_{x}\bar{\chi}_{0}^{(1)}=0,\\\
q_{y}+N\partial_{y}\bar{\chi}_{0}^{(1)}=0.\end{cases}$ (32)

The given system of equations demonstrates that according to the quasi-
classical approximation, the first derivatives of the
$\bar{\chi}_{0}^{(1)}$-function determine the momentums transferred along the
corresponding coordinate axes (we denote maximum absolute values of
$q_{x},q_{y}$ appropriate for this system as $q_{x\ max},q_{y\ max}$). We will
solve this system numerically. Let us illustrate the solution of this system
graphically with the figure 4: the solution is the point of intersection of
the curves which correspond to the fixed level lines of
$\partial_{x}\bar{\chi}_{0}^{(1)},\partial_{y}\bar{\chi}_{0}^{(1)}$. We can
see that for some values of the momentum transfer the system of equations has
no solution.

Let us obtain results using the formula (31) for the region of momentum
transfer, for which the system of equations (32) has a solution. We will also
calculate the differential scattering cross section by numerical integration
of the formula (26). Differential cross sections obtained by these methods are
presented in figures 5-7.

(a) $q_{y}=0.001\ q_{y\ max}$

(b) $q_{y}=0.3\ q_{y\ max}$

(c) $q_{y}=0.78\ q_{y\ max}$

(d) $q_{y}=0.89\ q_{y\ max}$

Figure 6: The differential cross-section for scattering on the flat beam for
$N=1\ 000$ particles in the beam and various $q_{y}$ values

(a) $q_{y}=0.001\ q_{y\ max}$

(b) $q_{y}=0.3\ q_{y\ max}$

(c) $q_{y}=0.56\ q_{y\ max}$

(d) $q_{y}=0.89\ q_{y\ max}$

Figure 7: The differential cross-section for scattering on the flat beam for
$N=10\ 000$ particles in the beam and various $q_{y}$ values

##  6 Discussion

So, we obtained the differential cross section for fast charged particles
scattering on a flat beam of relativistic charged particles in the eikonal
approximation of quantum electrodynamics. The eikonal approximation, compared
to the Born approximation, has a wider application region, which justifies its
use. Also, to simplify the calculations, we used the continuous potential
approximation, which followed naturally from the approach suggested in [5, 6].
This approximation does not take into account incoherent scattering, but
reflects scattering in general terms.

We note that due to the long-range potential, the functions determining the
scattering in our case have different properties compared to the functions
associated with the short-range potential (which, for example, occurs in a
crystal). The $\bar{\chi}_{0}^{(1)}$-function has a logarithmic divergence
which is usual for the Coulomb potential, due to which we needed to implement
the cutoff for physical reasons. The cutoff is justified by the fact that in
real experimental facilities the beam potential is shielded at certain
distances either by oppositely charged particles of the incident beam or by
the walls of the facility. The $\partial_{x}\bar{\chi}_{0}^{(1)}$ function,
which according to the quasi-classical approach determines the momentum
transferred along the $x$-axis, is proportional to $x^{-1}$ at infinity, which
means that it decreases much slowlier than a similar function for scattering
in a crystal. Such behavior of the functions determines the significant
influence of the long-range impact parameters on the scattering. It is unlike
the scattering in short-range potentials, where the long-range impact
parameters almost do not contribute to the scattering cross section.

We also pay attention to the notion that the oscillations of the differential
cross section obtained by direct numerical integration can be explained by the
cutoff that we made and might have no physical sense. This issue requires
additional consideration.

The manifestation of the rainbow scattering effect for this problem also needs
additional research, but we assume that for fixed transferred momentum $q_{y}$
and fixed impact parameters in the $y$-axis direction, the rainbow scattering
is possible. This assumption is based on the form of the
$\partial_{x}\bar{\chi}_{0}^{(1)}$-function (Fig. 8), which determines the
transferred momentum $q_{x}$ according to a quasi-classical approach.

Figure 8: The $\partial_{x}\bar{\chi}_{0}$-function for the flat beam for
various $y$ values

## Acknowledgements

The work was partially supported by the National Academy of Sciences of
Ukraine (project 0123U103077) and by Agence Universitaire de la Francophonie.
We also acknowledge the fruitful discussions with I.V. Kyryllin, S.V.
Trofymenko, M.V. Bondarenko, I.V. Omelchenko.

## References

  * [1] M.L. Ter-Mikaelian, High-Energy Electromagnetic Processes in Condenced Media, Wiley, New York, 1972
  * [2] A.I. Akhiezer, V.B. Berestetsky. Quantum Electrodynamics. Interscience Publishers. 1965
  * [3] A.I. Akhiezer, N.F. Shul’ga, High-Energy Electrodynamics in Matter. Amsterdam, Gordon and Breach Pub., 1996
  * [4] J. Lindhard, K. Dan. Vidensk. Selsk. Mat.Fys. Medd. 31 (1965) 14
  * [5] N.F. Shul’ga, V.D. Koriukina. On coherent and incoherent scattering of fast charged particles in ultrathin crystals. Problems of Atomic Science and Technology. 2020. V. 127. pp. 120-125. https://doi.org/10.46813/2020-127-120
  * [6] Shul’ga N.F., Koriukina V.D. The Eikonal Approximation of the Scattering Theory for Fast Charged Particles in a Thin Layer of Crystalline and Amorphous Media. Nucl. Instr. Meth B. 2021. Vol. 487. P. 25-29. https://doi.org/10.1016/j.nimb.2020.09.014


11institutetext: School of Computer Science and Engineering, Kyungpook
National University, Daegu, Republic of Korea,
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Department of Computer Science, Information Technology
University (ITU), Lahore, Pakistan,
22email<EMAIL_ADDRESS>33institutetext: Department of Computer
Science, University of Warwick, United Kingdom,
33email<EMAIL_ADDRESS>
# Unsupervised RGBD Video Object Segmentation Using GANs

Maryam Sultana 11 0000-0002-8831-843X Arif Mahmood 22 0000-0001-5986-9876
Sajid Javed 33 Soon Ki Jung 11 0000-0003-0239-6785

###### Abstract

Video object segmentation is a fundamental step in many advanced vision
applications. Most existing algorithms are based on handcrafted features such
as HOG, super-pixel segmentation or texture-based techniques, while recently
deep features have been found to be more efficient. Existing algorithms
observe performance degradation in the presence of challenges such as
illumination variations, shadows, and color camouflage. To handle these
challenges we propose a fusion based moving object segmentation algorithm
which exploits color as well as depth information using GAN to achieve more
accuracy. Our goal is to segment moving objects in the presence of challenging
background scenes, in real environments. To address this problem, GAN is
trained in an unsupervised manner on color and depth information independently
with challenging video sequences. During testing, the trained GAN generates
backgrounds similar to that in the test sample. The generated background
samples are then compared with the test sample to segment moving objects. The
final result is computed by fusion of object boundaries in both modalities,
RGB and the depth. The comparison of our proposed algorithm with five state-
of-the-art methods on publicly available dataset has shown the strength of our
algorithm for moving object segmentation in videos in the presence of
challenging real scenarios.

###### Keywords:

Video object segmentation Generative Adversarial Networks Background
Estimation.

##  1 Introduction

A fundamental step in many computer vision and artificial intelligence
applications involves objects segmentation, for various tasks such as object
detection [27, 1], visual object tracking [30, 26], video surveillance [3,
11], salient motion detection [5] and image inpainting [25]. For object
segmentation, background modeling plays a crucial role as it is a key process,
which describes a scene without the presence of any foreground objects.
However, foreground object detection is the process of extracting moving
objects with prior knowledge of the background. Foreground object detection
becomes more challenging in real-time environments for instance illumination
variations, out of range scenarios and color camouflage of foreground objects
concerning the background information [17]. There are other challenging
conditions as well, such as sequences with shadows or ghosting artifacts and
also bootstrapping in which foreground objects are present in almost all
frames of a video sequence. Over the past decade, many techniques have been
proposed in the literature to address the problems of these challenging
background scenes for the tasks of foreground detection and evaluation [2, 12,
23].

To address the challenges mentioned above for object segmentation, we present
a fusion based deep learning method. Our proposed method is based on
Generative Adversarial Network (GAN) [8] working on the idea of back
propagation steps to generate specific kind of data. In this study, our
primary focus is foreground object segmentation in the presence of various
challenging conditions in background scenes. To address this problem, we
present a solution based on GANs [22] which works on the principle of
generating background samples with scene specific information. Our proposed
GAN model trains in an unsupervised manner on all video sequences irrespective
of the presence of foreground objects in them as scene specific model. The
primary purpose of training scene-specific GAN model is that the network
should be able to learn the semantics of the scene containing various
foreground background objects in the presence of challenging conditions such
as illumination variations, shadows, and dynamic background information. The
key idea behind training our GAN model with various challenging scenes is that
during testing, our model will be able to generate the background scene
according to our given test sample information. Since our network is trained
on the data containing both background as well as foreground information, we
need to eliminate moving foreground objects from our test sample in order to
generate the exact background image via our trained GAN model. This can be
done by multiplying test samples containing foreground objects with their
motion masks evaluated by using optical flow [15]. This step helps our GAN
model to generate exact image sample containing only background information
similar in semantics as test sample via back-propagation technique. The
generated background sample is then subtracted from the given test sample to
detect foreground objects as shown in Figure 1.

##  2 Related Work

In the last two decades, several methods have proposed for foreground
detection by exploiting color and depth information. There are various
traditional background subtraction methods, for instance, a very popular and
classic method for background subtraction is GMM [24]. The idea of this
proposed method is, modeling of each pixel is done with a mixture of Gaussian,
and it works on color information processing only. However, SOBS algorithms,
such as [18] and [16], shows improvement in performance by using color and
depth features. They are based on self-organizing neural networks, which
achieves better results in various challenging environments. Another efficient
background subtraction model with a fusion of depth information embedded in
the basic structure of Robust Principal Component Analysis (RPCA) is proposed
by Javed et al. [9, 10]. It is a hybrid RPCA model with spatial and temporal
information handling mechanism. Therefore, it is called Spatiotemporal Robust
Principal Component Analysis (SRPCA). SRPCA algorithm is a graph-regularized
method, which preserves the spatiotemporal information of background that is
low-rank matrix formation in the form of dual spectral graphs. To address the
challenges in background subtraction problem, Bo Xin et al. [28] presented a
technique called Background Subtraction via Generalized Fused Lasso Foreground
Modeling (BS-GFL). Their primary objective is to address the problem of
missing information in foreground detection because of various challenges like
illumination variations and dynamic backgrounds. To solve this problem, they
consider generalized fused lasso regularization to search for intact
structured foregrounds to recover the missing content during the background
subtraction process. Although BS-GFL is an efficient algorithm with a very
good performance, however, it is an offline and partially supervised method.
To solve the same problem of the missing content of foreground regions during
background subtraction process, Xiaowei Zhou et al. [31] proposed a method
called Moving Object Detection by Detecting Contiguous Outliers in the Low-
Rank Representation (DECOLOR). The authors of this proposed method exploited
Markov random field technique to recover the missing content of foreground
detection during background subtraction process. For foreground detection by
exploiting color as well as depth information, Massimo De Gregorio et al. [6]
presented an algorithm known as cwisardH+ working on the idea of decoupling
the color information from the pixel depth information. The two video
sequences runs are synchronously in this algorithm but independently and
modeled by weightless neural networks at each pixel value.

##  3 Proposed Method

In this section, we describe each step of our proposed algorithm in detail.
The workflow diagram of our proposed method ForeGAN_RGBD is presented in
Figure 1. Our proposed method aims to perform foreground object segmentation
by adding RGB and depth information into DCGAN model [21]. ForeGAN_RGBD has
two phases: Phase 1.) Training of the RGB and depth video sequences with two
independent models as shown in Figure 1. Phase 2.) Testing of the two trained
models with video sequences including foreground objects or background objects
via backpropagation technique. It means that in phase 1 we learn two models
based on GAN [29] representing indoor and outdoor scenes with various
challenges on color as well as depth information. Each GAN train two models
simultaneously, a discriminator model and a generator model, to differentiate
between real and fake generated data containing RGB or depth information. The
objective of the model is to achieve equilibrium of costs, increasing the
ability of the model to generate data that is more accurate. In the following,
we explain our proposed ForeGAN_RGBD model in detail.

Figure 1: Workflow of the proposed ForeGAN_RGBD model. Phase 1: Training has
two independent models one to learn color information of different scenes,
called ForeGAN_RGB model and the other to learn depth information which is
called ForeGAN_D model. Phase 2: Testing also has two independent trained
models for evaluation. (1) Input test RGB image sample, (2) Motion masks of
test sample via optical flow, (3) Input RGB image sample multiplied by motion
mask, (4) Generated RGB background image sample, (5) RGB foreground detection
by equation 3, (6) Input test depth image sample, (7) Generated depth
background image sample, (8) Depth foreground detection also by equation 3,
(9) Motion mask, (10) Depth foreground detection after masking to eliminate
noisy background pixels, and (11) Final output by fusion of results from (10)
and (5) which is simply pixel by pixel addition. The samples presented in this
Figure are from SBM-RGBD dataset [4] video sequence ’genSeq2’ category
’Shadows’.

###  3.1 ForeGAN_RGBD Phase 1

A GAN model has two adversarial modules, a discriminator $D$ and a generator
$G$. The main objective of generator $G$ is to learn a distribution $p_{gen}$
over input data $X_{i}$ via mapping of $z$ samples through $G(z)$. This
mapping facilitates the $1D$ vectors of input noise which is uniformly
distributed and sampled from latent space $Z$ to the $2D$ image
representation. In a basic GAN model discriminator, $D$ is a Convolutional
Neural Network (CNN) model that maps a $2D$ image representation to a single
value $D(\cdot)$. This single value $D(\cdot)$ of discriminator’s output can
be explicated as the probability that whether the input given to the
discriminator $D$ was a fake generated image $G(z)$ by the generator $G$ or a
real image $X$ sampled from training data $X_{i}$. The discriminator and the
generator are simultaneously optimized via cross entropy loss functions in a
following two-player minimax game with $\Gamma(D,G)$ as a value function:

$\begin{split}\min\limits_{G}\max\limits_{D}~{}\Gamma(G,D)=\mathbb{E}_{x\sim
p_{data}(x)}[log(D(x))]\\\ +\mathbb{E}_{z\sim
p_{z}(z)}[log(1-D(G(z)))].\end{split}$ (1)

The discriminator in GAN model is a decision maker entity which is trained to
maximize the probability of assigning real training sample to actual input
data and samples from $p_{gen}$ to the fake generated data. The purpose of the
generator $G$ is to fool $D$ via minimizing the cost function $\Gamma(G)$ =
$log(1-D(G(z)))$, which is basically to maximize the following equation:

$\Gamma(G)=D(G(z)).$ (2)

During the training process, the generator tries to improve itself by
generating realistic images, and the discriminator tries to identify the real
and fake generated images. Once the training is done then the next step is
testing which is discussed in detail in the next section.

###  3.2 Inverse Mapping of Test Images to Latent Space Representation, Phase
2

During the phase of adversarial training, the generator learns the mapping
from latent space representations $z$ which is random noise to more realistic
images, $G(z)=z\mapsto x$. However inverse mapping $\theta(x)=x\mapsto z$ in
GAN is not a straightforward process, rather we need a different mechanism for
this purpose. To achieve inverse mapping Yeh et al. [29] presented a back-
propagation method to input data. This technique was already used by Google’s
DeepDream to create dreamlike images [20]. Moreover, the back-propagation
method has also been used to understand and visualize neural network’s learned
features by inverting the network by updating gradients at the level of input
layer [29], [19], [14], [7]. Every back-propagation based method requires
specific loss functions for inverse mapping process in a neural network.
Therefore our method has two loss functions as well. The purpose of the
inverse mapping in our proposed method is to generate the specific kind of
data during testing via back propagation strategy.

####  3.2.1 Foreground Object Segmentation Loss

Given a test image $x$ we aim to find that particular random noise $z$ in the
latent space that was mapped to generate image $G(z)$ via back-propagation
method. To find that specific $z$, we have to select an initial random sample
$z_{o}$, from the latent space and reinforce it to the trained generator
network to generate $G(z_{o})$. The loss functions are designed by the
generated image $G(z_{o})$, which provides significant update information of
the coefficients of $z_{o}$ leading our back-propagation step to be shifted
from $z_{o}$ to $z_{1}$ in the latent space distribution. The most similar
generated image $G(z_{\beta})$ can be found by various back-propagation steps
$\beta=0,1,2,...,\phi$ by minimizing the following loss function:

$\Upsilon_{F}(z_{\beta})=\sum|x-G(z_{\beta})|,$ (3)

$\Upsilon_{F}(z_{\beta})$ function measures the visual dissimilarity between
the test image $x$ and generated image $G(z_{\beta})$ via various back-
propagation steps. In an ideal scenario if test and generated images are
identical, $\Upsilon_{F}(z_{\beta})$ will be zero.

####  3.2.2 Feature Matching Loss

Now the next issue is how to generate those images $G(z_{\beta})$ which
exactly makes the best match with test image $x$. To solve this problem, we
need another loss term which helps the generator to generate similar images as
test images and minimize the loss function in equation $\eqref{FDL}$. In
contrast to the loss functions for best matching technique defined in [29],
Thomas Schlegl et al. [22] adopted an interesting method to lead generator to
generate the specific images as test images. This method addresses the problem
arises due to over-training of the discriminator, which causes the instability
of GANs. To improve the inverse mapping of test image $x$ to that specific
random noise $z$, feature matching technique is designed to force the
generator to generate the data with similar statistics as test data. Since
discriminator feeds its gradients to the generator, it is best to design the
loss function on discriminator. This loss function is defined in such a way
that the intermediate feature layer of the discriminator is feed with the
generated image $G(z_{\beta})$:

$\Upsilon_{M}(z_{\beta})=\sum|l(x)-l(G(z_{\beta}))|,$ (4)

Where $l(\cdot)$ represents the output of the intermediate layer of the
discriminator, which describes the test image $x$. Based on this loss
function, the discriminator is now used as a feature extractor rather than a
decision maker for real or fake image representations during testing process
only.
The overall loss functions can be represented as the weighted sum of both loss
terms defined in equations (3) and (4):

$\Upsilon_{M}(z_{\beta})=(1-\eta)\Upsilon_{F}(z_{\beta})+\eta\Upsilon_{M}(z_{\beta}).$
(5)

The back-propagation method is only applied to the coefficients of $z$, while
other hyper parameters of the trained GAN model remain unchanged in phase 2
which is Testing.

###  3.3 Motion Masks via Optical Flow

To identify fast moving objects in video frames, we use optical flow [15]
which creates a motion mask to capture moving objects. This technique
calculates motion between each pair of consecutive frames in the given input
video sequence $I$. Motion mask $M$ is estimated by using motion information
from consecutive frames of video sequences. Suppose $I_{t}$ and $I_{t-1}$ be
the two consecutive video frames in $I$ at time any instant $t$ and $t-1$,
respectively. Taking $U^{y}_{t,p}$ be the horizontal component and
$V^{x}_{t,p}$ be the vertical component of the motion vector at position $p$
computed between consecutive frames, the corresponding motion mask,
$M_{t}\in\\{0,1\\}$ will be estimated as:

$M_{t,p}=\begin{cases}1,~{}~{}if~{}~{}\sqrt{(U^{y}_{t,p})^{2}+(V^{x}_{t,p})^{2}}<T,\\\
0,~{}~{}otherwise.\end{cases}$ (6)

In the above equation, $T$ is a threshold of motion mask magnitude, and it is
computed by taking the average of all pixels in the estimated motion field.

##  4 Implementation

Our work is inspired by [29] and we have adapted the network [22] which is
based on DCGAN-tensorflow implementation [21]. Our proposed ForeGAN_RGBD model
has a generative network, G, takes a random noise vector drawn from a uniform
distribution and generates an image with dimension $64\times 64\times 3$. The
discriminator model, however, runs in reverse order of the generator. The
input images are fixed at $64\times 64\times 3$ running through a series of
five convolution layers with down sampling and the number of channels are
double the size of the previous convolution layer. Furthermore, the last layer
of the network is a two class softmax. For training the ForeGAN_RGBD model, we
use Adam [13] for optimization. The back-propagation steps in our model are
$2000$ during testing to generate the background image similar to the test
image. The hyper-parameters are discussed in detail in section 5.1 and 5.2.
Experiments are conducted using TitanX GPU with fixed input image samples of
size in training as well as testing. Since all the categories have a limited
number of video sequences, we have increased the training data by data
augmentation technique which is translation and rotation of all the training
images. Note that only those samples are considered in training whose ground
truth information is not available in the dataset. Similarly only those video
frames are considered for evaluation/testing purpose whose ground truth
information is available in the dataset.

##  5 Experiments

We have evaluated our proposed approach on $6$ different categories of SBM-
RGBD dataset [4] containing $27$ video sequences. The video sequences are very
challenging as they contain indoor and outdoor background scenes for
foreground detection in the presence of color as well as depth information.
The experimental process is divided into two phases, training and testing. The
detail explanation of both phases is presented as follows:

###  5.1 Phase 1: Training

The training is performed on $6$ categories of the SBM-RGBD dataset using the
proposed technique. We have trained two models independently one to learn of
color information (ForeGAN_RGB) of the video sequences and the other to learn
the depth information (ForeGAN_D) with $200$ epochs each. During the training
phase, the goal of our GAN network is to learn the color as well as depth
information of the background scene. However, the only difference between two
models is that the RGB model is trained on all kinds of video sequences
regardless of the presence of foreground objects while the depth model of GAN
is trained with the video sequences containing only depth information of
background scene. Both of the models are trained individually on scene-
specific information based on each category in SBM-RGBD dataset.

Figure 2: Visual results of each step of our propose ForeGAN_RGBD method. (a)
Input test RGB sequences $I_{t}$ as mentioned in section 3.3 (b) Input test
RGB sequences $I_{t-1}$ (c) Input test RGB sequence $I_{t-1}$ multiplied by
motion masks to identify moving objects (d) generated background samples by
our $ForeGAN\\_RGB$ model (e) Foreground objects detected by subtracting (b)
and (d) which is equation (3), (f) Ground truths (g) Test depth sequences (h)
Generated background depths by $ForeGAN\\_D$ model. The sequences represented
in $row~{}1$ is from ’Bootstrapping_ds’ category ’Bootstrapping’, $row~{}2$ is
from ’fall01cam1’ category ’Shadows’, $row~{}3$ is from ’Hallway ’ category
’Color Camouflage’ and $row~{}4$ is from ’Despatx_ds’ category ’Depth
Camouflage’.

###  5.2 Phase 2: Testing

The testing of the proposed models is performed independently on both trained
models, and the output of both models is fused to get the final results. All
the testing samples are also fixed as mentioned in section 4 and given as
input to the models individually for validation with back-propagation steps.
We have evaluated our proposed model by using $F$-measure score which is
calculated as follows:

$F=2\frac{P\times R}{P+R},$ (7) $P=\frac{T_{p}}{T_{p}+F_{p}},$ (8)
$R=\frac{T_{p}}{T_{p}+F_{n}},$ (9)

where $T_{p}$ is True positives, $T_{n}$ is True negatives, $F_{p}$ is False
positives, $F_{n}$ is False negatives, $F$ is $F$-Measure, $P$ is Precision,
and $R$ is Recall. For better foreground detection, the aim is to maximize $F$
score. Since the two models, $ForeGAN\\_RGB$ and $ForeGAN\\_D$ are trained in
way different way; they are also tested discretely, to improve the performance
of the foreground object detection. The trained model for color information
$ForeGAN\\_RGB$ is tested in such a way that the input test sequences are
first multiplied with motion mask to eliminate the moving foreground object
information (see Figure 2 (3)). This step helps our trained network to
generate the background image sample with similar information by back-
propagation steps. Later this generated background image is subtracted with
input test sample to extract the foreground objects (Figure 2 (a), (c) and
(d), Figure 1 (5)). The $ForeGAN\\_D$ is also evaluated in the same way, but
the only difference is that the testing sampling is given as it is to the
trained model which generates the similar depth information as test sample but
it will be a background scene. After that pixel by pixel subtraction of
generated background depth information and testing sample gives us foreground
object segmentation with a lot of noise in it mostly as shown in Figure 1 (8).
So motion masks are multiplied with foreground objects detected via depth
information to eliminate the noisy background depth pixels as shown in Figure
1 (8), (9) and (10) respectively. The last step is to fuse the information of
foregrounds detected by following fusion technique from both models to get the
final output

$ForeGAN\\_RGB+ForeGAN\\_D=ForeGAN\\_RGBD$

Table 1: Comparison of average $F$-measure score on SBM-RGBD dataset [4]
concerning depth and color features. The first and second best performing
methods are shown in red and blue colors respectively.

Methods | Depth features | RGB features |
---|---|---|---
BS-GFL | 0.21 | 0.63 |
DECOLOR | 0.20 | 0.56 |
cwisardH+ | — | — |
RGB-SOBS | — | 0.70 |
SRPCA | 0.24 | 0.71 |
ForeGAN_RGBD | 0.45 | 0.51 |

Table 2: Comparison of average $F$-measure score on SBM-RGBD dataset [4].
Figure 3 represents visual comparisons with other methods. The first and
second best performing methods are shown in red and blue colors respectively.

Categories | BS-GFL | DECOLOR | cwisardH+ | RGB-SOBS | SRPCA | ForeGAN_RGBD |
---|---|---|---|---|---|---|---
Illumination Changes | 0.2991 | 0.4861 | 0.4581 | 0.4527 | 0.4454 | 0.8158 |
Color Camouflage | 0.7333 | 0.7030 | 0.9510 | 0.4864 | 0.8329 | 0.9635 |
Depth Camouflage | 0.7540 | 0.7252 | 0.7648 | 0.8935 | 0.8083 | 0.9360 |
Out of Range | 0.7182 | 0.5874 | 0.8987 | 0.8527 | 0.8011 | 0.8726 |
Shadows | 0.5869 | 0.9051 | 0.9264 | 0.9218 | 0.7591 | 0.9271 |
Bootstrapping | 0.5711 | 0.7601 | 0.5669 | 0.8007 | 0.8098 | 0.8646 |
Average | 0.6104 | 0.6482 | 0.7609 | 0.7075 | 0.7698 | 0.8966 |

###  5.3 Foreground Object Segmentation by ForeGAN_RGBD

To highlight the significance of our proposed GAN model, we compared it with
$5$ state-of-the-art methods by $F$-measure as shown in Table 2. By using
original implementations of the authors, we have compared our proposed method
with BS-GFL [28], DECOLOR [31], cwisardH+ [6], RGB-SOBS [16] and SRPCA [9].
Qualitative comparisons of our proposed method are presented in Figure 3 and
$F$-measure score comparison with $5$ state-of-the-art methods are presented
in Tables 1 and 2. It can be seen in Table 1 that our proposed method has
performed well regarding foreground object segmentation by using depth
features as compared to all methods. However, by using only RGB features the
$F-measure$ is quite low which means our ForeGAN_RGB model does not detect the
foreground objects very well. The main reason behind this fact is that
sometimes there are too much background noisy pixel values in frame difference
detection of background-foreground images. Nevertheless, this problem is
eliminated in ForeGAN_D model because after frame difference of background
depth image generated by our network and test sample; it is masked with
optical flow motion information as shown in Figure 1. But upon combining the
features of depth as well as color information the $F-measure$ is increased
significantly in the case of our proposed algorithm as compared to all methods
as shown in Table 1.
The detail explanation of category wise comparison of our proposed method with
$5$ state-of-the-art methods is as follows:
Category: Illumination Changes contains 4 video sequences in SBM-RGBD dataset
[4]. The average $F$-measure scores among all four video sequences are shown
in Table 2. It can be seen in Table 2 that our proposed method ForeGAN_RGBD
has achieved maximum $F$-measure score among all the compared methods.
However, the second best performing method is cwisardH+ among all the compared
methods. Since this category represents challenging video sequences with
several illumination variations it can be seen in Table 2 that all compared
methods except our proposed method have a very low $F$-measure score. It is
because illumination variations pose a lot of challenges to all compared
methods, but our proposed method can generate the background sequence with
exact illumination condition. This aspect of our proposed method favors
precise foreground detection with challenging illumination conditions. Visual
results are represented in Figure 3.

Figure 3: Performance comparison of our proposed method with 5 state-of-the-
art methods. From left to right: each input sequence is selected from
different categories (a) ’fall01cam0’ category ’Bootstrapping’, (b)
’colorCam1’ category ’Color Camouflage’, (c) ’Wall’ category ’Depth
Camouflage’, (d) ’genSeq1’ category ’Illumination Changes’, (e) ’MultiPeople1’
category ’Out of Range’ and (f) ’shadows1’ category ’Shadows’. $row~{}1$: 6
images from input test video sequences $row~{}2$: Ground truth, $row~{}3$:
foreground detection by our proposed ForeGAN_RGBD method, $row~{}4$: SOBS-RGB,
$row~{}5$: DECOLOR, $row~{}6$: SRPCA and $row~{}7$: BS-GFL.

Category: Color Camouflage also contains 4 video sequences including
foreground objects that are very close in color to the background in SBM-RGBD
dataset [4]. The average $F$-measure scores among all four video sequences are
shown in Table 2 and it can be seen that only our proposed method and
cwisardH+ method has performed well. This category contains challenging video
sequences in which foreground objects are very close in color to the
background. Still, our proposed method achieved the highest $F$-measure score.
The main reason is that as our algorithm extracts foreground objects from RGB
and depth models independently, so if even one of the model suffers
performance degradation the other model can detect the foreground objects.
Qualitative results are shown in Figure 2 and 3 in comparison to ground truths
and state-of-the-art methods respectively.
Category: Depth Camouflage also contains 4 videos including foreground objects
very close in depth to the background in SBM-RGBD dataset [4]. It can be seen
Table 2 that our proposed method has outperformed all the compared methods and
RGB-SOBS method has achieved the second best score. The top performance of our
proposed method is due to the same reason mentioned previously in color
camouflage challenge discussion. The visual results are represented in Figure
2 (d) row 4 and Figure 3 (c). It can be seen in the qualitative analysis that
our proposed method has better results than all the compared methods except
RGB-SOBS method which has performed second best in this category.
Category: Out of Range has 5 videos sequences including foreground or
background objects that are too close to/far from the sensor [4]. It can be
seen in the Table 2 that our proposed method has achieved second best
$F$-measure score 0.8726 however cwisardH+ method achieved best score 0.8987.
The reason is that sometimes our proposed method cannot generate the perfect
background sequence for foreground objects which are too far from the sensor.
Visual results presented in Figure 3 (e) which shows the comparison of our
proposed method with all compared methods and ground truth information as
well.
Category: Shadows also contains 5 videos showing shadows caused by foreground
objects. These can be visible-light shadows in the RGB channels or IR shadows
in the depth channel [4]. In this challenging category, our proposed method
has also achieved the highest $F$-measure score 0.9271 but the cwisardH+
method also achieved the best score with a minimal difference as shown in
Table 2. The visual results shown in Figure 2 (d) row 2 and Figure 3 (f)
represents that our proposed method can detect foreground objects even in the
presence of shadows.
Category: Bootstrapping also contains 5 video sequences including foreground
objects in almost all their frames [4]. It can be seen in Table 2 that our
proposed method ForeGAN_RGBD has also achieved highest $F$-measure score and
SRPCA has achieved second highest $F$-measure score. It is a very challenging
category for all compared methods because it contains video sequences
including foreground objects in almost all their frames but still our proposed
method outperformed all compared methods. The main reason behind this fact is
that our proposed method can learn the background scene information even if
the foreground objects have left the scene once. This aspect of our proposed
GAN network improves its efficiency in foreground object detection with
challenging scenarios like bootstrapping. For instance visual results are
shown in Figure 2 (c) and (d) row 1 and Figure 3 (a).

##  6 Conclusion

In this study, we present the foreground detection algorithm based on
generative adversarial network (GAN) exploiting color as well as depth
information. Our goal is to do foreground segmentation in the presence of
various major challenges in background scenes in real environments. To handle
these challenges we present a fusion based foreground detection algorithm
which exploits color as well as depth information, with the help of Generative
Adversarial Network (GAN). Our proposed GAN model is trained in an
unsupervised fashion on color and depth information independently with various
challenges. After that for testing, the GAN model has to generate the same
background samples as test samples with similar statistics via back-
propagation technique. The generated background samples are then subtracted
from the given test samples to detect foreground objects from color as well as
depth information, and final results are calculated by the fusion of both
model’s segmented outputs. The comparison of our proposed method with five
state-of-the-art methods highlights the strength of our algorithm for
foreground object segmentation in the presence of various challenging
conditions.

##  7 acknowledgements

This research was supported by Development project of leading technology for
future vehicle of the business of Daegu metropolitan city(No. 20181030).

## References

  * [1] Bouwmans, T., Javed, S., Zhang, H., Lin, Z., Otazo, R.: On the applications of robust pca in image and video processing. Proceedings of the IEEE 106(8), 1427–1457 (2018)
  * [2] Bouwmans, T., Maddalena, L., Petrosino, A.: Scene background initialization: a taxonomy. Pattern Recognition Letters 96, 3–11 (2017)
  * [3] Bouwmans, T., Zahzah, E.H.: Robust pca via principal component pursuit: A review for a comparative evaluation in video surveillance. Computer Vision and Image Understanding 122, 22–34 (2014)
  * [4] Camplani, M., Maddalena, L., Gabriel, M., Petrosino, A., Salgado, L.: Rgb-d dataset: Background learning for detection and tracking from rgbd videos. In: IEEE ICIAP-Workshops (2017)
  * [5] Chen, Y., Zou, W., Tang, Y., Li, X., Xu, C., Komodakis, N.: Scom: Spatiotemporal constrained optimization for salient object detection. IEEE Transactions on Image Processing 27(7), 3345–3357 (2018)
  * [6] De Gregorio, M., Giordano, M.: Cwisardh+: Background detection in rgbd videos by learning of weightless neural networks. In: International Conference on Image Analysis and Processing. pp. 242–253. Springer (2017)
  * [7] Dosovitskiy, A., Brox, T.: Inverting visual representations with convolutional networks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 4829–4837 (2016)
  * [8] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. In: Advances in neural information processing systems. pp. 2672–2680 (2014)
  * [9] Javed, S., Bouwmans, T., Sultana, M., Jung, S.K.: Moving object detection on rgb-d videos using graph regularized spatiotemporal rpca. In: International Conference on Image Analysis and Processing. pp. 230–241. Springer (2017)
  * [10] Javed, S., Mahmood, A., Al-Maadeed, S., Bouwmans, T., Jung, S.K.: Moving object detection in complex scene using spatiotemporal structured-sparse rpca. IEEE Transactions on Image Processing (2018)
  * [11] Javed, S., Mahmood, A., Bouwmans, T., Jung, S.K.: Spatiotemporal low-rank modeling for complex scene background initialization. IEEE Transactions on Circuits and Systems for Video Technology (2016)
  * [12] Javed, S., Mahmood, A., Bouwmans, T., Jung, S.K.: Background–foreground modeling based on spatiotemporal sparse subspace clustering. IEEE Transactions on Image Processing 26(12), 5840–5854 (2017)
  * [13] Kinga, D., Adam, J.B.: A method for stochastic optimization. In: International Conference on Learning Representations (ICLR). vol. 5 (2015)
  * [14] Linden, A., Kindermann, J.: Inversion of multilayer nets. In: Proc. Int. Joint Conf. Neural Networks. vol. 2, pp. 425–430 (1989)
  * [15] Liu, C., et al.: Beyond pixels: exploring new representations and applications for motion analysis. Ph.D. thesis, Massachusetts Institute of Technology (2009)
  * [16] Maddalena, L., Petrosino, A.: The sobs algorithm: what are the limits? In: Computer Vision and Pattern Recognition Workshops (CVPRW), 2012 IEEE Computer Society Conference on. pp. 21–26. IEEE (2012)
  * [17] Maddalena, L., Petrosino, A.: Background subtraction for moving object detection in rgbd data: A survey. Journal of Imaging 4(5), 71 (2018)
  * [18] Maddalena, L., Petrosino, A., et al.: A self-organizing approach to background subtraction for visual surveillance applications. IEEE Transactions on Image Processing 17(7), 1168 (2008)
  * [19] Mahendran, A., Vedaldi, A.: Understanding deep image representations by inverting them. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 5188–5196 (2015)
  * [20] Mordvintsev, A., Olah, C., Tyka, M.: Inceptionism: Going deeper into neural networks. Google Research Blog. Retrieved June 20(14), 5 (2015)
  * [21] Radford, A., Metz, L., Chintala, S.: Unsupervised representation learning with deep convolutional generative adversarial networks. arXiv preprint arXiv:1511.06434 (2015)
  * [22] Schlegl, T., Seeböck, P., Waldstein, S.M., Schmidt-Erfurth, U., Langs, G.: Unsupervised anomaly detection with generative adversarial networks to guide marker discovery. In: International Conference on Information Processing in Medical Imaging. pp. 146–157. Springer (2017)
  * [23] Shimada, A., Nagahara, H., Taniguchi, R.i.: Background modeling based on bidirectional analysis. In: Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference on. pp. 1979–1986. IEEE (2013)
  * [24] Stauffer, C., Grimson, W.E.L.: Adaptive background mixture models for real-time tracking. In: Computer Vision and Pattern Recognition, 1999. IEEE Computer Society Conference on. vol. 2, pp. 246–252. IEEE (1999)
  * [25] Sultana, M., Mahmood, A., Javed, S., Jung, S.K.: Unsupervised deep context prediction for background foreground separation. arXiv preprint arXiv:1805.07903 (2018)
  * [26] Vaswani, N., Bouwmans, T., Javed, S., Narayanamurthy, P.: Robust pca and robust subspace tracking. arXiv preprint arXiv:1711.09492 (2017)
  * [27] Wu, Y., He, X., Nguyen, T.Q.: Moving object detection with a freely moving camera via background motion subtraction. IEEE Transactions on Circuits and Systems for Video Technology 27(2), 236–248 (2017)
  * [28] Xin, B., Tian, Y., Wang, Y., Gao, W.: Background subtraction via generalized fused lasso foreground modeling. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 4676–4684 (2015)
  * [29] Yeh, R.A., Chen, C., Lim, T.Y., Schwing, A.G., Hasegawa-Johnson, M., Do, M.N.: Semantic image inpainting with deep generative models. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 5485–5493 (2017)
  * [30] Zhang, T., Liu, S., Ahuja, N., Yang, M.H., Ghanem, B.: Robust visual tracking via consistent low-rank sparse learning. International Journal of Computer Vision 111(2), 171–190 (2015)
  * [31] Zhou, X., Yang, C., Yu, W.: Moving object detection by detecting contiguous outliers in the low-rank representation. IEEE T-PAMI 35(3), 597–610 (2013)