# Gravitational Harmonium: Gravitationally Induced Entanglement in a Harmonic
Trap

Jackson Yant<EMAIL_ADDRESS>Department of Physics and
Astronomy, Dartmouth College, Hanover, New Hampshire 03755, USA  Miles
Blencowe<EMAIL_ADDRESS>Department of Physics and Astronomy,
Dartmouth College, Hanover, New Hampshire 03755, USA

###### Abstract

Recent work has shown that it may be possible to detect gravitationally
induced entanglement in tabletop experiments in the not-too-distant future.
However, there are at present no thoroughly developed models for this type of
experiment where the entangled particles are treated more fundamentally as
excitations of a relativistic quantum field, and with the measurements modeled
using expectation values of field observables. Here we propose a thought
experiment where two particles (i.e., massive scalar field quanta) are
initially prepared in a superposition of coherent states within a common
three-dimensional (3D) harmonic trap. The particles then develop entanglement
through their mutual gravitational interaction, which can be probed through
particle position detection probabilities. The present work gives a non-
relativistic quantum mechanical analysis of the gravitationally induced
entanglement of this system, which we term the ‘gravitational harmonium’ due
to its similarity to the harmonium model of approximate electron interactions
in a helium atom; the entanglement is operationally determined through the
matter wave interference visibility. The present work serves as the basis for
a subsequent investigation, which models this system using quantum field
theory, providing further insights into the quantum nature of gravitationally
induced entanglement through relativistic corrections, together with an
operational procedure to quantify the entanglement.

##  I  Introduction

Quantum mechanics and general relativity are our two most fundamental,
experimentally tested, and predictive theories. However, these theories are
mutually incompatible and no comparably successful theory unifying the two has
been formulated, while little experimental progress has been made. Recently,
Bose et al. [1] and Marletto and Vedral [2] (BMV) proposed in-principle-
realizable experiments in which two particles could become entangled by their
mutual gravitational interaction [1, 2]. Detection of gravitationally induced
entanglement would then provide novel empirical evidence that might guide the
development of new theories aiming at resolving the tensions between classical
gravity and quantum mechanics.

The BMV proposals have inspired a great deal of interest in the quantum
information science community, with ongoing discussions about the theoretical
implications of detecting gravitationally induced entanglement. Some contend
that a positive detection would provide definitive evidence for the quantum
nature of the gravitational field [1, 2, 3, 4, 5, 6, 4, 7, 8, 9]. Others have
pointed out possible ways in which a classical model of gravity could induce
entanglement in such an experiment [10, 11, 12, 13, 14, 15, 16]. Some have
raised points about how several of the principles of quantum information
theory used to assert that gravitationally induced entanglement necessitates
quantum gravity, are not clearly satisfied in a quantum field theoretic
setting [17, 18, 19, 20]. In response, there have also been several
investigations into what implications can be drawn from gravitationally
induced entanglement without necessarily relying on such assumptions [21, 22,
23, 24].

Several variations of a gravitationally induced entanglement experiment have
been proposed [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. To
date, none of these proposed experiments have been modelled using a fully
quantum field theoretic description of both gravity, the particles involved
and the measures of entanglement generated between the latter. Such an ab
initio approach would allow one to explore the quantum nature of the
gravitational field, including relativistic effects [39]. In this work, we
seek to develop a gravitationally induced entanglement thought experiment that
adapts the BMV proposal to enable such a description.

Figure 1: Scheme of the thought experiment for generating and probing
gravitationally-induced entanglement between two massive particles. Snapshots
(a)-(e) show the coherent state position maxima of the two particles, each
initially in a superposition state, at subsequent times. When the particle
wavefunction components in a given superposition state overlap, interference
fringes develop in the spatial dependence of the particle detection
probability along the oscillation axis of the particle. The fringe visibility
decreases over time, corresponding to increasing entanglement between the two
particles.

Our thought experiment considers two particles that are located within a 3D
harmonic trap. The initial condition is a product state wherein each particle
is prepared in a displaced coherent state superposition along mutually
orthogonal axes. As the coherent state wavefunction components of one of the
particles evolve, they will overlap periodically, producing interference
fringes in the spatial probability distribution for detecting the particle.
The gravitationally-induced entanglement that develops between the two
particles will cause the visibility of these fringes to decrease. Figure 1
illustrates this scheme. We term this model the ‘gravitational harmonium’
after its similarity to a model for approximating the Helium electrons’
Coulombic trapping potentials with harmonic potentials [40].111While we view
such a scheme as a thought experiment that allows for a relativistic quantum
field theoretic analysis, it does resemble recent optical trapping schemes for
massive particles that share a similar goal: to probe gravitational
entanglement in quantum systems [41].

The present work analyzes the gravitational harmonium using non-relativistic
quantum mechanics, serving as the basis for a future investigation that will
apply a quantum field theory analysis to the same system [42]. In particular,
the gravitational harmonium is readily suited for a full, quantum field
theoretic treatment. Under low energy laboratory conditions, we can model such
a thought experiment by a massive real scalar field $\phi$ coupled to gravity
through the Einstein Hilbert action, expanded to second order in
$\kappa=\sqrt{32\pi G}$ (with units $\hbar=c=1$) [43, 44]:

$S[\phi,h_{\mu\nu}]=S_{M}[\phi]+S_{E}[h_{\mu\nu}]+S_{I}[\phi,h_{\mu\nu}].$ (1)

The action for the scalar field, $\phi$, contains a harmonic potential that is
spatially centered, e.g., at rectilinear coordinate location ${\bf{r}}=0$ and
given by

$S_{M}[\phi]=-\frac{1}{2}\int
d^{4}x\,\left[\eta^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+m^{2}\left(1+\omega^{2}r^{2}\right)\phi^{2}\right],$
(2)

while the gravity and interaction actions are respectively

$\begin{split}S_{E}[h_{\mu\nu}]=&\int
d^{4}x\left(-\frac{1}{2}\partial^{\rho}h^{\mu\nu}\partial_{\rho}h_{\mu\nu}+\partial_{\nu}h^{\mu\nu}\partial^{\rho}h_{\mu\rho}-\partial_{\mu}h\partial_{\nu}h^{\mu\nu}+\frac{1}{2}\partial^{\mu}h\partial_{\mu}h\right),\\\
S_{{I}}=&\int
d^{4}x\left(\frac{\kappa}{2}T^{\mu\nu}\left(\phi\right)h_{\mu\nu}+\frac{\kappa^{2}}{4}U^{\mu\nu\rho\sigma}\left(\phi\right)h_{\mu\nu}h_{\rho\sigma}\right),\end{split}$
(3)

where $T_{\mu\nu}\left(\phi\right)$ is the scalar field energy-momentum
tensor, and $U_{\mu\nu\rho\sigma}\left(\phi\right)$ is a quadratic tensor in
$\phi$ [44]. Decoherence aspects of certain states within such a framework
have been considered in Ref. [45], and in a 0D model [46]. We may then
describe the fringe visibility by using the expectation value of the particle
detection probability given in the Schrödinger picture,

${\mathrm{Tr}}\left[\rho\left(t\right)\left(V^{-1}\int_{V}d{\bf{x}}\phi\left({\bf{x}},0\right)\right)^{2}\right]=V^{-2}\int_{V}d{\bf{x}}d{}{\bf
x}^{\prime}{\mathrm{Tr}}\left[\rho\left(t\right)\phi\left({\bf{x}},0\right)\phi\left({\bf{x}}^{\prime},0\right)\right],$
(4)

where $\rho\left(t\right)$ is the density operator for the evolving two-
particle state and $V$ is some small coordinate averaging volume corresponding
to the spatial region occupied by the particle detector.

In Section II we solve the corresponding non-relativistic Schrödinger equation
to derive the approximate time evolution of the two-particle, initial
superposition pair state of this system as gravitational entanglement between
the two particles develops, with an accompanying fringe visibility reduction
in the spatial dependence of the single particle detection probability;
further details are given in three appendices. As discussed above, this
analysis will be used as a reference for future work where the corresponding
quantum dynamics will be obtained within the above-described quantum field
theoretic description [42]; we hope to gain new insights into how gravity
induces quantum entanglement by comparing the approximate non-relativistic
model developed here with the more fundamental one involving relativistic
quantum fields. Section III considers the effect that alternative, modified
gravity interaction potentials would have on the non-relativistic quantum
mechanical model for entanglement generation and the corresponding reduction
in visibility. Section IV gives some concluding remarks.

##  II NRQM Gravitational Harmonium

Consider two particles with identical mass $m$ that are confined to a harmonic
trap with frequency $\omega$ and interacting through gravity. The non-
relativistic quantum dynamics of this system is described by the Schrödinger
equation with Hamiltonian given by

$H=\frac{p_{1}^{2}}{2m}+\frac{p_{2}^{2}}{2m}+\frac{1}{2}m\omega^{2}x_{1}^{2}+\frac{1}{2}m\omega^{2}x_{2}^{2}-\frac{Gm^{2}}{|\bm{x}_{1}-\bm{x}_{2}|},$
(5)

where $\bm{x}_{1}$ and $\bm{x}_{2}$ are the position coordinates of particles
1 and 2, respectively. This Hamiltonian is qualitatively similar to that of
the Helium atom [with the substitution $Gm^{2}\rightarrow-e^{2}$ in Eq. (5)],
where the electrons’ Coulomb confining potentials are approximately replaced
by harmonic confining potentials, and is known as ‘Hooke’s atom’ or the
‘Harmonium’ [40]. Because of this similarity, we name this model the
‘gravitational harmonium.’ Hamiltonian (5) can be conveniently expressed in
dimensionless form with the position and time coordinate subsitutions
$x_{i}\rightarrow\sqrt{\frac{m\omega}{\hbar}}x_{i}$, $t\rightarrow\omega t$,
and $H\rightarrow H/\hbar\omega$, giving:

$H=\frac{1}{2}\left(p_{1}^{2}+p_{2}^{2}+x_{1}^{2}+x_{2}^{2}\right)-\frac{\sqrt{2}g_{0}}{|\bm{x}_{1}-\bm{x}_{2}|},$
(6)

where $g_{0}=G\sqrt{\frac{m^{5}}{2\hbar^{3}\omega}}$ is a dimensionless
parameter characterizing the Newtonian gravitational interaction between the
two particles. In these scaled, dimensionless coordinates, the Schrödinger
equation is given as $i\partial_{t}\psi=H\psi$.

The initial state of our system is taken to be a product state wherein
particle 1 is in a superposition of coherent states $|\bm{\alpha}_{1}\rangle$
displaced along the $x$ axis direction [i.e.,
$\bm{\alpha}_{1}=(\alpha_{1x},0,0)$, with $\alpha_{1x}\neq 0$], while particle
2 is in a superposition of coherent states $|\bm{\alpha}_{2}\rangle$ displaced
along the orthogonal $y$ axis direction [i.e.,
$\bm{\alpha}_{2}=(0,\alpha_{2y},0)$, with $\alpha_{2y}\neq 0$] (Fig. 1):

$\begin{split}\ket{\psi}=&N\left(\ket{\bm{\alpha}_{1a}}+\ket{\bm{\alpha}_{1b}}\right)\left(\ket{\bm{\alpha}_{2a}}+\ket{\bm{\alpha}_{2b}}\right)\\\
=&N\left(\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2a}}+\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2b}}+\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2a}}+\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2b}}\right),\end{split}$
(7)

where $N$ is a normalization factor. Here, the single particle coherent states
are each parameterized by a complex vector
$\bm{\alpha}=\left(\expectationvalue{\bm{x}}+i\expectationvalue{\bm{p}}\right)/\sqrt{2}$.
For the initial state (7), both particles subsequently remain undisplaced in
the $z$-coordinate axis direction and we therefore suppress further mention of
the $z$ coordinate and treat $\bm{x}_{1}$ and $\bm{x}_{2}$ as two dimensional
vectors from now on; the coherent state parameter vectors $\bm{\alpha}$ are
then given by a pair of $x$ and $y$ component complex numbers for each
particle [$\bm{\alpha}_{1}=(\alpha_{1x},0)$ and
$\bm{\alpha}_{2}=(0,\alpha_{2y})$].

In order to determine the subsequent Schrödinger evolution of the initial
state in Eq. (7) with Hamiltonian (6), we first consider the evolution of one
of the individual coherent state products
$\ket{\bm{\alpha}_{1},\,\bm{\alpha}_{2}}$ for particles 1 and 2 that make up
the superposition on the second line of Eq. (7). We can then determine the
full state evolution by linearity. It is convenient to work in terms of
symmetric center of mass (SCOM) coordinates given by

$\bm{R}=\frac{1}{\sqrt{2}}\left(\bm{x}_{1}+\bm{x}_{2}\right)\ \mathrm{and,}\
\bm{r}=\frac{1}{\sqrt{2}}\left(\bm{x}_{1}-\bm{x}_{2}\right),$ (8)

with $\bm{R}$ referred to as the center of mass coordinate and $\bm{r}$ the
relative difference coordinate. In terms of these coordinates, Hamiltonian (6)
decomposes into the sum of the center of mass and relative difference
coordinate Hamiltonians:

$H=\frac{1}{2}\left(p_{R}^{2}+R^{2}\right)+\frac{1}{2}\left(p_{r}^{2}+r^{2}-\frac{g_{0}}{r}\right).$
(9)

The ket $\ket{\bm{\alpha}_{1},\,\bm{\alpha}_{2}}$ can equivalently be
described by coherent state parameter vectors for the SCOM coordinates as
follows:

$\bm{\alpha}_{R}=\frac{1}{\sqrt{2}}\left(\bm{\alpha}_{1}+\bm{\alpha}_{2}\right)\
\mathrm{and,}\
\bm{\alpha}_{r}=\frac{1}{\sqrt{2}}\left(\bm{\alpha}_{1}-\bm{\alpha}_{2}\right).$
(10)

These parameter vectors are nonzero in both their $x$ and $y$ components. From
Eq. (9), the Hamiltonian is separable in the SCOM coordinates, and the center
of mass coordinate evolves as a simple, 2D quantum harmonic oscillator.

In order to approximate the evolution of the relative difference coordinate,
we Taylor expand the potential terms in our Hamiltonian to second order about
the position expectation values. For coherent states, whose wavefunctions are
given by Gaussians of the form

$\psi\left(\bm{r}\right)=\exp
i\left[\left(\bm{r}-\expectationvalue{\bm{r}}\right)^{T}\Omega\left(\bm{r}-\expectationvalue{\bm{r}}\right)+\expectationvalue{\bm{p}_{r}}\left(\bm{r}-\expectationvalue{\bm{r}}\right)+\gamma\right],$
(11)

the Gaussian parameters then evolve as follows [47]:

$\displaystyle\expectationvalue{\dot{\bm{r}}}$ $\displaystyle=\frac{\partial
H}{\partial\expectationvalue{\bm{p}_{r}}},\quad\expectationvalue{\dot{\bm{p}_{r}}}=-\frac{\partial
H}{\partial\expectationvalue{\bm{r}}},$ (12) $\displaystyle\dot{\Omega}$
$\displaystyle=-2\Omega^{2}-\frac{1}{2}\frac{\partial^{2}V}{\partial
r_{i}\partial r_{j}}\Big{|}_{\bm{r}=\expectationvalue{\bm{r}}},$ (13)
$\displaystyle\dot{\gamma}$
$\displaystyle=i\,\mathrm{Tr}\,\Omega+\mathcal{L}.$ (14)

Here, $\mathcal{L}$ is the system’s classical Lagrangian evaluated for the
solutions $\expectationvalue{\bm{r}}$ and $\expectationvalue{\bm{p}_{r}}$, and
$V$ includes the harmonic and gravitational potential terms. We consider
initial coherent state parameter values $\bm{\alpha}_{1}=(\alpha,0)$ and
$\bm{\alpha}_{2}=(0,i\alpha)$ for some real number $\alpha$. This corresponds
to the condition that the average relative distance
$\expectationvalue{\bm{r}}$ between particles 1 and 2 remains fixed in the
absence of their mutual gravitational attraction ($g_{0}=0$). Approximate
solutions to Eqs. (12)-(14) are derived in Appendix A, and their numerical
validation is given in Appendix C.

As an aside, it is interesting to note that the initial product state
$\ket{\bm{\alpha}_{1},\bm{\alpha}_{2}}$ evolves into a state where particles 1
and 2 are entangled due to nonzero off-diagonal terms in $\Omega$. The
entanglement entropy is approximately

$\displaystyle
S\approx\frac{9g_{0}^{2}t^{2}}{64\alpha^{6}}\left[1-\log\left(\frac{9g_{0}^{2}t^{2}}{128\alpha^{6}}\right)\right],$
(15)

for $\frac{g_{0}t}{\alpha^{3}}\ll 1$. This is in contrast to other
gravitational entanglement proposals, where the approximate time-evolved state
only develops entanglement due to the different gravitationally-induced
relative phases between various localized position state components resulting
from their initial spatial superpositions. Such entanglement results in our
method of approximation since it keeps track of the motion of the particles’
expectation values subject to the gravitational interaction between them.
However, the gravitationally-induced relative phase terms will provide the
dominant contribution to the entanglement between the two particles when we
consider superpositions of spatially localized states as we show in the
following.

In order to isolate the entanglement arising from the gravitationally-induced
relative phases, we only keep the time-dependent phase contributions to each
initial product state component of the full superposition state that arise
from evaluating the action along the gravitationally perturbed two-particle
trajectories:

$\displaystyle\ket{\psi}=$ $\displaystyle N\left(e^{i\delta
S_{aa}}\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2a}}+e^{i\delta
S_{ab}}\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2b}}+e^{i\delta
S_{ba}}\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2a}}+e^{i\delta
S_{bb}}\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2b}}\right),$ (16)

where $N$ is the normalization. The $\delta S$ terms are given by the
differences between the action along the classical path for particles subject
to both the harmonic and gravitational interactions and the action for the
harmonic interaction alone, where the initial conditions are given by the
position and momentum expectation values for the coherent product state
component the $e^{i\delta S}$ phase terms are each multiplying.

We choose the following initial coherent state parameter values related to the
individual particle (not SCOM) coordinates:
$\alpha_{1ax}=\alpha,\alpha_{1bx}=\alpha e^{-i\pi/4},\alpha_{2ay}=\alpha
e^{i\pi/2},$ and $\alpha_{2by}=\alpha e^{i\pi/4}$, with $\alpha$ a real number
[see Fig. 1(a)]. This corresponds to components
$\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2a}}$ and
$\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2b}}$ both being in the configuration
described for the single coherent product state component above where the
particles’ relative difference coordinate is constant in our approximation.
These two pairs are out of phase with each other by a phase factor of $\pi/4$.
This implies that the components $\ket{\bm{\alpha}_{1a},\,\bm{\alpha}_{2b}}$
and $\ket{\bm{\alpha}_{1b},\,\bm{\alpha}_{2a}}$ have coherent state parameter
values for the $r_{x}$ and $r_{y}$ components that are out of phase by
$\pm\pi/4$, which does not correspond exactly to the evolution considered for
the single component case above, but may be approximated by the same methods
as described in Appendix A. The entanglement entropy for the state (16) is
approximated in Appendix B and we obtain

$\displaystyle S\approx$
$\displaystyle\left(\frac{c\,g_{0}}{\alpha}t\right)^{2}\left\\{1-\log\left[\left(\frac{c\,g_{0}}{\alpha}t\right)^{2}\right]\right\\},$
(17)

where $c$ is an $\mathcal{O}(10^{-1})$ parameter defined in Appendix B.
Comparing the entropy (15) arising for the individual initial two-particle
product state components (7) with the relative phase entanglement entropy
(17), we see that the former scales as $\alpha^{-6}$ while the latter scales
as $\alpha^{-2}$. Recall that $\alpha$ is chosen as a real number that can be
interpreted as the initial maximum displacement of the coherent state position
space wave function in units of the harmonic oscillator zero point
uncertainty. In any conceivable, realistic implementation, we expect
$\alpha\ggg 1$, so that the entanglement entropy for the relative phase terms
dominates over the entropy arising for the single component terms. This
justifies the approximation of keeping only the relative phase terms in the
time-evolved state for the purposes of determining the gravitationally-induced
entanglement between the two harmonically trapped particles.

While the growth of entanglement entropy as given by Eq. (17) demonstrates
that the evolving state [resulting from the initial state (7)] develops
genuine entanglement between the two particles, the entanglement entropy is
not itself an observable; we still require a means for measuring the
entanglement. Consider the probability distribution for detecting particle 1
on the $x$ axis. In terms of the approximate evolved state (16), we have

$\begin{split}\left|\bra{x_{1}}\ket{\psi}\right|^{2}=&|N|^{2}\,\left\\{2\left[1+\mathrm{Re}\left(e^{i\left(\delta
S_{aa}-\delta
S_{ab}\right)}\bra{\bm{\alpha}_{2b}}\ket{\bm{\alpha}_{2a}}\right)\right]|\psi_{1a}\left(x_{1}\right)|^{2}\right.\\\
&\left.+2\left[1+\mathrm{Re}\left(e^{i\left(\delta S_{ba}-\delta
S_{bb}\right)}\bra{\bm{\alpha}_{2b}}\ket{\bm{\alpha}_{2a}}\right)\right]|\psi_{1b}\left(x_{1}\right)|^{2}\right.\\\
&\left.+2\mathrm{Re}\left[\psi_{1a}\left(x_{1}\right)\psi^{*}_{1b}\left(x_{1}\right)\left(e^{i\left(\delta
S_{aa}-\delta S_{ba}\right)}+e^{i\left(\delta S_{ab}-\delta
S_{bb}\right)}\right.\right.\right.\\\ &\left.\left.\left.+e^{i\left(\delta
S_{aa}-\delta
S_{bb}\right)}\bra{\bm{\alpha}_{2b}}\ket{\bm{\alpha}_{2a}}+e^{i\left(\delta
S_{ab}-\delta
S_{ba}\right)}\bra{\bm{\alpha}_{2a}}\ket{\bm{\alpha}_{2b}}\right)\right]\right\\},\end{split}$
(18)

where $\psi_{1a}\left(x_{1}\right)=\bra{x_{1}}\ket{\bm{\alpha}_{1a}}$ etc.
When the average positions for $\psi_{1a}$ and $\psi_{1b}$ overlap, which
corresponds to their classical paths crossing in configuration space, we have
$|\psi_{1a}|=|\psi_{1b}|$. The third term in the above expression will then
have spatially oscillatory terms:

$\begin{split}\psi_{1a}\psi_{1b}^{*}=&|\psi_{1a}|^{2}\exp\left[i\sqrt{2}\left(\mathrm{Im}\,\alpha_{1ax}-\mathrm{Im}\,\alpha_{1bx}\right)\left(x_{1}-\frac{1}{\sqrt{2}}\mathrm{Re}\,\alpha_{1ax}\right)\right]=|\psi_{1a}|^{2}e^{i\Phi(x_{1})},\end{split}$
(19)

where here $\Phi$ denotes the phase term in the middle expression. These
oscillatory terms appear as fringes modulated by a Gaussian envelope, as
indicated schematically in Fig. 1. Expressing the probability distribution
(18) in the form
$|\bra{x_{1}}\ket{\psi}|^{2}=\chi+2[\mathrm{Re}\,\mu\cos\Phi(x_{1})-\mathrm{Im}\,\mu\sin\Phi(x_{1})]$,
the extrema of the distribution occur at $\chi\pm 2|\mu|$. The fringe
visibility $\mathcal{V}$, defined as the ratio of the difference between the
extrema to their sum, is given approximately as

$\mathcal{V}=\frac{2|\mu|}{\eta}\approx
1-\frac{1}{8}\left(\frac{cg_{0}t}{\alpha}\right)^{2}$ (20)

for $\alpha\gg 1$; the constant $c$ is the same as that appearing in Eq. (17).
Comparing Eq. (17) with Eq. (20), we can see that the visibility scales
monotonically with the entanglement entropy (Fig. 2). As the entanglement
entropy increases, the reduced state for either particle behaves increasingly
as an effective classical mixture of coherent states and so the quantum
interference in the particle position detection probability distribution
decreases correspondingly. Therefore, the fringe visibility can be employed as
an operational measure of the entanglement entropy for probing gravitationally
induced entanglement in our harmonic particle trap scheme. Using, for example,
masses $m\sim 10^{-14}\,\mathrm{kg}$ and a separation distance $\sim
10^{2}\,\mu\mathrm{m}$ as in Ref. [1], and an optical trapping potential
frequency $\omega\sim 10^{2}\,\mathrm{kHz}$ as in Ref. [41], gives a value
$\frac{g_{0}}{\alpha}\sim 10^{-6}$ which would result in percentage level
changes in the fringe visibility in times of order one second.

Recently, a rather general quantum interferometric identity was obtained that
relates the concurrence ($\mathcal{C}$), visibility ($\mathcal{V}$), and
which-path information (“particleness” $\mathcal{P}$) as complementary
measures for certifying the quantum aspects of gravitational entanglement
interactions [48]:

$\mathcal{V}^{2}+\mathcal{P}^{2}+\mathcal{C}^{2}=1.$ (21)

For the entangled state considered in this work, and in the limit of large
$\alpha$ where we can treat the coherent state component pairs as mutually
orthogonal to a good approximation, our state can be alternatively interpreted
as that arising for the two interferometer setup described in Ref. [48], and
using the forms for $\mathcal{P,V,}$ and $\mathcal{C}$ given there. For our
system state, we have: $\mathcal{P}\approx|(|e^{i\delta
S_{aa}}|^{2}+|e^{i\delta S_{ab}}|^{2})-(|e^{i\delta S_{ba}}|^{2}+|e^{i\delta
S_{bb}}|^{2})|=0$, $\mathcal{V}^{2}\approx 2|e^{i\delta S_{aa}-i\delta
S_{ba}}+e^{i\delta S_{aa}-i\delta S_{bb}}|\approx
1-\frac{1}{4}\left(\frac{cg_{0}}{\alpha}t\right)^{2}$ and
$\mathcal{C}^{2}\approx 2|e^{i\delta S_{aa}+i\delta S_{bb}}-e^{i\delta
S_{ab}+i\delta
S_{ba}}|\approx\frac{1}{4}\left(\frac{cg_{0}}{\alpha}t\right)^{2}$, hence
obeying the identity (21) to the leading order in $\alpha^{-1}$ with which we
calculate these measures.

Figure 2: Entanglement entropy $S$ (units of nits, where 1 nit $=\frac{1}{\ln
2}$ bits) for the reduced state of one of the particles and the fringe
visibility $1-{\mathcal{V}}$ for particle 1 as a function of the scaled time
parameter $\frac{g_{0}}{\alpha}t$.

##  III Alternative Potentials

The method of approximating the evolution of the state considered in this work
is not unique to the Newtonian gravitational potential; it can also be applied
to other weak interactions that are treated as perturbations to the harmonic
trap. Of potential interest is a Coulomb type potential generalized to
arbitrary dimensions or a Yukawa potential [49]. Coulomb potentials of various
dimensions may be useful if practical applications of this proposal use
extended objects to effectively reduce the dimension of the interaction.
Yukawa potentials arise in theories with massive gravitons or in theories
containing light scalar particles such as dilatons, and both could be used as
alternatives to test for possible deviations from the inverse square law [50].
Here, we calculate the effect these different potentials would have on the
entanglement entropy or fringe visibility.

Changing the potential results in different values for the first order
corrections to the classical trajectories used in the calculation of the
entanglement entropy given by Eq. (17), which in turn results in different
forms for the $\delta S$ terms that appear in the phase factors of the time
evolved state. Since the entanglement entropy and fringe visibility can be
expressed as functions of these $\delta S$ terms, it suffices to compute the
modified forms of the latter. A new interaction potential equates to changing
the Kameltonian $K$ in Eq. (34), which will result in a modification to the
values

$\displaystyle\dot{\bar{\eta}}_{1x}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}-\frac{\partial
K}{\partial\phi_{x}}\Big{|}_{0}dt=\frac{1}{2\pi}\int_{0}^{2\pi}-\frac{\partial
K}{\partial
r}\frac{\eta_{x0}\sin\left(t+\phi_{x0}\right)\cos\left(t+\phi_{x0}\right)}{r}\Big{|}_{0}dt$
(22) $\displaystyle\dot{\bar{\phi}}_{1x}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial
K}{\partial\eta_{x}}\Big{|}_{0}dt=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial
K}{\partial r}\frac{\sin^{2}\left(t+\phi_{x0}\right)}{r}\Big{|}_{0}dt,$ (23)

and similarly for the $y$ components. Here, $\eta$ and $\phi$ are the Hamilton
Jacobi coordinates for the classical harmonic oscillator which are used in our
approximation in Appendix A. For the initial conditions considered above,
i.e., $\alpha_{1ax}=\alpha,\alpha_{1bx}=\alpha e^{-i\pi/4},\alpha_{2ay}=\alpha
e^{i\pi/2},$ and $\alpha_{2by}=\alpha e^{i\pi/4}$, these modifications result
in phase terms given by

$\displaystyle\delta S_{aa}=\delta S_{bb}$
$\displaystyle=\left(\frac{1}{2}\alpha\frac{\partial K}{\partial
r}-K\right)\Big{|}_{r=\alpha}t,$ (24) $\displaystyle\delta S_{ab}+\delta
S_{ba}$
$\displaystyle=2\left(\alpha^{2}\dot{\bar{\phi}}_{1x}-\bar{K}\right)\Big{|}_{r=r_{ab}}t.$
(25)

Here, $\delta S_{aa}$ and $\delta S_{bb}$ are evaluated along $r=\alpha$ which
is a constant in the first order classical canonical perturbation theory
approximation. The terms $\delta S_{ab}$ and $\delta S_{ba}$ can be evaluated
along the trajectory for $r$ in the ab component and we only use the
perturbation to $\phi_{x}$, and not $\phi_{y}$, due to the symmetries of the
initial conditions. The expressions for the entanglement entropy and fringe
visibility can be recovered from Eqs. (17) and (20), but now with the constant
$\frac{cg_{0}}{\alpha}$ replaced by

$C=2\left(\alpha^{2}\dot{\bar{\phi}}_{1x}-\bar{K}\right)\Big{|}_{r=r_{ab}}+\left(2K-\alpha\frac{\partial
K}{\partial r}\right)\Big{|}_{r=\alpha}.$ (26)

We calculate this quantity for a Coulomb potential in arbitrary dimension,
i.e.

$K=\begin{cases}-g_{0}\,r^{2-d}&d>2\\\ g_{0}\log{r}&d=2\end{cases},$ (27)

and for a Yukawa potential with

$K=g_{0}\frac{e^{-\mu r}}{r},$ (28)

where $\mu$ is related to the mass of the mediating particle. The values of
$C$ for the Yukawa potential with various values of the mass $\mu$ and for the
Coulomb potential in various dimensions $d$ are shown in Fig. 3. These
different values of the parameter $C$ would result in different scalings for
the entanglement entropy and fringe visibility. With precise enough visibility
measurements, such differences could serve as a probe for possible new quantum
modified gravity theories through their resulting particle entanglement.

(a) Yukawa

(b) Coulomb

Figure 3: Value of the constant $C$ given in the text which controls the
scaling of entanglement entropy and fringe visibility for (a) Yukawa potential
with various values of the Yukawa mass $\mu$, and (b) Coulomb potential with
dimension, $d=2,3$, as a function of the initial coherent state eigenvalue
$\alpha$, with $g_{0}=1$.

##  IV Conclusion

Merging quantum mechanics and gravity is one of the loftiest goals of modern
physics. Recent work investigating gravitationally induced quantum
entanglement provides a pathway to new insights into this problem.

In this paper, we have developed a thought experiment model to detect
gravitationally induced entanglement at the level of non-relativistic quantum
mechanics, similar to other proposals in the literature. Advantageously, our
setup readily lends itself to a more fundamental quantum field theoretic
description, to be presented in a subsequent work, wherein the massive
particles which become entangled through gravity are excitations of a quantum
scalar field, and the fringe visibility is modeled by an expectation value of
field observables.

This quantum field theoretic approach will allow us to go beyond the
predictions of our non-relativistic quantum mechanical model in several ways.
First of all, we will be able to systematically account for relativistic
corrections to our predicted measure of entanglement, corresponding to the
limit of a large harmonic oscillator trapping potential frequency. Secondly,
since different methods of carrying out this approach, such as the closed time
path formalism [51, 52], or the Dirac constraint quantization procedure [53,
45], treat gauge fixing in different ways, we may be able to better understand
more comprehensively the gauge aspects of gravitational entanglement
generation. Apart from these insights, there may be other, less expected,
lessons to be learned. As mentioned in the introduction, the description of
gravitational entanglement in terms of non-relativistic quantum mechanics is
not sensitive to the spin (i.e., helicity) of the mediating massless boson
[39]; a relativistic, quantum field theoretic approach that accounts for
retardation effects, such as the present proposed approach, is required in
order to reveal the distinguishing helicity two, graviton nature of the
entanglement.

## Acknowledgements

We thank Charis Anastopoulos, Markus Aspelmeyer, Bei-Lok Hu, David Mattingly,
Sougato Bose, Alexander Smith, Qidong Xu, and Shadi Ali Ahmad for very helpful
conversations This work was supported by the NSF under Grant No. PHY-2011382.

## References

  * Bose  _et al._ [2017] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. A. Geraci, P. F. Barker, M. Kim, and G. Milburn, Spin entanglement witness for quantum gravity, Phys. Rev. Lett. 119, 240401 (2017).
  * Marletto and Vedral [2017] C. Marletto and V. Vedral, Gravitationally-induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity, Phys. Rev. Lett. 119, 240402 (2017), arXiv:1707.06036 [quant-ph] .
  * Marshman  _et al._ [2020] R. J. Marshman, A. Mazumdar, and S. Bose, Locality and entanglement in table-top testing of the quantum nature of linearized gravity, Phys. Rev. A 101, 052110 (2020).
  * Marletto and Vedral [2019] C. Marletto and V. Vedral, Answers to a few questions regarding the BMV experiment, arXiv e-prints , arXiv:1907.08994 (2019), arXiv:1907.08994 [quant-ph] .
  * Marletto and Vedral [2018] C. Marletto and V. Vedral, When can gravity path-entangle two spatially superposed masses?, Phys. Rev. D 98, 046001 (2018).
  * Galley  _et al._ [2022] T. D. Galley, F. Giacomini, and J. H. Selby, A no-go theorem on the nature of the gravitational field beyond quantum theory, Quantum 6, 779 (2022), arXiv:2012.01441 [gr-qc, physics:quant-ph].
  * Marletto and Vedral [2020] C. Marletto and V. Vedral, Witnessing nonclassicality beyond quantum theory, Phys. Rev. D 102, 086012 (2020).
  * Christodoulou and Rovelli [2019] M. Christodoulou and C. Rovelli, On the possibility of laboratory evidence for quantum superposition of geometries, Physics Letters B 792, 64 (2019).
  * Danielson  _et al._ [2022] D. L. Danielson, G. Satishchandran, and R. M. Wald, Gravitationally mediated entanglement: Newtonian field versus gravitons, Phys. Rev. D 105, 086001 (2022).
  * Hall and Reginatto [2018] M. J. W. Hall and M. Reginatto, On two recent proposals for witnessing nonclassical gravity, Journal of Physics A Mathematical General 51, 085303 (2018), arXiv:1707.07974 [quant-ph] .
  * Reginatto and Hall [2018] M. Reginatto and M. J. W. Hall, Entanglement of quantum fields via classical gravity (2018), arXiv:1809.04989 [gr-qc, physics:quant-ph].
  * Rydving  _et al._ [2021] E. Rydving, E. Aurell, and I. Pikovski, Do Gedanken experiments compel quantization of gravity?, Physical Review D 104, 086024 (2021).
  * Anastopoulos and Hu [2018] C. Anastopoulos and B.-L. Hu, Comment on “a spin entanglement witness for quantum gravity” and on “gravitationally induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity”, arXiv preprint arXiv:1804.11315 (2018).
  * Anastopoulos  _et al._ [2021] C. Anastopoulos, M. Lagouvardos, and K. Savvidou, Gravitational effects in macroscopic quantum systems: a first-principles analysis, Classical and Quantum Gravity 38, 155012 (2021).
  * Ma  _et al._ [2022] Y. Ma, T. Guff, G. Morley, I. Pikovski, and M. S. Kim, Limits on inference of gravitational entanglement, Physical Review Research 4, 013024 (2022), arXiv:2111.00936 [quant-ph].
  * Fragkos  _et al._ [2022] V. Fragkos, M. Kopp, and I. Pikovski, On inference of quantization from gravitationally induced entanglement (2022), arXiv:2206.00558 [gr-qc, physics:quant-ph].
  * Altamirano  _et al._ [2018] N. Altamirano, P. Corona-Ugalde, R. B. Mann, and M. Zych, Gravity is not a pairwise local classical channel, Classical and Quantum Gravity 35, 145005 (2018).
  * Anastopoulos and Savvidou [2021] C. Anastopoulos and N. Savvidou, Quantum Information in Relativity: The Challenge of QFT Measurements, Entropy 24, 4 (2021).
  * Anastopoulos  _et al._ [2022] C. Anastopoulos, B.-L. Hu, and K. Savvidou, Quantum Field Theory based Quantum Information: Measurements and Correlations (2022), arXiv:2208.03696 [gr-qc, physics:hep-th, physics:quant-ph].
  * Anastopoulos and Hu [2022] C. Anastopoulos and B.-L. Hu, Gravity, quantum fields and quantum information: Problems with classical channel and stochastic theories, Entropy 24, 10.3390/e24040490 (2022).
  * Christodoulou  _et al._ [2022a] M. Christodoulou, A. Di Biagio, M. Aspelmeyer, v. Brukner, C. Rovelli, and R. Howl, Locally mediated entanglement through gravity from first principles (2022a), arXiv:2202.03368 [gr-qc, physics:quant-ph].
  * Christodoulou  _et al._ [2022b] M. Christodoulou, A. Di Biagio, R. Howl, and C. Rovelli, Gravity entanglement, quantum reference systems, degrees of freedom (2022b), arXiv:2207.03138 [gr-qc, physics:quant-ph].
  * Martín-Martínez and Perche [2022] E. Martín-Martínez and T. R. Perche, What gravity mediated entanglement can really tell us about quantum gravity (2022), arXiv:2208.09489 [gr-qc, physics:hep-th, physics:quant-ph].
  * DeLisle [2022] C. L. DeLisle, Quantum Information in Electromagnetism and Gravity (2022).
  * Nguyen and Bernards [2019] H. C. Nguyen and F. Bernards, Entanglement dynamics of two mesoscopic objects with gravitational interaction (2019), arXiv:1906.11184 [quant-ph] .
  * Tilly  _et al._ [2021] J. Tilly, R. J. Marshman, A. Mazumdar, and S. Bose, Qudits for witnessing quantum-gravity-induced entanglement of masses under decoherence, Phys. Rev. A 104, 052416 (2021).
  * Krisnanda  _et al._ [2019] T. Krisnanda, G. Y. Tham, M. Paternostro, and T. Paterek, Observable quantum entanglement due to gravity (2019), arXiv:1906.08808 [quant-ph] .
  * van de Kamp  _et al._ [2020] T. W. van de Kamp, R. J. Marshman, S. Bose, and A. Mazumdar, Quantum gravity witness via entanglement of masses: Casimir screening, Phys. Rev. A 102, 062807 (2020).
  * Matsumura and Yamamoto [2020] A. Matsumura and K. Yamamoto, Gravity-induced entanglement in optomechanical systems, Phys. Rev. D 102, 106021 (2020).
  * Chevalier  _et al._ [2020] H. Chevalier, A. J. Paige, and M. S. Kim, Witnessing the nonclassical nature of gravity in the presence of unknown interactions, Phys. Rev. A 102, 022428 (2020).
  * Carney  _et al._ [2021] D. Carney, H. Müller, and J. M. Taylor, Using an atom interferometer to infer gravitational entanglement generation, PRX Quantum 2, 030330 (2021).
  * Carlesso  _et al._ [2019] M. Carlesso, A. Bassi, M. Paternostro, and H. Ulbricht, Testing the gravitational field generated by a quantum superposition, New Journal of Physics 21, 093052 (2019).
  * Feng and Vedral [2022] T. Feng and V. Vedral, Amplification of Gravitationally Induced Entanglement (2022), arXiv:2202.09737 [quant-ph] .
  * Zhou  _et al._ [2022a] R. Zhou, R. J. Marshman, S. Bose, and A. Mazumdar, Catapulting towards massive and large spatial quantum superposition (2022a), arXiv:2206.04088 [quant-ph].
  * Marshman  _et al._ [2022] R. J. Marshman, A. Mazumdar, R. Folman, and S. Bose, Constructing nano-object quantum superpositions with a Stern-Gerlach interferometer, Physical Review Research 4, 023087 (2022).
  * Zhou  _et al._ [2022b] R. Zhou, R. J. Marshman, S. Bose, and A. Mazumdar, Gravito-diamagnetic forces for mass independent large spatial quantum superpositions (2022b), arXiv:2211.08435 [quant-ph].
  * Schut  _et al._ [2022] M. Schut, J. Tilly, R. J. Marshman, S. Bose, and A. Mazumdar, Improving resilience of quantum-gravity-induced entanglement of masses to decoherence using three superpositions, Physical Review A 105, 032411 (2022).
  * Guff  _et al._ [2022] T. Guff, N. Boulle, and I. Pikovski, Optimal fidelity witnesses for gravitational entanglement, Physical Review A 105, 022444 (2022).
  * Carney [2022] D. Carney, Newton, entanglement, and the graviton, Phys. Rev. D 105, 024029 (2022).
  * Cioslowski and Pernal [2000] J. Cioslowski and K. Pernal, The ground state of harmonium, The Journal of Chemical Physics 113, 8434 (2000).
  * Delić  _et al._ [2020] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367, 892 (2020), https://www.science.org/doi/pdf/10.1126/science.aba3993 .
  * Yant and Blencowe [2023] J. Yant and M. P. Blencowe (2023), in preparation.
  * Donoghue [1994] J. F. Donoghue, General relativity as an effective field theory: The leading quantum corrections, Phys. Rev. D 50, 3874 (1994).
  * Arteaga  _et al._ [2004] D. Arteaga, R. Parentani, and E. Verdaguer, Propagation in a thermal graviton background, Phys. Rev. D 70, 044019 (2004).
  * Oniga and Wang [2017] T. Oniga and C. H.-T. Wang, Quantum coherence, radiance, and resistance of gravitational systems, Physical Review D 96, 084014 (2017).
  * Xu and Blencowe [2022] Q. Xu and M. P. Blencowe, Zero-dimensional models for gravitational and scalar QED decoherence, New Journal of Physics 24, 113048 (2022), arXiv:2005.02554 [gr-qc, physics:quant-ph].
  * Heller [1975] E. J. Heller, Time‐dependent approach to semiclassical dynamics, The Journal of Chemical Physics 62, 1544 (1975).
  * Maleki and Maleki [2022] Y. Maleki and A. Maleki, Complementarity-entanglement tradeoff in quantum gravity, Physical Review D 105, 10.1103/physrevd.105.086024 (2022).
  * Barker  _et al._ [2022] P. Barker, S. Bose, R. J. Marshman, and A. Mazumdar, Entanglement based tomography to probe new macroscopic forces, Physical Review D 106, L041901 (2022).
  * Mostepanenko and Klimchitskaya [2020] V. M. Mostepanenko and G. L. Klimchitskaya, The State of the Art in Constraining Axion-to-Nucleon Coupling and Non-Newtonian Gravity from Laboratory Experiments, Universe 6, 147 (2020).
  * Campos and Hu [1998] A. Campos and B. Hu, Nonequilibrium dynamics of a thermal plasma in a gravitational field, Phys. Rev. D 58, 125021 (1998).
  * Calzetta and Hu [2008] E. A. Calzetta and B.-L. Hu,  _Nonequilibrium Quantum Field Theory_ (Cambridge University Press, Cambridge, England, 2008).
  * Oniga and Wang [2016] T. Oniga and C. H.-T. Wang, Quantum gravitational decoherence of light and matter, Phys. Rev. D 93, 044027 (2016).
  * Demarie [2012] T. F. Demarie, Pedagogical introduction to the entropy of entanglement for Gaussian states (2012), arXiv:1209.2748 [quant-ph].
  * Colavita and Hacyan [2005] E. Colavita and S. Hacyan, Coherent states with elliptical polarization, Physics Letters A 337, 183 (2005).

##  Appendix A Approximate Evolution of Gaussian Parameters

In order to find the solutions to Eqs. (12)-(14), restated here,

$\displaystyle\expectationvalue{\dot{\bm{r}}}$ $\displaystyle=\frac{\partial
H}{\partial\expectationvalue{\bm{p}_{r}}},\quad\expectationvalue{\dot{\bm{p}_{r}}}=-\frac{\partial
H}{\partial\expectationvalue{\bm{r}}},$ (29) $\displaystyle\dot{\Omega}$
$\displaystyle=-2\Omega^{2}-\frac{1}{2}\frac{\partial^{2}V}{\partial
r_{i}\partial r_{j}}\Big{|}_{\bm{r}=\expectationvalue{\bm{r}}},$ (30)
$\displaystyle\dot{\gamma}$
$\displaystyle=i\,\mathrm{Tr}\,\Omega+\mathcal{L}.$ (31)

we introduce a bookkeeping parameter $\lambda$ into the Hamiltonian (9),
multiplying the gravitational potential term

$H=\frac{1}{2}\left(p_{r}^{2}+r^{2}-\lambda\frac{g_{0}}{r}\right).$ (32)

Thus we treat the gravitational term as a perturbation to the harmonic
oscillator, and at the end of the calculation take $\lambda\rightarrow 1$;
this procedure is validated a posteriori with the relative magnitudes of the
terms in the resulting perturbation series being sufficiently small. For
$\expectationvalue{\bm{r}}$ and $\expectationvalue{\bm{p}_{r}}$, we use
classical canonical perturbation theory, moving to the Hamilton-Jacobi
coordinates of the 2D harmonic oscillator given by

$\expectationvalue{r_{x}}=\sqrt{2\eta_{x}}\sin\left(t+\phi_{x}\right),\
\expectationvalue{r_{y}}=\sqrt{2\eta_{y}}\sin\left(t+\phi_{y}\right).$ (33)

The Hamilton-Jacobi coordinates are those generated by a canonical
transformation where the unperturbed Hamiltonian is zero and therefore these
coordinates are constants. However, the gravitational term remains and so the
transformed Hamiltonian (Kameltonian) in these coordinates is given by

$K=-\lambda\frac{g_{0}}{\sqrt{2}}[\eta_{x}\sin^{2}\left(t+\phi_{x}\right)+\eta_{y}\sin^{2}\left(t+\phi_{y}\right)]^{-1/2}.$
(34)

We assume a perturbative solution and to first order the equations of motion
are given by

$\displaystyle\dot{\eta}_{x1}$ $\displaystyle=-\frac{\partial
K}{\partial\phi_{x}}\Big{|}_{0}=-\lambda
g_{0}\frac{\eta_{x0}\sin\left(t+\phi_{x0}\right)\cos\left(t+\phi_{x0}\right)}{\sqrt{2}[\eta_{x0}\sin^{2}\left(t+\phi_{x0}\right)+\eta_{y0}\sin^{2}\left(t+\phi_{y0}\right)]^{3/2}},$
(35) $\displaystyle\dot{\phi}_{x1}$ $\displaystyle=\frac{\partial
K}{\partial\eta_{x}}\Big{|}_{0}=\lambda
g_{0}\frac{\sin^{2}\left(t+\phi_{x0}\right)}{2\sqrt{2}[\eta_{x0}\sin^{2}\left(t+\phi_{x0}\right)+\eta_{y0}\sin^{2}\left(t+\phi_{y0}\right)]^{3/2}},$
(36)

and similarly for the $y$ components. The right hand side of these equations
is periodic and we approximate the solutions by their time averages over an
oscillator period as follows:

$\displaystyle\eta_{x1}$
$\displaystyle\approx\dot{\overline{\eta}}_{x1}t=-\lambda\frac{g_{0}t}{2\sqrt{2}\pi}\int_{0}^{2\pi}\frac{\eta_{x0}\sin\left(t+\phi_{x0}\right)\cos\left(t+\phi_{x0}\right)}{[\eta_{x0}\sin^{2}\left(t+\phi_{x0}\right)+\eta_{y0}\sin^{2}\left(t+\phi_{y0}\right)]^{3/2}}dt,$
(37) $\displaystyle\phi_{x1}$
$\displaystyle\approx\dot{\overline{\phi}}_{x1}t=\lambda\frac{g_{0}t}{4\sqrt{2}\pi}\int_{0}^{2\pi}\frac{\sin\left(t+\phi_{x0}\right)^{2}}{[\eta_{x0}\sin^{2}\left(t+\phi_{x0}\right)+\eta_{y0}\sin^{2}\left(t+\phi_{y0}\right)]^{3/2}}dt,$
(38)

and similarly for the $y$ components. For initial conditions given by
$\alpha_{1}=(\alpha,0)$ and $\alpha_{2}=(0,\alpha e^{i\pi/2})$ for some
$\alpha>0$, the relative difference coordinate coherent state parameters are
$\alpha_{rx}=\frac{\alpha}{\sqrt{2}}$ and
$\alpha_{ry}=-i\frac{\alpha}{\sqrt{2}}$. This leads to the following first
order solutions:

$\expectationvalue{r_{x}}=\alpha\cos\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right],\
\expectationvalue{r_{y}}=-\alpha\sin\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right].$
(39)

Such a time average has the advantage that $\expectationvalue{r}=\alpha$ is
constant.

We can only approximate $\Omega$ analytically for these particular initial
conditions where $\expectationvalue{r}=\alpha$ is constant. Assuming a
solution of the form $\Omega\approx\Omega_{0}+\lambda\Omega_{1}$ to Eq. (13)
for the initial conditions $\alpha_{1}=(\alpha,0)$ and $\alpha_{2}=(0,\alpha
e^{i\pi/2})$ this results in the following equations of motion:

$\displaystyle\dot{\Omega}_{0}=$ $\displaystyle-2\Omega_{0}^{2}-\mathbb{1},$
(40) $\displaystyle\dot{\Omega}_{1}=$
$\displaystyle-2(\Omega_{0}\Omega_{1}+\Omega_{1}\Omega_{0})$
$\displaystyle-\frac{g_{0}}{2\alpha^{3}}\begin{bmatrix}1-3\cos^{2}\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]&&3\cos\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]\sin\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]\\\
3\cos\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]\sin\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]&&1-3\sin^{2}\left[\left(1+\lambda\frac{g_{0}}{2\alpha^{3}}\right)t\right]\end{bmatrix}.$
(41)

The initial condition for all coherent state parameters is given by
$\Omega=\frac{i}{2}\mathbb{1}$. This gives a solution to the equation for
$\Omega_{1}$, which we then expand as a series in $\lambda$, but keep terms
where it appears in the phases of oscillating terms. We find that the solution
has a lowest order term proportional to $\lambda^{-1}$. This is due to the
fact that we are taking an ‘improper’ series expansion by leaving in the phase
terms. However if we assume an approximate solution for $\Omega$ given by
$\Omega_{0}$ and the $\lambda^{-1}$ order term of $\Omega_{1}$ we find that
this is also a solution of the zeroth order equation of motion:

$\Omega=\frac{i}{2}\mathbb{1}+\frac{3}{8}\left[e^{-i\left(2+\frac{g_{0}}{\alpha^{3}}\right)t}-e^{-2it}\right]\begin{bmatrix}i&&1\\\
1&&-i\end{bmatrix},$ (42)

where we have taken $\lambda\rightarrow 1$ and note that this expression is
still valid for $\alpha\gg 1$. The off-diagonal terms clearly make the state
nonseparable in the $r_{x}$ and $r_{y}$ coordinates, and upon transforming
back to the particle coordinates, the scaled covariance matrix which we will
denote $\Omega_{12}$ is given by

$\Omega_{12}=\frac{i}{2}\mathbb{1}+\frac{3}{16}\left[e^{-i\left(2+\lambda\frac{g_{0}}{\alpha^{3}}\right)t}-e^{-2it}\right]\begin{bmatrix}i&&1&&i&&-1\\\
1&&-i&&-1&&i\\\ i&&-1&&i&&1\\\ -1&&i&&1&&-i\end{bmatrix},$ (43)

which implies that this is an entangled state between particles 1 and 2.

Thus, gravity gives rise to entanglement without considering a superposition
of coherent state products. In order to quantify this entanglement, we
evaluate the entanglement entropy using the formalism for Gaussian states
[54]. The Wigner function for this state is given by a Gaussian distribution
on phase space whose covariance matrix, in the two particle phase space
coordinates basis, is given by
$\left(\bm{x}_{1},\bm{x}_{2},\bm{p}_{1},\bm{p}_{2}\right)$, is

$\Sigma=\begin{bmatrix}\frac{1}{2}\left(\mathrm{Im}\,\Omega\right)^{-1}&&-\left(\mathrm{Im}\,\Omega\right)^{-1}\mathrm{Re}\,\Omega\\\
Re\,\Omega\left(\mathrm{Im}\,\Omega\right)^{-1}&&2\left[\mathrm{Im}\,\Omega+\mathrm{Re}\,\Omega\left(\mathrm{Im}\,\Omega\right)^{-1}\mathrm{Re}\,\Omega\right]\end{bmatrix}.$
(44)

The symplectic eigenvalues, $\nu_{i}$, of one of the subsystems are obtained
by determining the eigenvalues of $iJ\Sigma_{1}$, where $\Sigma_{1}$ is the
submatrix of the covariance matrix corresponding to particle 1 and $J$ is the
symplectic form given by

$J=\begin{bmatrix}0&&\mathbb{1}\\\ -\mathbb{1}&&0\end{bmatrix}.$ (45)

The entanglement entropy is then given by

$\displaystyle S=$
$\displaystyle\sum_{i}\frac{1}{2}\left(\nu_{i}+1\right)\log\left[\frac{1}{2}\left(\nu_{i}+1\right)\right]-\frac{1}{2}\left(\nu_{i}-1\right)\log\left[\frac{1}{2}\left(\nu_{i}-1\right)\right]$
$\displaystyle\approx\frac{9g_{0}^{2}t^{2}}{64\alpha^{6}}\left[1-\log\left(\frac{9g_{0}^{2}t^{2}}{128\alpha^{6}}\right)\right],$
(46)

where we sum over the symplectic eigenvalues for one of the particles. The
bottom expression here contains the lowest order terms in $\alpha^{-1}$.

##  Appendix B Approximate Phase Dependent Entanglement Entropies

For initial conditions given by $\alpha_{1ax}=\alpha,\alpha_{1bx}=\alpha
e^{-i\pi/4},\alpha_{2ay}=\alpha e^{i\pi/2},$ and $\alpha_{2by}=\alpha
e^{i\pi/4}$, we find the terms $\delta S$ for the various components by using
the classical canonical perturbation theory method described in appendix A to
solve for the classical paths. We then use these trajectories to calculate the
difference in action between the perturbed and unperturbed (pure harmonic
oscillator) solutions. We also take $\lambda\rightarrow 1$ here and find that
our approximation are still valid for $\alpha\gg 1$. The $\delta S$ terms are
given by

$\displaystyle\delta S_{aa}=$ $\displaystyle\delta
S_{bb}=\frac{3g_{0}}{2\alpha}t$ (47) $\displaystyle\delta S_{ab}=$
$\displaystyle\frac{\alpha^{2}}{4}\left\\{\cos\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]-\sin\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]-\cos\left(2t\right)+\sin\left(2t\right)\right\\}$
$\displaystyle+\frac{g_{0}}{\alpha}t\left(C_{\phi}+\frac{1}{2}C_{\lambda}\left\\{\cos\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]+\sin\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]\right\\}+C_{gab}\right)$
(48) $\displaystyle\delta S_{ba}=$
$\displaystyle-\frac{\alpha^{2}}{4}\left\\{\cos\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]-\sin\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]-\cos\left(2t\right)+\sin\left(2t\right)\right\\}$
$\displaystyle+\frac{g_{0}}{\alpha}t\left(C_{\phi}-\frac{1}{2}C_{\lambda}\left\\{\cos\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]+\sin\left[2\left(1+C_{\phi}\frac{g_{0}}{\alpha^{3}}\right)t\right]\right\\}+C_{gba}\right),$
(49)

where $C_{\phi},C_{\lambda},C_{gab},$ and $C_{gba}$ are given respectively by

$\displaystyle C_{\phi}=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}dt\frac{\cos^{2}\left(t+\frac{\pi}{4}\right)}{\left(\sin^{2}(t)+\cos^{2}\left(t+\frac{\pi}{4}\right)\right)^{3/2}}$
(50) $\displaystyle C_{\lambda}=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}dt\frac{\sin(t)\cos(t)}{\left(\sin^{2}(t)+\cos^{2}\left(t+\frac{\pi}{4}\right)\right)^{3/2}}$
(51) $\displaystyle C_{gab}=$
$\displaystyle\frac{1+C_{\phi}\frac{g_{0}}{\alpha^{3}}}{\pi}\int_{0}^{\frac{\pi}{1+C_{\phi}\frac{g_{0}}{\alpha^{3}}}}dt\left(\sin^{2}\left(t+\frac{\pi}{4}\right)+\cos^{2}(t)\right)^{-1/2}$
(52) $\displaystyle C_{gba}=$
$\displaystyle\frac{1+C_{\phi}\frac{g_{0}}{\alpha^{3}}}{\pi}\int_{0}^{\frac{\pi}{1+C_{\phi}\frac{g_{0}}{\alpha^{3}}}}dt\left(\sin^{2}(t)+\cos^{2}\left(t+\frac{\pi}{4}\right)\right)^{-1/2}.$
(53)

The reduced density matrix for particle 1 is given by

$\begin{split}\rho_{1}=&|N|^{2}\left\\{\ket{\bm{\alpha}_{1a}}\bra{\bm{\alpha}_{1a}}\left[2+2e^{\alpha^{2}\left(\frac{1}{\sqrt{2}}-1\right)}\cos\left(\delta
S_{aa}-\delta S_{ab}+\frac{\alpha^{2}}{\sqrt{2}}\right)\right]\right.\\\
&+\ket{\bm{\alpha}_{1a}}\bra{\bm{\alpha}_{1b}}\left[e^{i\left(\delta
S_{ab}-\delta S_{bb}\right)}+e^{i\left(\delta S_{aa}-\delta
S_{ba}\right)}\right.\\\
&\left.+e^{\alpha^{2}\left(\frac{1}{\sqrt{2}}-1\right)}\left(e^{i\left(\delta
S_{aa}-\delta S_{bb}+\frac{\alpha^{2}}{\sqrt{2}}\right)}+e^{i\left(\delta
S_{ab}-\delta S_{ba}-\frac{\alpha^{2}}{\sqrt{2}}\right)}\right)\right]\\\
&+\ket{\bm{\alpha}_{1b}}\bra{\bm{\alpha}_{1a}}\left[e^{-i\left(\delta
S_{ab}-\delta S_{bb}\right)}+e^{-i\left(\delta S_{aa}-\delta
S_{ba}\right)}\right.\\\
&\left.+e^{\alpha^{2}\left(\frac{1}{\sqrt{2}}-1\right)}\left(e^{-i\left(\delta
S_{aa}-\delta S_{bb}+\frac{\alpha^{2}}{\sqrt{2}}\right)}+e^{-i\left(\delta
S_{ab}-\delta S_{ba}-\frac{\alpha^{2}}{\sqrt{2}}\right)}\right)\right]\\\
&\left.\ket{\bm{\alpha}_{1b}}\bra{\bm{\alpha}_{1b}}\left[2+2e^{\alpha^{2}\left(\frac{1}{\sqrt{2}}-1\right)}\cos\left(\delta
S_{ba}-\delta
S_{bb}+\frac{\alpha^{2}}{\sqrt{2}}\right)\right]\right\\}.\end{split}$ (54)

We shall denote the four outer product terms in the curly brackets of this
expression respectively as $\rho_{aa}$, $\rho_{ab}$, $\rho_{ba}$, and
$\rho_{bb}$. The eigenvalues of the density matrix $\sigma_{i}$ are given by

$\begin{split}\sigma_{i}=&\frac{|N|^{2}}{2}\left[\rho_{aa}+\rho_{bb}+2\mathrm{Re}\,\left(e^{\alpha^{2}\left(e^{i\pi/4}-1\right)}\rho_{ab}\right)\pm\left(\left(\rho_{aa}-\rho_{bb}\right)^{2}+4\left\\{e^{\alpha^{2}(\sqrt{2}-2)}\rho_{aa}\rho_{bb}+|\rho_{ab}|^{2}\right.\right.\right.\\\
&\left.\left.\left.+\left(\rho_{aa}+\rho_{bb}\right)\mathrm{Re}\,\left(e^{\alpha^{2}(e^{i\pi/4}-1)}\rho_{ba}\right)-\left[\mathrm{Im}\,\left(e^{\alpha^{2}(e^{i\pi/4}-1)}\rho_{ba}\right)\right]^{2}\right\\}\right)^{1/2}\right],\end{split}$
(55)

where $i=\\{1,2\\}$. The coherent state components of particle 1 are
orthogonal to a good approximation for $\alpha\gg 1$. In this limit, the
reduced state for particle 1 can be written as

$\begin{split}\rho_{1}\approx&\frac{1}{4}\\{2\ket{\bm{\alpha}_{1a}}\bra{\bm{\alpha}_{1a}}+\left(e^{i(\delta
S_{ab}-\delta S_{bb})}+e^{i(\delta S_{aa}-\delta
S_{ba})}\right)\ket{\bm{\alpha}_{1a}}\bra{\bm{\alpha}_{1b}}\\\
&+\left(e^{-i\left(\delta S_{ab}-\delta S_{bb}\right)}+e^{-i\left(\delta
S_{aa}-\delta
S_{ba}\right)}\right)\ket{\bm{\alpha}_{1b}}\bra{\bm{\alpha}_{1a}}+2\ket{\bm{\alpha}_{1b}}\bra{\bm{\alpha}_{1b}}\\},\end{split}$
(56)

and the eigenvalues can be approximated as

$\sigma_{i}\approx\frac{1}{2}\pm\frac{1}{4}|\rho_{ab}|,$ (57)

with

$\begin{split}|\rho_{ab}|\approx&\left[2+2\cos\left(\delta S_{ab}-\delta
S_{bb}-\delta S_{aa}+\delta
S_{ba}\right)\right]^{1/2}=\left\\{2+2\cos\left[\left(2C_{\phi}+C_{gab}+C_{gba}-3\right)\left(\frac{g_{0}t}{\alpha}\right)\right]\right\\}^{1/2}\\\
\approx&2-\frac{1}{4}\left(2C_{\phi}+C_{gab}+C_{gba}-3\right)^{2}\left(\frac{g_{0}t}{\alpha}\right)^{2}.\end{split}$
(58)

The entanglement entropy is then approximately

$\begin{split}S=&-\sigma_{1}\log\sigma_{1}-\sigma_{2}\log\sigma_{2}\\\
\approx&\frac{1}{16}\left(2C_{\phi}+C_{gab}+C_{gba}-3\right)^{2}\left(\frac{g_{0}t}{\alpha}\right)^{2}\\\
&\left\\{1-\log\left[\frac{1}{16}\left(2C_{\phi}+C_{gab}+C_{gba}-3\right)^{2}\left(\frac{g_{0}t}{\alpha}\right)^{2}\right]\right\\}.\end{split}$
(59)

This entanglement grows faster than that due to the covariance terms in one
component of our full state by a factor of $\alpha^{4}$, which will be very
large in any realistic implementation, so our approximation where the local
entanglement terms are ignored is justified. Noting this, we can further
simplify the above expression by approximating $C_{gab}$, and $C_{gba}$ to
find that they only have corrections proportional to $\alpha^{-3}$ which are
negligible. Therefore, we have

$\displaystyle C_{gab}\approx$
$\displaystyle\,C_{gba}\approx\frac{1}{2\pi}\int_{0}^{2\pi}dt\left\\{\sin^{2}\left(t+\frac{\pi}{4}\right)+\cos^{2}t\right\\}^{-1/2},$
(60)

giving the following expression for the entanglement entropy:

$\displaystyle S\approx$
$\displaystyle\left(\frac{cg_{0}}{\alpha}t\right)^{2}\left\\{1-\log\left[\left(\frac{cg_{0}}{\alpha}t\right)^{2}\right]\right\\},$
(61)

where $c=\frac{1}{4}\left(2C_{\phi}+2C_{gab}-3\right)\approx 0.1004$.

##  Appendix C Numerical Validation of Analytical Approximation to Time
Evolution

In order to validate the analytical approximation to the time evolved state of
the gravitational harmonium described above, we also numerically approximate
the time evolution of this state and compare the numerical and analytical
approximations. In order to do this we choose a basis for our state space
given by simultaneous eigenstates of the two dimensional harmonic oscillator
Hamiltonian and the angular momentum operator, whose position space wave
functions are given by

$\bra{r,\phi}\ket{n,\pm
l}=\left(\frac{n!}{\pi\left(n+l\right)!}\right)^{1/2}e^{\pm
il\phi}e^{-r^{2}/2}r^{l}L_{n}^{l}(r^{2})$ (62)

where $l>0$, and $L_{n}^{l}$ are the generalized Laguerre polynomials [55].
The energy of these states is $E_{n,l}=2n+l+1$ for a simple harmonic
oscillator, in the absence of gravity.

The gravitational potential term can be evaluated in this basis with matrix
elements given by

$\displaystyle V_{nmjl}=$
$\displaystyle\bra{n,l}-\frac{g_{0}}{r}\ket{m,j}=-2g_{0}\delta_{jl}\left(\frac{n!m!}{(n+l)!(m+j)!}\right)^{1/2}\int_{0}^{\infty}e^{-r^{2}}r^{j+l}L_{n}^{l}(r^{2})L_{m}^{j}(r^{2})dr$
(63) $\displaystyle=$ $\displaystyle-
g_{0}\sum_{i=0}^{N}\frac{\Gamma(N+l+1/2)x_{i}}{N!(N+1)^{2}(L_{N+1}^{l-1/2}(x_{i}))^{2}}L_{n}^{l}(x_{i})L_{m}^{l}(x_{i}),$

where $N>n,m$, the $x_{i}$ are the zeros of $L_{N}^{l-1/2}$, and the last
equality is found by Gauss-Laguerre quadrature. Therefore the Hamiltonian for
the gravitational harmonium has matrix elements in this basis given by
$H_{nmlj}=\left(2n+l+1\right)\delta_{nm}\delta_{jl}+V_{nmjl}$. By evaluating a
finite number of these basis elements, we approximate the Hamiltonian on the
subspace spanned by the finite number of basis states corresponding to these
matrix elements. We then numerically diagonalize this matrix and approximate
the propagator on this subspace by $U(t)=Pe^{-iEt}P^{\dagger}$, where $P$ is
the change of basis matrix to the eigenbasis of $H$, and $E$ is the diagonal
matrix of eigenvalues of $H$.

We act with this propagator on an initial state corresponding to one component
of our initial state considered above, i.e., a product of coherent states of
particles 1 and 2. This corresponds to a coherent state of a 2D harmonic
oscillator in the separation coordinate we are considering here. In this
basis, coherent states are given by

$\ket{\alpha,\beta}_{\pm}=\sum_{n,l=0}^{\infty}\left(\frac{(n+l)!}{n!}\right)^{1/2}\frac{1}{l!}\alpha^{l}\beta^{n}\ket{n,\pm
l}.$ (64)

The parameters $\alpha,\beta$ correspond roughly to the semi-major axis and
eccentricity of the trajectory of the position expectation value respectively,
and the label $\pm$ corresponds to counter-clockwise or clockwise rotation
[55]. For the component of our state corresponding to particles 1 and 2 having
constant separation between position expectation values, we find $\beta=0$ and
$\alpha$ is the same as considered in previous sections.

In order to determine if the numerically approximated time evolved state
validates our analytical approximation, we make use of a distance measure
between states given by

$d(\ket{\psi_{1}},\ket{\psi_{2}})=1-\mathrm{Re}\,\bra{\psi_{1}}\ket{\psi_{2}}$
(65)

and evaluate the distance between our numerical approximation and either the
analytical approximation of the gravitational harmonium, or a pure harmonic
oscillator without gravity, both projected onto the subspace spanned by the
basis elements we used here, i.e., $d_{GH}\left(U(t)\ket{\alpha},e^{i\delta
S}\ket{\alpha(t)}\right)$ and
$d_{SHO}\left(U(t)\ket{\alpha},\ket{\alpha(t)}\right)$, respectively. We find
that for $1<\alpha<10$ and $\frac{g_{0}}{\alpha}=10^{-6}$, and $t$ up to $1$s
(which corresponds to $t=10^{5}$ in units $\omega^{-1}$) as in Sec. II, we
have $d_{GH}<d_{SHO}$ (see Fig. 4). This establishes that our time evolved
state is closer to that of our analytical approximation of the gravitational
harmonium than to the evolving harmonic oscillator state with gravity
neglected, thus validating our approximation.

Figure 4: Distance measures $d_{SHO}$ and $d_{GH}$ for various values of the
parameters $g_{0}$ and $\alpha$, with $\frac{g_{0}}{\alpha}=10^{-6}$ as a
function of time (in units $\omega^{-1}$). With $d_{SHO}>d_{GH}$, we have that
the numerical approximation is closer to our analytical approximation than to
a pure harmonic oscillator as quantified by this distance measure.


# Targetless Extrinsic Calibration of Stereo Cameras, Thermal Cameras, and
Laser Sensors in the Wild

Taimeng Fu, Huai Yu, Wen Yang, Yaoyu Hu and Sebastian Scherer This work has
been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be
accessible.Taimeng Fu is with the School of Electronic Information, Wuhan
University, Wuhan 430072, China, and also with the School of Data Science, The
Chinese University of Hong Kong, Shenzhen, 2001 Longxiang Boulevard, Longgang
District, Shenzhen, China, 518712. E-mail<EMAIL_ADDRESS>Yu
and Wen Yang are with the School of Electronic Information, Wuhan University,
Wuhan 430072, China. E-mail<EMAIL_ADDRESS>Hu and Sebastian
Scherer are with the AirLab, Carnagie Mellon University, Pittsburgh, PA 15213
E-mail: {yaoyuh<EMAIL_ADDRESS>
###### Abstract

The fusion of multi-modal sensors has become increasingly popular in
autonomous driving and intelligent robots since it can provide richer
information than any single sensor, enhance reliability in complex
environments. Multi-sensor extrinsic calibration is one of the key factors of
sensor fusion. However, such calibration is difficult due to the variety of
sensor modalities and the requirement of calibration targets and human labor.
In this paper, we demonstrate a new targetless cross-modal calibration
framework by focusing on the extrinsic transformations among stereo cameras,
thermal cameras, and laser sensors. Specifically, the calibration between
stereo and laser is conducted in 3D space by minimizing the registration
error, while the thermal extrinsic to the other two sensors is estimated by
optimizing the alignment of the edge features. Our method requires no
dedicated targets and performs the multi-sensor calibration in a single shot
without human interaction. Experimental results show that the calibration
framework is accurate and applicable in general scenes.

##  I INTRODUCTION

Extrinsic calibration between multi-modal sensors is a fundamental task for a
variety of computer vision and robot applications, such as visual sensor
fusion [1] and object fusion [2], 3D reconstruction [3], and SLAM [4]. The
calibration process is to estimate the rigid-body transformation between each
two sensors, which serves as a bridge to connect sensors, providing fusion
information for the perception system, and enhancing the robot’s ability to
perceive the environment. Over the past several years, substantial works have
been proposed to solve this problem either by manual calibration or using
specific calibration targets, such as checkerboard and AprilTag. However,
these methods are suffering from limited flexibility and are inaccurate on
dynamic systems in which the extrinsic may slightly change after the
calibration. In this paper, we mainly focus on the automatic extrinsic
calibration between stereo cameras, thermal cameras, and laser sensors in the
wild without any specific target. It has the great potential to overcome these
limitations of existing calibration methods.

Figure 1: Thermal mapping on a dense reconstructed point cloud with the
thermal extrinsic calibrated by our method.

We will review the recent works on the calibration of monocular/stereo visual
cameras to laser sensors and visual to thermal cameras, then propose our
contributions beyond these works.

###  I-A Calibration of monocular camera and laser sensor

Existing literature on the calibration of monocular cameras and laser sensors
can be classified according to 1) whether they require a specific calibration
target (e.g. checkerboard), 2) whether they require manually pairing the laser
points and camera pixels. Scaramuzza _et al._ [5] demonstrated a laser point
cloud visualization technique that enables manually labeling the
correspondences between laser points and image pixels. They employed the
perspective-from-n-points (PnP) algorithm to calculate the transformation
based on these matches. Their method does not depend on fixed calibration
targets, however, the manual matching is time-consuming and inflexible. On the
other side, Nunez _et al._ [6] chose to increase the degree of automation with
a fixed checkerboard. They developed an algorithm to automatically detect the
checkerboard in laser’s view and align it to camera image to get the extrinsic
parameters. To improve the accuracy, multiple laser frames were aggregated
with the help of an inertial measurement unit (IMU). Although their
calibration process does not require much manual intervention, the usage of
checkerboard still limits its flexibility, since the specific calibration
scene is not always available.

There are also some attempts to get rid of both fixed scene and manual
manipulation. Pandey _et al._ [7] provided a solution of automatic targetless
extrinsic calibration by maximizing mutual information (MI) of laser
reflectivity and image intensity. They found that after aggregating 10 scans,
their MI cost function becomes convex and thus easy to optimize. However, the
laser reflectivity and image intensity might not be strongly correlated in
some scenarios. This might result in inaccurate calibrations. Levinson and
Thrun [8] developed an online system that automatically aligns laser edges to
image edges to correct sensor drifting. They projected laser edge points onto
the generated cost maps to get the edge alignment cost and tried to adjust the
extrinsic parameters to lower it. However, since the laser scans are
relatively sparse, some of the laser edge points may not precisely lie on the
boundary of the objects, which affects the calibration accuracy. Besides, they
used a greedy approach to reduce the cost, which is less efficient than a
gradient-based optimizer.

###  I-B Calibration of stereo cameras and laser sensor

The stereo-laser calibration is more straightforward than that of monocular
cameras and lasers since both stereo cameras and lasers have depth
information. Guindel _et al._ [9] calibrated stereo cameras and lasers based
on a four-hole calibration board. They developed a segmentation pipeline for
extracting boundaries of the holes from stereo point cloud and laser point
cloud respectively and calculated the extrinsic parameters by minimizing the
registration distance between the clustered centroids of the four holes. Dhall
_et al._ [10] reported a method that relies on a calibration board with a
visual tag on it. They obtained the 3D position of the board corners in both
stereo and laser’s view, and solve a set of equations to get the
transformation that minimizes the distance between corresponding corner
points. However, these methods still need known calibration targets to
formulate the easy feature correspondences, which have the same problem of
limited flexibility.

###  I-C Calibration of visual and thermal cameras

In recent years, with the development and wide use of thermal cameras, visual-
thermal calibration becomes an unavoidable task for image fusion and object
fusion [1, 2]. Most of the attempts tried to build calibration targets that
are distinguishable in both visual and thermal views. Li _et al._ [11] placed
LED bulbs on their calibration board. The bulbs have high intensity on both
RGB and thermal images since they emit both light and heat. They developed an
algorithm to localize the bulbs on images and calibrate the extrinsic
parameters by minimizing the reprojection error. Shivakumar _et al._ [12]
mounted aluminum squares on black acrylic background to form a checkerboard.
Since the checkerboard blocks have different color and thermal reflectivity,
it is easy to recognize the pattern on both RGB and thermal images. They then
employed OpenCV’s [13] camera calibration toolbox to estimate the extrinsic
parameters.

###  I-D Our contribution framework

Most of the calibration methods mentioned above require specific calibration
targets. This configuration limits their flexibility and may results in
inaccuracy on dynamic systems with unpredictable extrinsic drifts after the
calibration. On the other side, the scene-independent ones either require
manual matching or are not robust enough. Besides, all of them only focus on
two-modal calibration. It’s inconvenient to use them sequentially when types
of sensors increase.

To solve these problems, we propose an automatic, targetless, all-in-one
calibration framework for stereo cameras, thermal cameras, and laser sensors
in the wild. For the stereo-laser calibration, we develop a multi-frame ICP to
register stereo and laser point clouds. To calibrate the thermal extrinsic, we
optimize the alignment of the edge features in stereo and laser point clouds
and thermal images. Since the edge features are commonly available, our method
works in arbitrary environments. Besides, it calibrates the three sensors with
a single round of data collection. This greatly simplifies the cross-modal
multi-sensor calibration process.

##  II METHODOLOGY

The goal of our algorithm is to take a series of $n$ synchronized stereo image
pairs $I^{left}_{1:n},I^{right}_{1:n}$, thermal images $I^{thermal}_{1:n}$,
and laser point clouds $C^{laser}_{1:n}$, captured in arbitrary scenes, and
automatically optimize the initial guess of the 6-DoF rigid-body
transformations to get accurate calibrations. The transformation is defined by
six parameters $\xi=\\{x,y,z,row,pitch,yaw\\}$, where $x,y,z$ are
translations, and $roll,pitch,yaw$ are Euler angle rotations. We take the
stereo left camera’s coordinate system $\hat{S}$ as the base coordinate
system, and calibrate other sensors to $\hat{S}$. As the stereo right to left
$4\times 4$ transformation matrix $T_{LR}$ can be easily obtained with tools
such as OpenCV [13], there are two transformations remaining to be estimated:
laser to stereo transformation matrix $T_{SL}$, and thermal to stereo
transformation matrix $T_{ST}$. We assume that the stereo and thermal cameras’
lens distortions have been calibrated such that the pinhole camera model is
applicable, and let $K^{left},K^{right},$ and $K^{thermal}$ be the intrinsic
matrices of these cameras, respectively.

Our calibration framework has three steps. First, it generates stereo point
clouds by image feature matching and triangulation. Then, we calibrate the
stereo and laser by point cloud registration. Finally, our model optimizes the
thermal-stereo transformation by minimizing the alignment error of the edge
features detected in the three sensors. The flow of the system is shown in
Fig. 2.

Figure 2: Flow chart of our calibration framework. It takes synchronized
stereo image pairs, thermal images, and laser point clouds as input, and
automatically calibrates the laser-stereo and thermal-stereo transformations
in one system.

###  II-A Stereo point cloud generation

We generate point clouds from stereo image pairs to make full use of their 3D
information in later calibrations. The Scale Invariant Feature Transform
(SIFT) feature detection algorithm [14] is employed to extract key points and
compute descriptors on $n$ stereo image pairs
$(I^{left}_{i},I^{right}_{i}),i=1\cdots n$. A RegionsMatcher implemented in
OpenMVG [15] is utilized to conduct the feature matching. Then the matched
features are triangulated to get the 3D positions $p^{stereo}_{i,j}$, where
$j$ is the index of matched features in image pair $i$. We use the OpenCV’s
[13] triangulation function here. The stereo point clouds
$C^{stereo}_{i}=\\{p^{stereo}_{i,j}\\},i=1\cdots n$. An example of stereo
point cloud generation is shown in Fig. 3.

Figure 3: Example of the matched feature points in a stereo image pair (left)
and the generated stereo point cloud (right).

###  II-B Stereo-laser calibration

The aim of this part is to estimate the laser to stereo transformation
$T_{SL}$. We use the Iterative Closest Point (ICP) algorithm [16] to register
the stereo and laser point clouds. We modify the original algorithm to a
multi-frame ICP (MFICP) that considers point clouds from all frames in the
dataset to optimize a unified extrinsic. It takes an initial guess of the
laser extrinsic and iteratively optimizes the initial guess to minimize the
distance between the closest points in each corresponding point cloud pair
$(C^{stereo}_{i},C^{laser}_{i}),i=1\cdots n$, where $n$ is the number of
frames in the dataset. Specifically, in each iteration, suppose
$p^{i}_{k_{1}}\in C^{stereo}_{i}$, let $q^{i}_{k_{2}}$ be the closest point to
$T^{-1}_{SL}p^{i}_{k_{1}}$ in $C^{laser}_{i}$, the registration cost is
formulated as

$Cost_{MFICP}=\sum_{i=1}^{n}\sum_{k_{1},k_{2}}\left\|T^{-1}_{SL}p^{i}_{k_{1}}-q^{i}_{k_{2}}\right\|^{2}_{2}$
(1)

As a standard nonlinear optimization problem, we use Ceres solver [17] to
optimize the initial guess of $T_{SL}$.

###  II-C Thermal extrinsic calibration

Visual-thermal and laser-thermal calibrations are challenging because of the
appearance difference and modality gap [1]. Therefore, finding common features
in RGB images, thermal images, and laser point clouds is the key to build
connections. We observed that although the interior regions of textured object
surfaces are not very distinguishable in thermal images, the object shape and
edges are usually clear. Besides, the object outlines are also distinguishable
in RGB and laser data. This leads us to the solution of calibrating the
thermal extrinsic by aligning the edges in RGB images, thermal images, and
laser point clouds.

####  II-C1 Stereo edge points detection

Detecting stereo edge points is straightforward since we can first detect
edges on stereo images and then triangulate the matched 2D edge points to 3D
space (see Section II-A). We use Sobel operator [18] to detect edges on each
stereo image pair $(I^{left}_{i},I^{right}_{i}),i=1\cdots n$, then pick out
the image feature points on these edges, and mark their corresponding 3D
points as stereo edge points. Fig. 4 gives one example of the stereo edge
point detection. The edge points in the $i^{th}$ stereo cloud forms
$E^{stereo}_{i}$.

Figure 4: Example of the matched features on edges in a stereo image pair
(left) and the edge points (marked in red) detected in the stereo point cloud
(right).

####  II-C2 Laser edge points detection

Because the laser point density along the perpendicular direction of the
scanning lines is too small, we only use the horizontal depth difference
between adjacent scans to detect laser edges. For a laser point $p$, we
collect its $k$ relative left points and $k$ relative right points, written as
$N_{0}$ and $N_{1}$, as shown in Fig. 5 (left). Then we define binary states
$\alpha_{\lambda}$ as whether all points in $N_{\lambda}$ is close to $p$, and
$\beta_{\lambda}$ as whether all points in $N_{\lambda}$ is deeper than $p$,

$\displaystyle\alpha_{\lambda}$ $\displaystyle=\forall p^{\prime}\in
N_{\lambda},\left|\|p^{\prime}\|_{2}-\|p\|_{2}\right|\leq\epsilon$ (2)
$\displaystyle\beta_{\lambda}$ $\displaystyle=\forall p^{\prime}\in
N_{\lambda},\|p^{\prime}\|_{2}-\|p\|_{2}>\epsilon,$ $\displaystyle\lambda=0,1$

where $\epsilon$ is the depth difference threshold. Then the laser edge points
are picked by

$\displaystyle p\text{ is edge
point}\Leftarrow\left(\alpha_{0}\land\beta_{1}\right)\lor\left(\beta_{0}\land\alpha_{1}\right)$
(3)

This picks the points with smaller depth at laser depth discontinuity edges,
which will still be visible when viewed from the camera (while the deeper ones
may not), as shown in Fig. 5 (left). An example of the laser edge point
detection result is given in Fig. 5 (right). All edge points in the $i^{th}$
laser frame forms $E^{laser}_{i}$.

Figure 5: _Left:_ Schematic diagram of laser edge points detection with
sampling radius $k=3$. _Right:_ The detected edge points (marked in red) in a
laser point cloud.

####  II-C3 Thermal edge detection and attraction field map generation

We use the Canny edge detector [19] to detect edges on thermal images
$I^{thermal}_{i},i=1\cdots n$. The edges shorter than $50$ pixels and the
cluttered interior edges are filtered out to improve the edge quality. One
example of filtered thermal edge map is shown in Fig. 6 (middle).

Figure 6: Example of a thermal image (left), its edge map (middle), and its
attraction field map (right). Each pixel value in the attraction field map
represents the normalized L2 distance from itself to the nearest edge.

With the principle that there should be correspondences between stereo, laser,
and thermal edges, we project each stereo and laser edge point to the
synchronized thermal image and take the closest thermal edge as its
correspondence. However, since individual stereo/laser edge point to thermal
edge pixel correspondence is unstable and with many outliers, inspired by
Levinson and Thrun’s work [8], we build thermal edge attraction field maps to
find the correspondences globally. The attraction field maps, written as
$G_{i},i=1\cdots n$, are generated by applying distance transform on thermal
edge maps. According to the properties of distance transform, each pixel in
$G_{i}$ stores the L2 distance from itself to the nearest edge. Therefore, the
reversed gradient of attraction field map points to the nearest edge, which
implicitly indicates the edge correspondences and “attracts” the projected
point to achieve a better edge alignment. One example of attraction field map
is shown in Fig. 6 (right).

####  II-C4 Edge alignment optimization

In this step, we calibrate thermal extrinsic by minimizing the reprojection
error of aligned edges (REAE). We first define two projection functions, which
project points in stereo point clouds and laser frames onto thermal images,
respectively,

$\displaystyle proj^{stereo}(p)=\left[\begin{matrix}u\\\ v\\\
1\end{matrix}\right]_{p}=\frac{1}{Z_{p}}K^{thermal}\left(T^{-1}_{ST}\left[\begin{matrix}x_{p}\\\
y_{p}\\\ z_{p}\\\ 1\end{matrix}\right]\right)_{1:3}$ (4a) $\displaystyle
proj^{laser}(q)=\left[\begin{matrix}u\\\ v\\\
1\end{matrix}\right]_{q}=\frac{1}{Z_{q}}K^{thermal}\left(T^{-1}_{ST}\,T_{SL}\left[\begin{matrix}x_{q}\\\
y_{q}\\\ z_{q}\\\ 1\end{matrix}\right]\right)_{1:3}$ (4b)

where the laser extrinsic $T_{SL}$ has been estimated in laser calibration
(Section II-B); the thermal extrinsic $T_{ST}$ is the optimization target with
a rough initial guess. Then, we define the cost function REAE as

$\displaystyle REAE=$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\widehat{E}^{stereo}_{i}}G_{i}\left(proj^{stereo}(p)\right)$
(5) $\displaystyle+$
$\displaystyle\sum_{i=1}^{n}\sum_{q\in\widehat{E}^{laser}_{i}}G_{i}\left(proj^{laser}(q)\right)$

where $G_{i}$ is the attraction field map defined in Section II-C3,
$\widehat{E}^{stereo}_{i}$ and $\widehat{E}^{laser}_{i}$ are sets of inlier
edge points. The inliers are selected from the stereo and laser edge points
based on the distance of their projections to the nearest thermal edge.

$\displaystyle\widehat{E}^{stereo}_{i}$ $\displaystyle=\\{\;p\in
E^{stereo}_{i}\;|\;G_{i}\left(proj^{stereo}(p)\right)\leq th\;\\}$ (6a)
$\displaystyle\widehat{E}^{laser}_{i}$ $\displaystyle=\\{\;q\in
E^{laser}_{i}\;|\;G_{i}\left(proj^{laser}(q)\right)\leq th\;\\}$ (6b)

where $th$ is the distance threshold.

We can transform the optimization problem to an unconstrained form with Lie
algebra. For a stereo edge point $p\in\widehat{E}^{stereo}_{i}$, we define
$p^{\prime}=R_{ST}p+t_{ST}=\left[X\;Y\;Z\right]^{T}$, $u=proj^{stereo}(p)$,
the Lie algebra of the thermal extrinsic $T_{ST}$ is $\xi$ and its left
perturbation is $\delta\xi$. Then the Jacobian matrix of the REAE cost is

$\frac{\partial REAE}{\partial T_{ST}}=\frac{\partial G_{i}}{\partial
u}\frac{\partial u}{\partial p^{\prime}}\frac{\partial
p^{\prime}}{\partial\delta\xi}\delta\xi$ (7)

where $\frac{\partial G_{i}}{\partial u}$ is the image gradient of $G_{i}$ at
$u$, and the rest part is the standard 3-D to 2-D projection model [1],

$\small{\frac{\partial u}{\partial p^{\prime}}\frac{\partial
p^{\prime}}{\partial\delta\xi}=\left[\begin{matrix}\frac{f_{x}}{Z}&0&-\frac{f_{x}X}{Z^{2}}&-\frac{f_{x}XY}{Z^{2}}&f_{x}+\frac{f_{x}X^{2}}{Z^{2}}&-\frac{f_{x}Y}{Z}\\\
0&\frac{f_{y}}{Z}&-\frac{f_{y}Y}{Z^{2}}&-f_{y}-\frac{f_{y}Y^{2}}{Z^{2}}&\frac{f_{y}XY}{Z^{2}}&\frac{f_{y}X}{Z}\end{matrix}\right]}$
(8)

The Jacobian matrix for laser edge point $q\in\widehat{E}^{laser}_{i}$ is
similar as Eq. 7.

Although the cost function is not globally convex, in reality it is always
locally convex around the correct calibration [8], thus we can use Ceres
solver [17] to optimize it with a rough initialization.

####  II-C5 Rough calibration

The proposed optimization-based method is designed for achieving high accuracy
calibration with a rough initialization. To further improve the initialization
error tolerance capability, we introduced a rough calibration procedure before
the optimization. It consists of two consecutive grid searches on rotation
(grid size $None$) and translation (grid size $4cm$) respectively, over a
given range around the initial guess of the thermal extrinsic, and pick the
transformation which maximizes the number of inlier points,
$|\widehat{E}^{stereo}_{i}|+|\widehat{E}^{laser}_{i}|$. The rough calibration
result is then delivered to the optimizer as the initial value. With the rough
calibration, our framework can work properly under large initialization
errors. The user only needs to provide a rough initial guess of the thermal
extrinsic, which can be obtained with ruler and protractor, or even by visual
estimation.

##  III Experimental Results

In this section, we present the evaluation of our calibration framework on
real-world datasets. We first introduce our sensor setups and collected
datasets, then evaluate the performance of the laser extrinsic calibration and
thermal extrinsic calibration, respectively.

Figure 7: _Left:_ The entire UAV platform. The sensors are mounted under the
drone. _Right:_ The installation position and orientation of the stereo
cameras, thermal camera, and laser sensors. Notice that there are two laser
sensors on our platform, but we only use the horizontally installed one.
Figure 8: Stereo and laser edge projection (red dots) on thermal images before
and after thermal extrinsic calibration on the 4 datasets. _First row:_
Projection with rough initial extrinsic; _Second row:_ Projection with our
calibrated extrinsic.

###  III-A Sensor setups and data collection

The sensors we use include two Ximea MC124CG-SY RGB cameras, one FLIR Boson®
640 Longwave Infrared (LWIR) thermal camera, and one Ouster OS0-128 laser
sensor. The two RGB cameras are fixed on both sides of the laser sensor to
form a stereo pair with $baseline=22.270cm$; the thermal camera is installed
between the two RGB cameras, next to the left camera; and the laser sensor is
installed a little higher in the middle. All the sensors are hardware-
synchronized by a PPS trigger. The stereo cameras, laser, and thermal camera
are running at 2Hz, 10Hz, and 50Hz, respectively. The sensor group is mounted
under a UAV platform to move freely in space. Fig. 7 shows the overall
appearance of the UAV platform and the sensor installation schematic diagram.
The stereo and thermal cameras’ intrinsics $K^{left},K^{right},K^{thermal}$,
and their distortion models have been calibrated beforehand.

We collect 4 sets of data on the CMU campus. Each data set contains
synchronized stereo RGB, thermal, and laser data. Fig. 8 shows the 4 scenes
with one thermal frame, which contain both certain objects and background
structures.

###  III-B Laser extrinsic calibration results

In this section, we test the performance of our framework in calibrating the
laser extrinsic on the 4 datasets. To give a quantitative comparison, we
manually calibrate a laser extrinsic ground truth. It is obtained by a PnP
solver on manually matched laser points and RGB image pixels over 10 different
frames.

To validate the robustness of our laser calibration algorithm under different
initialization values, we randomly sample 20 initial guesses around the ground
truth (offset by $$$-$$$ for rotation, $16-24cm$ for translation). We conduct
our laser calibration algorithm with each of the 20 initialization values on
the 4 datasets respectively, totaling 80 tests. The extrinsic errors before
and after calibration are shown in Fig. 9. It can be observed that both
rotation and translation errors decrease significantly after the calibration,
and converge to almost the same value in each scene regardless of the initial
value distribution. The final rotation errors are less than $None$ and the
translation errors are less than $5cm$ among all datasets. Fig. 10 shows the
projection of laser near points ($depth<10m$) on RGB images with initial and
calibrated extrinsic parameters. By observing the alignment of the projected
points and the outline of the stone pier in the image, it can be seen that the
calibration is accurate. These results suggest that our laser extrinsic
calibration algorithm has good performance even with large initialization
errors.

Figure 9: Laser extrinsic rotation and translation errors of the rough initial
extrinsic and our calibrated extrinsic in all scene settings. Figure 10: Laser
near points projection (red dots) on RGB image with initial (left) and
calibrated laser extrinsic (right).

###  III-C Thermal extrinsic calibration result

In this section, we validate our thermal extrinsic calibration algorithm in
terms of cost convergence, calibration accuracy, and comparison with the
checkerboard-based method. The thermal extrinsic ground truth is obtained by
solving a PnP problem on the manually picked triangulated stereo points
(actually from stereo images) and corresponding thermal pixels. Like the laser
calibration validation, we randomly bias the ground truth to get 20 extrinsic
initialization values for each dataset with rotation error range at $$$-$$$
and translation error range at $8-12cm$. Thus totally we have 80 tests with
different rough extrinsic initialization values.

####  III-C1 Cost convergence validation

To verify the convergence of the cost function, we draw the cost values over
the number of iterations in the 4 datasets with variant initialization values,
as shown in Fig. 11. As can be seen, the cost values gradually decrease with
the number of iterations, and the variances also become smaller. All tests are
sufficiently converged after 20 iterations. This demonstrates the
effectiveness and robustness of our thermal calibration algorithm.

Figure 11: REAE costs over the number of iterations in all scene settings. The
solid lines and dots are the average costs of all 20 tests in each scene while
the shaded areas indicate the standard deviation.

####  III-C2 Calibration accuracy validation

To quantitatively evaluate the accuracy, we calculate the extrinsic errors of
the 80 tests relative to the ground truth, as shown in Fig. 12. It is clear
that both the rotation and translation errors are decreased after calibration,
and the calibration results converge to almost the same value in each scene,
regardless of the large variance of the initialization values. The final
rotation errors are less than $None$ and the translation errors are less than
$4cm$ among all tests. Edge points projection with the initial and calibrated
extrinsic are shown in Fig. 8. It is clear that the stereo, laser, and thermal
edges align well after the calibration.

Figure 12: Thermal extrinsic rotation and translation error of the rough
initial extrinsic and our calibrated extrinsic in all scene settings.

####  III-C3 Comparisons

To further demonstrate the effectiveness of our proposed method, we compare it
with a traditional checkerboard-based approach. Follow the method described in
[12], we make a special checkerboard that is distinguishable in both RGB and
thermal images. We use OpenCV’s [13] checkerboard detector to detect
checkerboard corner points in stereo and thermal images, match them according
to their geometric positions, and solve a PnP problem on triangulated stereo
corners and corresponding thermal corners. We also test our method with a
visual estimation initialization on this dataset. Table. I lists the
comparison results of the calibration errors. It can be seen that our method
achieves similar accuracy to the checkerboard-based one without using a
checkerboard, which demonstrates the flexibility of our proposed method.

TABLE I: Thermal extrinsic calibration comparison | Our Method | Checkerboard [12]
---|---|---
Rot. Err. (∘) | 0.15 | 0.17
Trans. Err. (cm) | 3.24 | 2.30

##  IV CONCLUSIONS

This paper presents a multi-modal extrinsic calibration framework among stereo
cameras, thermal cameras, and laser sensors. We first use multi-frame ICP to
calibrate the extrinsic between stereo cameras and laser sensors, then project
registered stereo and laser edge points onto thermal images and optimize the
edge alignment to estimate the thermal extrinsic. Our framework can work
automatically in arbitrary environments without relying on specific
calibration targets. Since the system is robust to large initialization
errors, the initial guess can be simply obtained by rulers, protractor
measurement or visual estimation, which makes our framework ease to use. The
system greatly reduces the complexity of the calibration process while
achieving the same level of calibration accuracy as traditional methods.

Future works will focus on exploiting the real-time performance potential of
the system by pre-filtering irrelevant information and improving the
parallelism capability. The real-time calibrator can be used to continuously
detect and correct the extrinsic parameters’ changes due to the sensor drift
or time misalignment. Besides, more kinds of edge information, such as the
edges in laser reflectivity maps, can be integrated while optimizing cloud-
image edge alignment to improve the system’s accuracy and robustness.

##  V Acknowledgement

The authors acknowledge the sponsorship of this work from the Shimizu
Institute of Technology (Tokyo).

## References

  * [1] H. Yu, W. Zhen, W. Yang, and S. Scherer, “Line-based 2d-3d registration and camera localization in structured environments,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 69, no. 11, pp. 8962–8972, 2020.
  * [2] X. Zhang, P. Ye, D. Qiao, J. Zhao, S. Peng, and G. Xiao, “Object fusion tracking based on visible and infrared images using fully convolutional siamese networks,” in _International Conference on Information Fusion (FUSION)_ , 2019, pp. 1–8.
  * [3] Y. Cao, B. Xu, Z. Ye, J. Yang, Y. Cao, C.-L. Tisse, and X. Li, “Depth and thermal sensor fusion to enhance 3d thermographic reconstruction,” _Optics express_ , vol. 26, no. 7, pp. 8179–8193, 2018.
  * [4] A. Beauvisage, K. Ahiska, and N. Aouf, “Multimodal tracking framework for visual odometry in challenging illumination conditions,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 11 133–11 139.
  * [5] D. Scaramuzza, A. Harati, and R. Siegwart, “Extrinsic self calibration of a camera and a 3d laser range finder from natural scenes,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2007, pp. 4164–4169.
  * [6] P. Núñez, P. Drews Jr, R. P. Rocha, and J. Dias, “Data fusion calibration for a 3d laser range finder and a camera using inertial data,” in _European Conference on Mobile Robots (ECMR)_ , 2009, pp. 31–36.
  * [7] G. Pandey, J. R. McBride, S. Savarese, and R. M. Eustice, “Automatic targetless extrinsic calibration of a 3d lidar and camera by maximizing mutual information,” in _AAAI Conference on Artificial Intelligence_ , 2012, pp. 2053–2059.
  * [8] J. Levinson and S. Thrun, “Automatic online calibration of cameras and lasers,” in _Robotics: Science and Systems_ , 2013, pp. 29–36.
  * [9] C. Guindel, J. Beltrán, D. Martín, and F. García, “Automatic extrinsic calibration for lidar-stereo vehicle sensor setups,” in _IEEE international conference on intelligent transportation systems (ITSC)_ , 2017, pp. 674–679.
  * [10] A. Dhall, K. Chelani, V. Radhakrishnan, and K. M. Krishna, “Lidar-camera calibration using 3d-3d point correspondences,” _arXiv preprint arXiv:1705.09785_ , 2017.
  * [11] Y. Li, J. Tan, Y. Zhang, W. Liang, and H. He, “Spatial calibration for thermal-rgb cameras and inertial sensor system,” in _IEEE International Conference on Pattern Recognition (ICPR)_ , 2018, pp. 2295–2300.
  * [12] S. S. Shivakumar, N. Rodrigues, A. Zhou, I. D. Miller, V. Kumar, and C. J. Taylor, “Pst900: Rgb-thermal calibration, dataset and segmentation network,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 9441–9447.
  * [13] G. Bradski, “The opencv library,” _Dr. Dobb’s Journal of Software Tools_ , 2000.
  * [14] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International Journal of Computer Vision_ , vol. 60, no. 2, pp. 91–110, 2004\.
  * [15] P. Moulon, P. Monasse, R. Perrot, and R. Marlet, “Openmvg: Open multiple view geometry,” in _International Workshop on Reproducible Research in Pattern Recognition_ , 2016, pp. 60–74.
  * [16] B. PaulJ and M. NeilD, “A method for registration of 3-d shapes,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 14, no. 2, pp. 239–256, 1992.
  * [17] S. Agarwal, K. Mierle, and Others, “Ceres solver,” http://ceres-solver.org.
  * [18] N. Kanopoulos, N. Vasanthavada, and R. L. Baker, “Design of an image edge detection filter using the sobel operator,” _IEEE Journal of Solid-state Circuits_ , vol. 23, no. 2, pp. 358–367, 1988.
  * [19] J. Canny, “A computational approach to edge detection,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , no. 6, pp. 679–698, 1986.