# Simulations of quantum dynamics with fermionic phase-space representations
using numerical matrix factorizations as stochastic gauges

François Rousse1, Massimiliano Fasi1,2, Andrii Dmytryshyn1, Mårten
Gulliksson1, Magnus Ögren1,3 1School of Science and Technology, Örebro
University, 70182 Örebro, Sweden
2Department of Computer Science, Durham University, Stockton Road, DH1 3LE, UK
3Hellenic Mediterranean University, P.O. Box 1939, GR-71004, Heraklion, Greece

###### Abstract

The Gaussian phase-space representation can be used to implement quantum
dynamics for fermionic particles numerically. To improve numerical results, we
explore the use of dynamical diffusion gauges in such implementations. This is
achieved by benchmarking quantum dynamics of few-body systems against
independent exact solutions. A diffusion gauge is implemented here as a so-
called noise-matrix, which satisfies a matrix equation defined by the
corresponding Fokker–Planck equation of the phase-space representation. For
the physical systems with fermionic particles considered here, the numerical
evaluation of the new diffusion gauges allows us to double the practical
simulation time, compared with hitherto known analytic noise-matrices. This
development may have far reaching consequences for future quantum dynamical
simulations of many-body systems.

##  I Introduction

Quantum dynamics of few-body systems are benchmarked in order to study a new
numerical method for diffusion gauges in the fermionic Gaussian Phase-Space
Representation (GPSR) Corney-fermionicII ; Corney-fermionicI  with the goal to
extend the practical usefulness of the method.

A diffusion gauge is manifested by a time-dependent matrix
$B\in\mathbb{C}^{n\times m}$ that satisfies the matrix equation

$D=BB^{T},$ (1)

where $D\in\mathbb{C}^{n\times n}$ represents the coefficients of second-order
derivatives in the corresponding Fokker–Planck Equation (FPE). By using the
operator mappings of the positive-$P$ Drummond_JoPA_1980  for bosons and the
fermionic Gaussian phase-space representation Corney-fermionicII ; Corney-
fermionicI , one can derive the FPE

$\frac{\partial P\left(\vec{z}\right)}{\partial
t}=-\sum_{j}\frac{\partial\left\\{a_{j}\left(\vec{z}\right)P\left(\vec{z}\right)\right\\}}{\partial
z_{j}}+\frac{1}{2}\sum_{j,k}\frac{\partial^{2}\left\\{D_{j,k}\left(\vec{z}\right)P\left(\vec{z}\right)\right\\}}{\partial
z_{j}\partial z_{k}},\ \vec{z}\in\mathbb{C}^{n},$ (2)

for the time-dependent probability distribution $P\left(\vec{z}\right)$ of
bosonic-, fermionic-, or mixed-bosonic-fermionic- systems containing terms
with up to four creation- and annihilation-operators in a second quantized
Hamiltonian.

We can rewrite the FPE (2) as the Stochastic Differential Equation (SDE) in
the Itô interpretation GardinerBook1

$\dot{\vec{z}}=\vec{a}\left(\vec{z}\right)+\vec{b}(\vec{z},\vec{\eta}),\ \ \
\vec{b}(\vec{z},\vec{\eta})=B\left(\vec{z}\right)\vec{\eta},$ (3)

where $\vec{\eta}$ is an $m\times 1$ Gaussian noise vector, and the dot
denotes the derivative with respect to time, $d/dt$. We need to construct a
matrix $B\in\mathbb{C}^{n\times m}$ that satisfies the matrix equation (1)
GardinerBook1 , where $D$ is the matrix in (2). The freedom in $B$ defined by
(1) is known as a diffusion gauge Plimak_PRA_2001 ; Deuar_PRA_2002 . This work
focuses on the use of diffusion gauges that improves the durability in real
time of the numerical implementation of GPSR. Note that, in principle, the
matrix $B$ can have any number of columns, and that a larger number of columns
will increase the number of noise components in $\vec{\eta}$.

The stochastic averages over the complex phase-space variables $z_{j}$ in (2)
and (3) are directly related to first-order physical quantum-operators,
$\hat{O}$, moments, i.e., $\langle
z_{j}\rangle\leftrightarrow\bigl{\langle}\hat{O}_{j}\bigr{\rangle}$, and
indirectly related to higher-order physical operator moments Corney-fermionicI
.

In the absence of boundary corrections Corney-fermionicI , and impractically
large sampling errors, the stochastic averages of the phase-space variables
will approach the exact quantum mechanical expectation values as we take more
and more trajectories, each trajectory being a particular realization of the
SDE (3). After a certain simulation time, however, boundary corrections may
appear, or individual trajectories may tend to infinity, causing the numerical
result to become unreliable from that point onward. We refer to the time
elapsed before the occurrence of these issues as the practical simulation
time.

When using imaginary time to represent the inverse temperature $\sim T^{-1}$,
similar boundary corrections hamper the calculations of fermionic
groundstates, i.e., the $T\rightarrow 0$ limit. This phenomenon was first
documented for bosonic systems Smith_PRA_1989 ; Gilchrist_PRA_1997 , and was
then also observed by researchers examining fermionic systems Assaad_PRB_2005
; Corboz_PRB_2008 , who managed to alleviate it by means of projection methods
Assaad_PRB_2005 ; Corboz_PRB_2008 ; Aimi_and_Imada_PRA_2007 .

However, for quantum dynamics Ogren_EPL_2010 ; Ogren_CPC_2011 ;
Corboz_bookchapter_2013 , i.e., real time evolution, boundary corrections have
always been observed to be accompanied by clear signals, known as spiking
trajectories.

We propose and study a new numerical method for diffusion gauges in the
fermionic Gaussian phase-space representation Corney-fermionicII ; Corney-
fermionicI . We assess the new technique by benchmarking quantum dynamics of
few-body systems against known exact solutions.

We also illustrate such spiking trajectories in calculations of fermionic
quantum dynamics. In particular, we show by numerical examples that the use of
numerical diffusion gauges can delay the onset of spiking trajectories, and
hence prolong the practical simulation time. A numerical diffusion gauge is
manifested by a matrix $B$, that for each discrete timestep in the
implementation of the SDE, is a numerical solution to the matrix equation (1).

Numerical diffusion gauges exhibit a second advantage: in practice, it may be
cumbersome to find an analytic solution $B$ to (1) for a general Hamiltonian,
although various specific examples are known Ogren_EPL_2010 ; Ogren_CPC_2011 ;
Corboz_bookchapter_2013 .

In the following sections, we discuss three examples of phase-space
representations. First, we introduce the well-known anharmonic oscillator for
bosons Quantum_Optics_Book , for which, as mentioned, the role of diffusion
gauges has been documented in the literature for at least 20 years. Then, we
describe a system governed by the Fermi–Hubbard Hamiltonian. Finally, we
discuss a mixed system that exemplifies the advantages of applying phase-space
methods to a fermionic-bosonic Hamiltonian. This is the so-called Fermi–Bose
model Ogren_JPA_2013 , which is here used to model the conversion from a
molecular BEC to pairs of fermionic atoms.

These three examples of a bosonic-, fermionic-, and a mixed-system are of
interest for example in quantum optics, solid state physics, and quantum atom-
optics, respectively. The method we consider, however, is more general and may
find use for the exact quantum dynamics of many other Hamiltonians in
different fields of applications.

Figure 1: (Colour online) Positive-$P$ dynamics for $20$ individual
trajectories with different diffusion gauges. Left: Standard Positive-$P$,
with the noise matrix $B$ of equation (23). Right: “optimal” noise matrix, $B$
of equation (24) with $\gamma\simeq 3.1985$, see Deuar_PRA_2002 . Note the
different scales on the two logarithmic y-axes. Parameters and initial
conditions are $\omega=0,\>N=10^{4},\>\chi=10^{-2}$ (i.e.,
$t_{\textnormal{opt}}=1/(\sqrt{N}\chi)=1$), respectively, according to
equation (6), as in Plimak_PRA_2001 .  Figure 2: (Colour online) Benchmarking
of the time-dependent correlation function
$|\langle\hat{a}^{\dagger}(t)\hat{a}(0)\rangle|/N=|\langle
z_{2}(t)\rangle|/\sqrt{N}$ as a function of dimensionless time Plimak_PRA_2001
. The dashed (black) curve describe the analytic result of equation (14) in
Plimak_PRA_2001 . Solid (blue) curve is numerical results from the phase-space
method for $t\leq 0.5$ which is close to the largest useful simulation time
with the specific implementation of (5), while the solid (green) curve for
$t\leq 3$ shows gauge-improved numerical results using instead a $B$ matrix
from equation (24). We have sampled $10^{4}$, green solid, and $10^{2}$, cyan
thin, to show the dependence on the number of trajectories. The dashed-dotted
(red) horizontal line describe the corresponding Gross-Pitaevskii mean-field
result, that do not vary with time in this case. Parameters and initial
conditions are the same as in figure 1.

##  II Examples of phase-space representations

The first example, A, which is meant as a simple introduction to diffusion
gauges, comes from early literature on the topic. The following two examples,
i.e., B and C, are representations for fermionic systems—we will use these to
evaluate the performance of our numerical diffusion gauges in section III.

###  II.1 A bosonic positive-$P$ phase-space representation, for an
interacting bosonic quantum field

As an introduction to the use of SDEs for exact quantum dynamics, we begin by
revising the well-known time evolution of the one-mode Bose–Hubbard-like
Hamiltonian (the Kerr oscillator Quantum_Optics_Book )

$\widehat{H}=\omega\hat{a}^{\dagger}\hat{a}+\frac{\chi}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}.$
(4)

Note that here and in the remainder of this work we choose our units so that
$\hbar=1$.

We formulate the equations for realizations with two $B$ matrices that are
known analytically. If we combine equations (20) and (23) in Appendix A, we
obtain a specific system of SDEs
$\dot{\vec{z}}=\vec{a}(\vec{z})+\vec{b}(\vec{z},\vec{\eta})$
Drummond_JoPA_1980 , which we can write elementwise as

$\left[\begin{array}[]{c}\dot{z}_{1}\\\
\dot{z}_{2}\end{array}\right]=-i\left[\begin{array}[]{c}\omega z_{1}+\chi
z_{2}z_{1}^{2}\\\ -\omega z_{2}-\chi z_{1}z_{2}^{2}\end{array}\right]\\\
+\sqrt{i\chi}\left[\begin{array}[]{c}-iz_{1}\eta_{1}\\\ z_{2}\eta_{2}\\\
\end{array}\right],$ (5)

where $\vec{\eta}^{T}=\left[\eta_{1},\eta_{2}\right]$ denotes uncorrelated
real Gaussian noises with zero mean and unit variance. We stress that in (5)
we follow the notation of Gilchrist_PRA_1997  and Steel_PRA_1998 , but the
complex variables are defined differently by other authors Plimak_PRA_2001 .
The specific noise terms in (5) may follow from the one-parameter analytic
diffusion gauge given in (24) of Appendix A.

If the noise terms in (5) are set to zero ($\eta_{j}\equiv 0$), then
$z_{2}=z_{1}^{*}$ and we obtain the deterministic equation
$i\dot{z}_{1}=\left(\omega+\chi|z_{1}|^{2}\right)z_{1}$, which is equivalent
to the time-dependent mean-field Gross–Pitaevskii formalism BEC_book .

For the initial condition

$\vec{z}(0)^{T}=\sqrt{N}[1,1],$ (6)

where $N$ is the number of bosons, we calculate the dynamics of (5) for
different stochastic trajectories to obtain the average $\langle
z_{2}(t)\rangle$.

In figure 1, we illustrate the different behavior of the dynamics for two
choices of diffusion gauges: the standard positive-$P$ (left panel) and an
improved (“optimal”) form Plimak_PRA_2001  (right panel). As is clearly seen
from the 20 stochastic trajectories we report, which are realizations of the
system of SDEs (5), the variance of the standard positive-$P$ formulation
(left figure 1), is several orders of magnitude larger than that of the
“optimal” form Plimak_PRA_2001  (right figure 1). This different behavior will
heavily influence the practical usefulness of averages obtained from such
different trajectories.

We remark that the “optimal” form of the $B$-matrix used in figure 1 was found
by relying upon the existence of analytic solutions to certain correlations in
the problem Plimak_PRA_2001 . Plimak, Olsen and Collett Plimak_PRA_2001
demonstrated that changing the $B$-matrix so to dramatically reduce noise was
enough to change the method from “computable in principle” to “computable in
practice”. Although dependent on an analytic solution, this early example
motivated an investigation of the role of the $B$-matrix in more difficult
problems, for which analytic solutions are not known—a setting in which
stochastic simulations can clearly be useful.

The stochastic average $\langle z_{2}(t)\rangle$ is here related to the
quantum-operator average $\langle\hat{a}^{\dagger}(t)\hat{a}(0)\rangle$, given
analytically in equation (14) of Plimak_PRA_2001 . In figure 2, we use the
analytic formula there to numerically illustrate the role of a diffusion gauge
in this initial bosonic example.

###  II.2 A fermionic Gaussian phase-space representation, for the
Fermi–Hubbard model

We study the time evolution of a Fermi–Hubbard Hamiltonian with nearest-
neighbour ($\left\langle i,j\right\rangle$) jumps and on-site interaction,

$\widehat{H}=-J\sum\nolimits_{{\left\langle
i,j\right\rangle},\sigma}\hat{c}_{{i},\sigma}^{\dagger}\hat{c}_{{j},\sigma}+U\sum\nolimits_{{j}}\hat{n}_{{j},\uparrow}\hat{n}_{{j},\downarrow},$
(7)

where $J$ is the hopping amplitude and $U$ the interaction strength. We focus
on a small two-site ($j=1$, $2$) chain in order to obtain detailed numerical
comparisons for fermionic quantum dynamics.

First, we formulate the equations for realizations of the GPSR for the system
in Corboz_bookchapter_2013  with different choices of diffusion gauges. We
compare the results with independent numerical solutions in the number state
representation (so-called exact diagonalization).

Combining equations (25) and (28) in Appendix A gives the following set of
SDEs, $\dot{\vec{z}}=\vec{a}(\vec{z})+\vec{b}(\vec{z},\vec{\eta})$
Corboz_bookchapter_2013 , or equivalently,

$\left[\begin{array}[]{c}\dot{z}_{1}\\\ \dot{z}_{2}\\\ \dot{z}_{3}\\\
\dot{z}_{4}\\\ \dot{z}_{5}\\\ \dot{z}_{6}\\\ \dot{z}_{7}\\\
\dot{z}_{8}\end{array}\right]=i\left[\begin{array}[]{c}J\left(z_{2}-z_{3}\right)\\\
J\left(z_{1}-z_{4}\right)+U\left(z_{5}-z_{8}\right)z_{2}\\\
J\left(z_{4}-z_{1}\right)+U\left(z_{8}-z_{5}\right)z_{3}\\\
J\left(z_{3}-z_{2}\right)\\\ J\left(z_{6}-z_{7}\right)\\\
J\left(z_{5}-z_{8}\right)+U\left(z_{1}-z_{4}\right)z_{6}\\\
J\left(z_{8}-z_{5}\right)+U\left(z_{4}-z_{1}\right)z_{7}\\\
J\left(z_{7}-z_{6}\right)\end{array}\right]+\sqrt{iU}\left[\begin{array}[]{c}i(\tilde{z}_{1}z_{1}\xi_{1}-z_{2}z_{3}\xi_{2})+\tilde{z}_{1}z_{1}\xi_{3}-z_{2}z_{3}\xi_{4}\\\
i(\tilde{z}_{4}z_{2}\xi_{2}-z_{1}z_{2}\xi_{1})+\tilde{z}_{1}z_{2}\xi_{3}-z_{2}z_{4}\xi_{4}\\\
i(\tilde{z}_{1}z_{3}\xi_{1}-z_{3}z_{4}\xi_{2})+\tilde{z}_{4}z_{3}\xi_{4}-z_{1}z_{3}\xi_{3}\\\
i(\tilde{z}_{4}z_{4}\xi_{2}-z_{2}z_{3}\xi_{1})+\tilde{z}_{4}z_{4}\xi_{4}-z_{2}z_{3}\xi_{3}\\\
i(\tilde{z}_{5}z_{5}\xi_{1}^{*}-z_{6}z_{7}\xi_{2}^{*})+\tilde{z}_{5}z_{5}\xi_{3}^{*}-z_{6}z_{7}\xi_{4}^{*}\\\
i(\tilde{z}_{8}z_{6}\xi_{2}^{*}-z_{5}z_{6}\xi_{1}^{*})+\tilde{z}_{5}z_{6}\xi_{3}^{*}-z_{6}z_{8}\xi_{4}^{*}\\\
i(\tilde{z}_{5}z_{7}\xi_{1}^{*}-z_{7}z_{8}\xi_{2}^{*})+\tilde{z}_{8}z_{7}\xi_{4}^{*}-z_{5}z_{7}\xi_{3}^{*}\\\
i(\tilde{z}_{8}z_{8}\xi_{2}^{*}-z_{6}z_{7}\xi_{1}^{*})+\tilde{z}_{8}z_{8}\xi_{4}^{*}-z_{6}z_{7}\xi_{3}^{*}\\\
\end{array}\right],$ (8)

where the notation $\tilde{z}_{j}\equiv 1-z_{j}$ is used for the so-called
hole-variables. In (8), we use the complex noise

$\xi_{j}=\dfrac{\eta_{1}^{\left(j\right)}+i\eta_{2}^{\left(j\right)}}{\sqrt{2}},\quad
j=1,2,3,4,$ (9)

where $\eta_{1}^{\left(j\right)}$ and $\eta_{2}^{\left(j\right)}$ denote real
Gaussian noises with zero mean and unit variance. The complex Gaussian noises
$\xi_{j}$ obey the correlations

$\langle\xi_{j}(t)\xi_{j^{\prime}}(t)\rangle=0,\qquad\langle\xi_{j}(t)\xi_{j^{\prime}}^{*}\
(t^{\prime})\rangle=\delta_{j,j^{\prime}}\delta_{t,t^{\prime}}.$

If the noise-terms in (8) are all set to zero ($\xi_{j}\equiv 0$), then we
have a deterministic system that is equivalent to the time-dependent
Hartree–Fock formalism RahavPRB2009  for the Hamiltonian (7).

###  II.3 A mixed Fermi–Bose phase-space representation, for dissociation of
molecules

As an example of a bosonic-fermionic system, we consider the dissociation of a
molecular BEC of dimers into pairs of fermionic atoms, with the Hamiltonian

$\widehat{H}=\hbar\sum\nolimits_{\mathbf{k},\sigma}\Delta_{\mathbf{k}}\hat{c}_{\mathbf{k},\sigma}^{\dagger}\hat{c}_{\mathbf{k},\sigma}-i\hbar\kappa\sum\nolimits_{\mathbf{k}}\left(\hat{a}^{\dagger}\hat{c}_{\mathbf{k},\uparrow}\hat{c}_{-\mathbf{k},\downarrow}-\hat{c}_{-\mathbf{k},\downarrow}^{\dagger}\hat{c}_{\mathbf{k},\uparrow}^{\dagger}\hat{a}\right),$
(10)

where $\Delta_{\mathbf{k}}=\Delta+\hbar{\mathbf{k}}^{2}/(2m)$ is a parameter
for the kinetic energy of the free atoms, translated with respect to the
energy for the bound molecular state Ogren_EPL_2010 , and $\kappa$ is the
strength of the atom-molecular coupling.

With only one momentum mode $\mathbf{k}=1$, and symmetric initial conditions
with respect to the spin-variable $\sigma$, we only need one phase-space
variable for the fermionic normal moment
$\langle\hat{c}_{{\mathbf{k}},\sigma}^{\dagger}\hat{c}_{{\mathbf{k}},\sigma}\rangle\leftrightarrow\langle
z_{1}\rangle$ and one for the anomalous moment
$\langle\hat{c}_{{\mathbf{k}},\uparrow}\hat{c}_{{-\mathbf{k}},\downarrow}\rangle\leftrightarrow\langle
z_{2}\rangle$, while the bosons are again mapped via the positive-$P$
representation (see section II.1).

Here, we write down the minimal system with the 5 phase-space variables
($z_{1}=n_{1}$, $z_{2}=m_{1}$, $z_{3}=m^{+}_{1}$, $z_{4}=\alpha$,
$z_{5}=\alpha^{+}$, compare to Ogren_EPL_2010 ) in explicit form. From
equations (29) and (32) in Appendix A, we have the system of SDEs
$\dot{\vec{z}}=\vec{a}(\vec{z})+\vec{b}(\vec{z},\vec{\eta})$ Ogren_EPL_2010 ,
where

$\left[\begin{array}[]{c}\dot{z}_{1}\\\ \dot{z}_{2}\\\ \dot{z}_{3}\\\
\dot{z}_{4}\\\
\dot{z}_{5}\end{array}\right]=\left[\begin{array}[]{c}\kappa\left(z_{3}z_{4}+z_{2}z_{5}\right)\\\
-2i\Delta_{1}z_{2}+\kappa z_{4}\left(1-2z_{1}\right)\\\
2i\Delta_{1}z_{3}+\kappa z_{5}\left(1-2z_{1}\right)\\\ -\kappa z_{2}\\\
-\kappa z_{3}\\\
\end{array}\right]+\sqrt{\kappa}\left[\begin{array}[]{c}z_{1}\left(z_{2}\xi_{1}^{*}+z_{3}\xi_{2}^{*}\right)\\\
z_{2}^{2}\xi_{1}^{*}-z_{1}^{2}\xi_{2}^{*}\\\
-z_{1}^{2}\xi_{1}^{*}+z_{3}^{2}\xi_{2}^{*}\\\ \xi_{1}\\\
\xi_{2}\end{array}\right].$ (11)

If the noise-terms in (11) are all set to zero ($\xi_{j}\equiv 0$), then we
have a deterministic system that is equivalent to the time-dependent pairing-
meanfield formalism PMFT_Jack_and_Pu_PRA_2005  for the Hamiltonian (10).

##  III Numerical benchmarking of fermionic systems

Now we describe the numerical experiments used to benchmark the fermionic
systems presented in sections II.2 and II.3 above, we discuss the results we
obtained, and we show how the numerical diffusion gauge based GPSR method we
propose improves on existing alternatives found in the literature.

###  III.1 Results for the Fermi–Hubbard model

Figure 3: (Colour online) The occurrence of spikes in the stochastic dynamics
for the Fermi–Hubbard model for individual trajectories with different
diffusion gauges. Left: Standard GPSR from the literature
Corboz_bookchapter_2013 , i.e., with the analytic noise matrix $B$ of equation
(28). Right: GPSR with numerical noise matrices $B$, i.e., dynamically solving
equation (1). Thin (cyan) curves are individual trajectories. Thick (blue)
curves are the average of $10^{3}$ such trajectories, surrounded by $\pm$ the
stochastic errors. Dashed-dotted (red) curves are from the corresponding
Hartree–Fock meanfield-method RahavPRB2009 , while dashed (black) curves are
from numerical exact diagonalization, for comparisons. Parameters and initial
conditions are $J=U=1$, respectively, according to (12).

As initial condition $\vec{z}(0)$ for the Hamiltonian (7), we choose

$\vec{z}(0)^{T}=\left[1,0,0,0,1,0,0,0\right],$ (12)

which corresponds to the initial ($t=0$) quantum-operator averages

$\langle\hat{c}_{i,\sigma}^{\dagger}\hat{c}^{\phantom{\dagger}}_{j,\sigma}\rangle=\begin{cases}1,\qquad&i=j=1,\\\
0,&\text{otherwise}.\end{cases}$ (13)

We begin by examining the dynamics of the stochastic average $\langle
z_{1}\rangle$. In figure 3, we show $20$ individual trajectories and the
stochastic averages over $10^{3}$ trajectories for the standard GPSR with
analytic noise matrices (left panel) and the GPSR with numerical noise
matrices (right panel). The numerical diffusion gauges clearly outperform the
analytic ones, yielding a practical simulation time that is twice as long.
Importantly, the average of the stochastic method follows the curve of the
exact diagonalization, which differs from that of the Hartree–Fock meanfield-
method. On the contrary, we can hardly distinguish the three curves in the
left panel of figure 3.

We then investigate quantum correlations that map ($\leftrightarrow$) to
stochastic averages over combinations of phase-space variables for the
Fermi–Hubbard model. In figure 4, we show in the top panel the total number of
particles

$\langle\widehat{N}_{\textnormal{tot}}\rangle=\sum\nolimits_{j,\sigma}\langle\hat{n}_{{j},\sigma}\rangle\leftrightarrow\langle
z_{1}\rangle+\langle z_{4}\rangle+\langle z_{5}\rangle+\langle z_{8}\rangle,$
(14)

and in the bottom panel the total energy

$\langle\widehat{H}\rangle=-J\sum\nolimits_{{\left\langle
i,j\right\rangle},\sigma}\langle\hat{c}_{{i},\sigma}^{\dagger}\hat{c}_{{j},\sigma}\rangle+U\sum\nolimits_{{j}}\langle\hat{n}_{{j},\uparrow}\hat{n}_{{j},\downarrow}\rangle\leftrightarrow-J(\langle
z_{2}\rangle+\langle z_{3}\rangle+\langle z_{6}\rangle+\langle
z_{7}\rangle)+U(\langle z_{1}z_{5}\rangle+\langle z_{4}z_{8}\rangle).$ (15)

We only report the real part of (14) and (15), as these two constants of
motion are equal to the expectation values of the initial state (12), which
are $\langle\widehat{N}_{\textnormal{tot}}\rangle=2$ and
$\langle\widehat{H}\rangle=1$, respectively, and will thus have no imaginary
part.

The Hartree–Fock method will produce a horizontal line (not shown) for the
total energy, but this will be the sum of two terms that both deviate
substantially from the correct values of kinetic- and interaction-energies.

Figure 4: (Colour online) Benchmarking of time-dependent correlation functions
for the Fermi–Hubbard model with different diffusion gauges. Upper plot, total
number of particles according to (14). Lower plot, total energy according to
(15). Standard GPSR from the literature Corboz_bookchapter_2013 , as in the
left figure 3, (blue) curves that here spike at $t\simeq 0.7$. GPSR with
numerical noise matrices $B$, as in the right figure 3, (green) curves that
here spike at $t\simeq 1.5$. All four stochastic curves are averages of
$10^{3}$ trajectories. Dashed (black) curves are from exact diagonalization
and are horisontal as expected. Parameters and initial conditions are the same
as in figure 3.

###  III.2 Results for the Fermi–Bose model

Figure 5: (Colour online) The occurrence of spikes in the stochastic dynamics
for the Fermi–Bose model for individual trajectories with different diffusion
gauges. Left: Standard mixed phase-space representation from the literature
Ogren_EPL_2010 , i.e., with the analytic noise matrix $B$ of equation (32).
Right: GPSR with numerical noise matrices $B$, i.e., dynamically solving
equation (1). Thick (blue) curves is the average of $10^{3}$ such
trajectories, surrounded by $\pm$ the stochastic errors. Dashed-dotted (red)
curves are from the corresponding pairing-meanfield-method
PMFT_Jack_and_Pu_PRA_2005 , while dashed (black) curves are from numerical
exact diagonalization, for comparisons. Parameters and initial conditions are
$N_{m}(0)=\Delta_{1}=1$, and $\kappa=2$, respectively, according to (16).

As initial condition $\vec{z}(0)$ for the Hamiltonian (10), we choose

$\vec{z}(0)^{T}=\left[0,0,0,1,1\right],$ (16)

i.e., $\sqrt{N_{m}(0)}=1$, which corresponds to the initial ($t=0$) quantum-
operator averages

$\langle\hat{c}_{{1},\sigma}^{\dagger}\hat{c}^{\phantom{\dagger}}_{{1},\sigma}\rangle=\langle\hat{c}^{\phantom{\dagger}}_{{1},\uparrow}\hat{c}^{\phantom{\dagger}}_{{-1},\downarrow}\rangle=0,\
\textnormal{and}\ \langle\hat{a}^{\dagger}\hat{a}\rangle=N_{m}(0),$ (17)

where the latter is the initial ($t=0$) number of dimer-molecules.

In figure 5, we show $20$ individual trajectories and the stochastic average
$N_{m}=\langle\hat{a}^{\dagger}\hat{a}\rangle\leftrightarrow\langle
z_{4}z_{5}\rangle$ for the quantum dynamics of the number of molecules. In
this case, we compare the GPSR with numerical noise matrices (right panel)
with the standard mixed-phase space representation from the literature (left
panel), which uses an analytic noise matrix. Again, we observe a doubling of
the practical simulation time, and we remark that the average of the GPSR
method in the right panel of figure 5 follows the curve of the exact
diagonalization, which differs from that obtained by the pairing-meanfield
method. Once again, the three curves are hardly distinguishable within the
useful simulation time of the standard method, as is evidenced in the left
panel of figure 5.

Next, we investigate the total number of particles

$\langle\widehat{N}_{\textnormal{tot}}\rangle=\sum_{k}\langle\hat{n}_{k}\rangle+\langle\hat{a}_{0}^{\dagger}\hat{a}_{0}\rangle\leftrightarrow\langle
z_{1}\rangle+\langle z_{4}z_{5}\rangle,$ (18)

and the total energy

$\langle\widehat{H}\rangle=2\sum_{k}\Delta_{k}\langle\hat{n}_{k}\rangle+i\kappa\sum_{k}\langle\hat{a}_{0}\hat{m}_{k}^{\dagger}-\hat{a}_{0}^{\dagger}\hat{m}_{k}\rangle\leftrightarrow
2\Delta_{1}\langle z_{1}\rangle+i\kappa(\langle z_{4}z_{3}\rangle-\langle
z_{5}z_{2}\rangle).$ (19)

The constants of motion are equal to the expectation values of the initial
state (16), which are $\langle\widehat{N}_{\textnormal{tot}}\rangle=1$ (the
atom with negative momentum is not included explicitly in the model) and
$\langle\widehat{H}\rangle=0$ here. We report the observed values of the real
parts of these two quantities in figure 6, as they are both expected to have
no imaginary part.

For the total energy, the pairing-meanfield method will result in a horizontal
line (not plotted), but this will be the sum of two terms that both deviate
substantially from the correct values of the kinetic- and conversion-energy.

Figure 6: (Colour online) Benchmarking of time-dependent correlation functions
for the Fermi–Bose model with different diffusion gauges. Upper plot, total
number of particles according to (18). Lower plot, total energy according to
(19). Standard mixed phase-space representation from the literature
Ogren_EPL_2010 , as in the left figure 5, (blue) curves that here spike at
$t\simeq 0.5$. GPSR with numerical noise matrices $B$, as in the right figure
5, (green) curves that here spike at $t\simeq 1.2$. Averages were calculated
from $10^{3}$ trajectories, parameters, and initial conditions are the same as
in figure 5.

##  IV Summary and outlook

Phase-space representations for quantum dynamics are an established tool for
numerical simulations, but, often hampered by limited simulation times, they
can only describe the initial dynamics of some large bosonic Wuster_PRE_2017
and fermionic Ogren_EPL_2010  quantum systems. In order to unravel the exact
quantum dynamics for larger times, it is essential to gradually increase the
practical simulation time of these numerical computations.

Inspired by early work on analytic diffusion gauges for simple bosonic
systems, we introduced a novel method based on the numerical evaluation of
diffusion gauges, and we benchmarked the new technique against existing
alternatives using two different complex fermionic systems. For both the
Fermi–Hubbard model (figure 3) and the Fermi–Bose model (figure 5), we
observed a doubling of the practical simulation time.

These numerical results are very encouraging, and they could, for example,
trigger an investigation of the exact quantum dynamics of large 2D
Fermi–Hubbard models. Such models are currently beyond the reach of exact
diagonalization, and have proven hard for matrix-product-states-based methods.

As for the reasons behind the prolonged simulation times observed, we can only
speculate here. It has been noted in the literature Gilchrist_PRA_1997 ;
Corboz_PRB_2008  that the emergence of spiking trajectories in the phase-space
variables, i.e. a $|z_{j}|$ approaching infinite, can be connected to the
build-up of fat tails in the distribution for $P\left(\vec{z}\right)$, and
hence possible occurrence of boundary corrections. Since the phase-space basis
is overcomplete, alternative distributions $P\left(\vec{z}\right)$, and
alternative dynamics of individual trajectories, may still describe the
correct state. With the introduced numerical diffusion gauges, which are
different at each time step, the stochastic dynamics related to the
distribution $P\left(\vec{z}\right)$ can be different, and in particular less
likely to grow in the same directions in consecutive time steps.

In order to further increase the practical simulation time, research on a more
fundamental level may be required. The better we understand how these
numerical diffusion gauges work, the more precise the constrained matrix
equations we can construct to obtain them will be.

##  V Acknowledgements

We thank Joel Corney for useful discussions.

##  Appendix A Motivations of the SDEs for the phase-space representations

In this appendix, we present one explicit instance of a system of SDEs for
each of the three Hamiltonians in the main body of the article. In each of the
three examples, we keep the number of variables to a minimum.

###  A.1 The positive-$P$ representation

Now we explicitly present one possible system of SDEs for the Hamiltonian (4)
with one bosonic mode. Using the operator mapping for the bosonic positive-$P$
representation GardinerBook1 , we obtain the $a_{j}$ and $D_{j,k}$ in equation
(2). With $n=2$ phase-space variables ordered as $z_{1}=\alpha$,
$z_{2}=\alpha^{+}$, the terms in the FPE (2) with the first-order derivatives
are Steel_PRA_1998

$\left[\begin{array}[]{c}a_{1}\left(\vec{z}\right)\\\
a_{2}\left(\vec{z}\right)\end{array}\right]=-i\left[\begin{array}[]{c}\omega
z_{1}+\chi z_{2}z_{1}^{2}\\\ -\omega z_{2}-\chi z_{1}z_{2}^{2}\\\
\end{array}\right].$ (20)

The second-order derivatives in the FPE (2) are defined by the symmetric
$D\in\mathbb{C}^{2\times 2}$ matrix

$D=i\chi\left[\begin{array}[]{cc}d_{11}&0\\\ 0&d_{22}\end{array}\right],$ (21)

where the two nonzero elements are

$d_{11}=-z_{1}^{2},\quad d_{22}=z_{2}^{2}.$ (22)

The matrix $B\in\mathbb{C}^{2\times 2}$

$B=\sqrt{i\chi}\left[\begin{array}[]{cc}iz_{1}&0\\\
0&-z_{2}\end{array}\right],$ (23)

then satisfies (1) for $D$ in (21). Hence, equation (20) together with
equation (23) give the SDE (5). However, the following one-parameter analytic
$B$-matrix also satisfies (1)

$B(\gamma)=-i\sqrt{i\chi}\left[\begin{array}[]{rr}\cosh(\gamma)z_{1}&i\sinh(\gamma)z_{1}\\\
-\sinh(\gamma)z_{2}&-i\cosh(\gamma)z_{2}\end{array}\right].$ (24)

The $B(\gamma)$ in (24) is equivalent to the trigonometric form presented in
Deuar_PRA_2002  and was suggested for the positive-$P$ method for the Kerr
oscillator in Plimak_PRA_2001 . For suitable values of the parameter, the
noise matrix in (24) shows a remarkable improvement in obtaining observables,
as shown in figure 2, with lower variance for longer times, as shown in figure
1. This improvement in the mean is due to the fact that the amplitudes
$|z_{1}|$ and $|z_{2}|$ for the different trajectories, as seen in figures 2
and 3 in Plimak_PRA_2001 . For $\gamma=0$, equations (23) and (24) are equal,
while for $\gamma\simeq 3.1985$, corresponding to $\cosh(\gamma)=\sqrt{301/2}$
(see Appendix 2 of Deuar_PRA_2002 ) we get the improved result reported in
figures 1 and 2 in section I.

###  A.2 The two-site Fermi–Hubbard model

In this subsection we explicitly present one possible system of SDEs for the
Hamiltonian (7) with two spatial sites. Using the operator mapping of the
Gaussian phase-space representation Corney-fermionicII  with only normal
operator moments (no anomalous operator moments, i.e., no $m$ or $m^{+}$
variables are used here), one can derive a Fokker–Planck equation for the
time-dependent probability distribution $P\left(\vec{z}\right)$ of the form
(2).

Here we use the notation $z_{1}=n_{11,\uparrow}$, $z_{2}=n_{12,\uparrow}$,
$z_{3}=n_{21,\uparrow}$, $z_{4}=n_{22,\uparrow}$, $z_{5}=n_{11,\downarrow}$,
$z_{6}=n_{12,\downarrow}$, $z_{7}=n_{21,\downarrow}$,
$z_{8}=n_{22,\downarrow}$ for the eight complex phase-space variables we use
for the two-site Fermi–Hubbard model. The phase-space variables are related to
the first-order physical moments
$\left\langle\hat{n}_{ij,\sigma}\right\rangle\leftrightarrow\langle
n_{ij,\sigma}\rangle$, and their stochastic averages approach the quantum
mechanical expectation values of those in the limit of many trajectories.

The terms in the FPE (2) with the first-order derivatives are in this case
Corboz_bookchapter_2013

$\left[\begin{array}[]{c}a_{1}\left(\vec{z}\right)\\\
a_{2}\left(\vec{z}\right)\\\ a_{3}\left(\vec{z}\right)\\\
a_{4}\left(\vec{z}\right)\\\ a_{5}\left(\vec{z}\right)\\\
a_{6}\left(\vec{z}\right)\\\ a_{7}\left(\vec{z}\right)\\\
a_{8}\left(\vec{z}\right)\end{array}\right]=i\left[\begin{array}[]{c}J\left(z_{2}-z_{3}\right)\\\
J\left(z_{1}-z_{4}\right)+U\left(z_{5}-z_{8}\right)z_{2}\\\
J\left(z_{4}-z_{1}\right)+U\left(z_{8}-z_{5}\right)z_{3}\\\
J\left(z_{3}-z_{2}\right)\\\ J\left(z_{6}-z_{7}\right)\\\
J\left(z_{5}-z_{8}\right)+U\left(z_{1}-z_{4}\right)z_{6}\\\
J\left(z_{8}-z_{5}\right)+U\left(z_{4}-z_{1}\right)z_{7}\\\
J\left(z_{7}-z_{6}\right)\end{array}\right].$ (25)

The second-order differential operator in (2) are defined by the symmetric
$D\in\mathbb{C}^{8\times 8}$ matrix has the form

$D=\frac{iU}{2}\left[\begin{array}[]{cccccccc}0&0&0&0&0&d_{16}&d_{17}&0\\\
0&0&0&0&d_{25}&d_{26}&d_{27}&d_{28}\\\ 0&0&0&0&d_{35}&d_{36}&d_{37}&d_{38}\\\
0&0&0&0&0&d_{46}&d_{47}&0\\\ 0&d_{25}&d_{35}&0&0&0&0&0\\\
d_{16}&d_{26}&d_{36}&d_{46}&0&0&0&0\\\ d_{17}&d_{27}&d_{37}&d_{47}&0&0&0&0\\\
0&d_{28}&d_{38}&0&0&0&0&0\end{array}\right],$ (26)

where the twelve different types of nonzero elements are

$\begin{cases}d_{16}=2z_{6}\left(z_{2}z_{3}+z_{1}\tilde{z}_{1}\right),\\\
d_{17}=-2z_{7}\left(z_{2}z_{3}+z_{1}\tilde{z}_{1}\right),\\\
d_{25}=2z_{2}\left(z_{6}z_{7}+z_{5}\tilde{z}_{5}\right),\\\
d_{26}=2z_{2}z_{6}\left(z_{4}+z_{8}-z_{1}-z_{5}\right),\\\
d_{27}=2z_{2}z_{7}\left(z_{1}+z_{8}-z_{4}-z_{5}\right),\\\
d_{28}=-2z_{2}\left(z_{6}z_{7}+z_{8}\tilde{z}_{8}\right),\\\
d_{35}=-2z_{3}\left(z_{6}z_{7}+z_{5}\tilde{z}_{5}\right),\\\
d_{36}=2z_{3}z_{6}\left(z_{4}+z_{5}-z_{1}-z_{8}\right),\\\
d_{37}=2z_{3}z_{7}\left(z_{1}+z_{5}-z_{4}-z_{8}\right),\\\
d_{38}=2z_{3}\left(z_{6}z_{7}+z_{8}\tilde{z}_{8}\right),\\\
d_{46}=-2z_{6}\left(z_{2}z_{3}+z_{4}\tilde{z}_{4}\right),\\\
d_{47}=2z_{7}\left(z_{2}z_{3}+z_{4}\tilde{z}_{4}\right).\end{cases}$ (27)

One example that satisfies (1) for (26) is Corboz_bookchapter_2013

$B=\sqrt{\frac{iU}{2}}\left[\begin{array}[]{cccccccc}i\tilde{z}_{1}z_{1}&-\tilde{z}_{1}z_{1}&-iz_{2}z_{3}&z_{2}z_{3}&\tilde{z}_{1}z_{1}&i\tilde{z}_{1}z_{1}&-z_{2}z_{3}&-iz_{2}z_{3}\\\
-iz_{1}z_{2}&z_{1}z_{2}&i\tilde{z}_{4}z_{2}&-\tilde{z}_{4}z_{2}&\tilde{z}_{1}z_{2}&i\tilde{z}_{1}z_{2}&-z_{2}z_{4}&-iz_{2}z_{4}\\\
i\tilde{z}_{1}z_{3}&-\tilde{z}_{1}z_{3}&-iz_{3}z_{4}&z_{3}z_{4}&-z_{1}z_{3}&-iz_{1}z_{3}&\tilde{z}_{4}z_{3}&i\tilde{z}_{4}z_{3}\\\
-iz_{2}z_{3}&z_{2}z_{3}&i\tilde{z}_{4}z_{4}&-\tilde{z}_{4}z_{4}&-z_{2}z_{3}&-iz_{2}z_{3}&\tilde{z}_{4}z_{4}&i\tilde{z}_{4}z_{4}\\\
i\tilde{z}_{5}z_{5}&\tilde{z}_{5}z_{5}&-iz_{6}z_{7}&-z_{6}z_{7}&\tilde{z}_{5}z_{5}&-i\tilde{z}_{5}z_{5}&-z_{6}z_{7}&iz_{6}z_{7}\\\
-iz_{5}z_{6}&-z_{5}z_{6}&i\tilde{z}_{8}z_{6}&\tilde{z}_{8}z_{6}&\tilde{z}_{5}z_{6}&-i\tilde{z}_{5}z_{6}&-z_{6}z_{8}&iz_{6}z_{8}\\\
i\tilde{z}_{5}z_{7}&\tilde{z}_{5}z_{7}&-iz_{7}z_{8}&-z_{7}z_{8}&-z_{5}z_{7}&iz_{5}z_{7}&\tilde{z}_{8}z_{7}&-i\tilde{z}_{8}z_{7}\\\
-iz_{6}z_{7}&-z_{6}z_{7}&i\tilde{z}_{8}z_{8}&\tilde{z}_{8}z_{8}&-z_{6}z_{7}&iz_{6}z_{7}&\tilde{z}_{8}z_{8}&-i\tilde{z}_{8}z_{8}\end{array}\right].$
(28)

We remind the reader of the notation $\tilde{z}_{j}\equiv 1-z_{j}$ for so-
called hole-variables. The $B\in\mathbb{C}^{8\times 8}$ matrix (28) is only
one out of many possible ways to represent the second-order derivatives from
the FPE (2) into an SDE (3).

###  A.3 The mixed Fermi–Bose system for molecular dissociation

In this subsection, we explicitly present one possible system of SDEs for the
Hamiltonian (10) with two fermionic momentum modes. Using the operator mapping
for the fermionic GPSR Corney-fermionicII  combined with the bosonic
positive-$P$ representation GardinerBook1 , we obtain the $a_{j}$ and
$D_{j,k}$ in equation (2). If we order the $n=5$ phase-space variables as
$z_{1}=n_{1}$, $z_{2}=m_{1}$, $z_{3}=m_{1}^{+}$, $z_{4}=\alpha$,
$z_{5}=\alpha^{+}$, then (2) reduces to Ogren_EPL_2010

$\left[\begin{array}[]{c}a_{1}\left(\vec{z}\right)\\\
a_{2}\left(\vec{z}\right)\\\ a_{3}\left(\vec{z}\right)\\\
a_{4}\left(\vec{z}\right)\\\
a_{5}\left(\vec{z}\right)\end{array}\right]=\left[\begin{array}[]{c}\kappa\left(z_{3}z_{4}+z_{2}z_{5}\right)\\\
-2i\Delta_{1}z_{2}+\kappa z_{4}\left(1-2z_{1}\right)\\\
2i\Delta_{1}z_{3}+\kappa z_{5}\left(1-2z_{1}\right)\\\ -\kappa z_{2}\\\
-\kappa z_{3}\\\ \end{array}\right].$ (29)

The second-order differential operator in the FPE (2) is defined by the
symmetric $D\in\mathbb{C}^{5\times 5}$ matrix

$D=\kappa\left[\begin{array}[]{ccccc}0&0&0&d_{14}&d_{15}\\\
0&0&0&d_{24}&d_{25}\\\ 0&0&0&d_{34}&d_{35}\\\ d_{14}&d_{24}&d_{34}&0&0\\\
d_{15}&d_{25}&d_{35}&0&0\end{array}\right],$ (30)

where the six different types of nonzero elements are

$\begin{cases}d_{14}=z_{1}z_{2},\\\ d_{15}=z_{1}z_{3},\\\ d_{24}=z_{2}^{2},\\\
d_{25}=-z_{1}^{2},\\\ d_{34}=-z_{1}^{2},\\\ d_{35}=z_{3}^{2}.\end{cases}$ (31)

With the following $B\in\mathbb{C}^{5\times 4}$ matrix

$B=\frac{\sqrt{\kappa}}{\sqrt{2}}\left[\begin{array}[]{cccc}z_{1}z_{2}&-iz_{1}z_{2}&z_{1}z_{3}&-iz_{1}z_{3}\\\
z_{2}^{2}&-iz_{2}^{2}&-z_{1}^{2}&iz_{1}^{2}\\\
-z_{1}^{2}&iz_{1}^{2}&z_{3}^{2}&-iz_{3}^{2}\\\ 1&i&0&0\\\
0&0&1&i\end{array}\right],$ (32)

we have one possible $n\times m$ matrix that solve equation (1). For a more
compact explicit notation of the SDEs (3), we can first define a matrix
$\widetilde{B}$ with half as many columns, $\widetilde{B}_{j,k}=B_{j,2k-1}$.

Then each of the $m/2$ columns of $\widetilde{B}$ is multiplied with by one of
the $m/2$ complex noises $\xi_{1},\ldots,\xi_{m/2}$, or by its complex
conjugate if $\textnormal{Im}[B_{j,2k}]/\textnormal{Re}[B_{j,2k-1}]=-1$. With
a noise vector written with only 2 complex noises $\xi_{1}$ and $\xi_{2}$,
according to equation (9), we have the noise terms from Ogren_EPL_2010 , that
is,

$\vec{b}(\vec{z},\vec{\xi})=\sqrt{\kappa}\left[\begin{array}[]{c}z_{1}z_{2}\xi_{1}^{*}+z_{1}z_{3}\xi_{2}^{*}\\\
z_{2}^{2}\xi_{1}^{*}-z_{1}^{2}\xi_{2}^{*}\\\
-z_{1}^{2}\xi_{1}^{*}+z_{3}^{2}\xi_{2}^{*}\\\ \xi_{1}\\\
\xi_{2}\end{array}\right].$ (33)

It is the same procedure, but using 4 complex noises, that have been used in
going from equation (28) to the SDE (8) for the Fermi–Hubbard model.

##  Appendix B Numerical computation of $B$-matrices using Takagi’s
decomposition

We begin by recalling Takagi’s decomposition and the algorithms that can be
used to compute it. We then describe how this decomposition can be used to
construct a $B$-matrix satisfying (1).

###### Theorem B.1

(Chapter 4, Corollary 4.4.4(c) of HoJo13 ) Let $D$ be a complex symmetric
$n\times n$ matrix. Then there exist a unitary matrix $U$ and a nonnegative
diagonal matrix $\varSigma=\mathrm{diag}(\sigma_{1},\dots,\sigma_{n})$ such
that $D=U\varSigma U^{T}$.

Takagi’s decomposition can be interpreted as a symmetric analog of a Singular
Value Decomposition (SVD) of $D$. Note that $\sigma_{1},\dots,\sigma_{n}$ are,
in fact, the singular values of $D$ Thom79 .

Several algorithms to compute Takagi’s decomposition have been presented in
the literature GuQi03 ; WaLC18 ; XuQi08 ; BuGr88 . In this paper, we use the
divide-and-conquer method developed by Xu and Qiao XuQi08 . As this algorithm
requires the input matrix be tridiagonal, we use the Lanczos
tridiagonalization with modified partial orthogonalization and restart to
reduce a complex symmetric matrix to a tridiagonal form GuQi03 .

By Theorem B.1, we have that $D=U\varSigma
U^{T}=\bigl{(}U\sqrt{\varSigma}\bigr{)}\bigl{(}U\sqrt{\varSigma}\bigr{)}^{T}$
for some unitary $U$ and nonnegative diagonal $\varSigma$. If we define
$B:=U\sqrt{\varSigma}$, where
$\sqrt{\varSigma}=\mathrm{diag}(\sqrt{\sigma_{1}},\dots,\sqrt{\sigma_{n}})$,
then it is easy to see that we have obtained the desired decomposition
$D=BB^{T}$. Note that this decomposition is not unique: for any orthogonal
matrix $Q$, for example, we have
$D=BQQ^{T}B^{T}=(BQ)(BQ)^{T}=:B^{\phantom{T}}_{1}B_{1}^{T}$.

##  Appendix C Itô calculus, to avoid Stratonovich corrections

It is common to convert the Itô stochastic differential equations to
Stratonovich form and to integrate the latter with semi-implicit methods
Drummond_Mortimer_1991 , as these sometimes have better convergence
properties. This in general means that nonzero Stratonovich corrections (SC)
to the drift-vector $\vec{a}(\vec{z})$, of the form

$SC_{i}=-\frac{1}{2}\sum_{j,k}B_{j,k}\frac{\partial B_{i,k}}{\partial z_{j}},$
(34)

need to be added GardinerBook1 . These have been tested for the specific
analytic examples of $B$-matrices here, i.e., equations (23), (28), and (32).
However, in the case of numerical $B$-matrices we want to avoid the numerical
differentiations needed in (34), thus our chosen strategy was to only use Itô
calculus when applying numerical $B$-matrices.

## References

  * (1) Corney J F and Drummond P D 2006 Gaussian phase-space representations for fermions Phys. Rev. B 73 125112
  * (2) Corney J F and Drummond P D 2006 Gaussian operator bases for correlated fermions J. Phys. A: Math. Gen. 39, 269
  * (3) Drummond P D and Gardiner C W 1980 Generalised P-representation in quantum optics J. Phys. A 13 2353
  * (4) Gardiner C W 2008 Handbook of Stochastic Methods. Springer, Berlin, 4th ed.
  * (5) Plimak L I, Olsen M K and Collett M J 2001 Optimization of the positive-P representation for the anharmonic oscillator Phys. Rev. A 64 025801
  * (6) Deuar P and Drummond P D 2002 Gauge P representations for quantum-dynamical problems: Removal of boundary terms Phys. Rev. A 66 033812
  * (7) Smith A M and Gardiner C W 1989 Simulations of nonlinear quantum damping using the positive $P$ representation Phys. Rev. A 39 3511
  * (8) Gilchrist A, Gardiner C W and Drummond P D 1997 Positive P representation: Application and validity Phys. Rev. A 55 3014
  * (9) Assad F F, Werner P, Corboz P, Gull E and Troyer M 2005 Symmetry projection schemes for Gaussian Monte Carlo methods Phys. Rev. B 72 224518
  * (10) Corboz P, Troyer M, Kleine A, McCulloch I P and Schollwöck U 2008 Systematic errors in Gaussian quantum Monte Carlo and a systematic study of the symmetry projection method Phys. Rev. B 77 085108
  * (11) Aimi T and Imada M 2007 Gaussian-basis monte carlo method for numerical study on ground states of itinerant and strongly correlated electron systems, J. Phys. Soc. Jpn. 76, 084709
  * (12) Ögren M, Kheruntsyan K V and Corney J F 2010 First-principles quantum dynamics for fermions: application to molecular dissociation Europhys. Lett. 92 36003
  * (13) Ögren M, Kheruntsyan K V and Corney J F 2011 Stochastic simulations of fermionic dynamics with phase-space representations Comp. Phys. Comm. 182 1999
  * (14) Corboz P, Ögren M, Kheruntsyan K V and Corney J F 2013 Phase-Space Methods for Fermions in Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics, Imperial College Press London, pp. 407-416. https://doi.org/10.1142/9781848168121_0027
  * (15) Walls D F and Milburn G J 2008 (2.nd edition) Quantum Optics (Springer-Verlag Berlin Heidelberg)
  * (16) Ögren M and Carlsson M 2013 On the dynamics of the Fermi-Bose model J. Phys. A:Math. Theor. 46, 015005
  * (17) Steel M J, Olsen M K, Plimak L I, Drummond P D, Tan S M, Collett M J, Walls D F and Graham R 1998 Dynamical quantum noise in trapped Bose-Einstein condensates Phys. Rev. A 58 4824
  * (18) Pethick C J and Smith H 2008 Bose–Einstein Condensation in Dilute Gases 2nd edn (Cambridge: Cambridge University Press)
  * (19) Rahav S and Mukamel S 2009 Gaussian phase-space representation of fermion dynamics: Beyond the time-dependent Hartree-Fock approximation Phys. Rev. B 79 165103
  * (20) Jack M W and Pu H 2005 Dissociation dynamics of a Bose-Einstein condensate of molecules Phys. Rev. A 72, 063625
  * (21) Wüster S, Corney J F, Rost J M and Deuar P 2017 Quantum dynamics of long-range interacting systems using the positive-P and gauge-P representations Phys. Rev. E 96, 013309
  * (22) Drummond P D and Mortimer I K 1991 Computer Simulations of Multiplicative Stochastic Differential Equations J. Comput. Phys. 93 144
  * (23) Horn R A and Johnson C R 2013 Matrix analysis. Cambridge University Press, Cambridge, second edition
  * (24) Guo C and Qiao S 2003 A stable Lanczos tridiagonalization of complex symmetric matrices, Technical report
  * (25) Thompson R C 1979 Singular values and diagonal elements of complex symmetric matrices Linear Algebra and its Applications, 26:65–106
  * (26) Wang X, Liang L, and Che M 2018 Iterative algorithms for computing the Takagi factorization of complex symmetric matrices IAENG International Journal of Applied Mathematics, 48:297–305, 08
  * (27) Xu W and Qiao S 2008 A divide-and-conquer method for the Takagi factorization SIAM Journal on Matrix Analysis and Applications, 30(1):142–153
  * (28) Bunse-Gerstner A and Gragg W B 1988 Singular value decompositions of complex symmetric matrices Journal of Computational and Applied Mathematics, 21(1):41–54


Generating spectral dental panoramic images from computed tomography volumes
Villseveri Somerkivi1,2,3, Anna Kolb1,2, Thorsten Sellerer1,2, Daniel
Berthe1,2, York Hämisch4, Tuomas Pantsar4, Henrik Lohman4, Franz
Pfeiffer1,2,5,6
1 Chair of Biomedical Physics, Department of Physics, School of Natural
Sciences, Technical University of Munich, 85748 Garching, Germany
2 Munich Institute of Biomedical Engineering, Technical University of Munich,
85748 Garching, Germany
3 Planmeca Oy, 00880 Helsinki, Finland
4 Varex Imaging Corp, 84104 Salt Lake City, USA.
5 Department of Diagnostic and Interventional Radiology, Klinikum rechts der
Isar, Technical University of Munich, 81677 Munich, Germany
6 Institute for Advanced Study, Technical University of Munich, 85748
Garching, Germany
Version typeset

Author to whom correspondence should be addressed. email:
<EMAIL_ADDRESS>

###### Abstract

Background: Spectral panoramic X-ray imaging is an emerging imaging modality
with possible clinical benefits. With preliminary experimental results
demonstrated, a method enabling verification and cross-comparison would
increase confidence in the experimental results.
Purpose: To generate synthetic spectral panoramic images and compare the
synthetic images against experimental spectral panoramic images.
Methods: A simulation framework for generating synthetic spectral panoramic
images from CT volumes is described. A cone beam CT scan of an anthropomorphic
head phantom is used as input. An experimental spectral panoramic image of the
same phantom is acquired. The framework also includes a tool for estimating
anatomical noise.
Results: The output of the framework of an anthropomorphic head phantom is
compared against an experimental spectral panoramic image of the same phantom.
The synthetic and experimental spectral panoramic images resemble each other
considerably, especially the bone features. In the soft tissue images, there
are some differences, which are a result of the differences in how each image
is generated. Finally, anatomical noise is revealed to be small, with the
spread in the values in the order of 1 cm for the soft tissue and 1 mm for the
cortical bone through most of the image.
Conclusions: It is demonstrated that generating synthetic spectral panoramic
images from CT volumes is possible. The synthetic images have many
similarities with the experimental results, increasing the confidence in the
correctness of the information contained within experimental spectral
panoramic images and indicating that the synthetic images could be useful in
further research.

The table of contents is for drafting and refereeing purposes only. Note that
all links to references, tables and figures can be clicked on and returned to
calling point using cmd[ on a Mac using Preview or some equivalent on PCs (see
View - go to on whatever reader).

###### Contents

  1. I. Introduction
  2. II. Methods
    1. II.A. Simulation framework
    2. II.B. Experimental measurements
      1. II.B.1. CBCT acquisition
      2. II.B.2. Spectral panoramic acquisition
  3. III. Results
  4. IV. Discussion
  5. V. Conclusion
  6. References

##  I. Introduction

Dental panoramic X-ray imaging is a common imaging modality around the world.
According to United Nations Scientific Committee on the Effects of Atomic
Radiation, the panoramic imaging examination frequency per 1000 population was
140 in Germany, 120 in Japan, 74 in the United Kingdom, and 64 in the United
States 1. It provides a low-dose radiographic image, which is used mainly for
a high- resolution view of the maxilla, the mandible, and the teeth and any
lesions in these structures 2. A panoramic image is generated by moving and
rotating the X-ray tube and the detector in parallel to the dental arch while
maintaining a constant magnification of the arch, as shown in fig. 1. In the
case of a digital detector, frames are captured at a high frame rate and
summed with a partial overlap to generate the panoramic image 3. Different
tissue types are superimposed in conventional, transmission-based radiography,
including panoramic imaging. This superimposition causes difficulties and
confusion in the interpretation of panoramic images. For example, the oral
cavity can be mistaken for a fracture in the mandible, even though the cavity
is in the soft tissue 4.

Figure 1: The working principle of panoramic imaging. Shown are the
trajectories of the image receptor and the X-ray source. The position of these
two decides the focal area (also called the sharp layer) and the rotation
center, which in this case is sliding. The focal area aligns with the teeth.

Spectral X-ray imaging has been studied extensively in the past 15 years 5 6 7
8. These studies have demonstrated that material decomposition provides
additional information on top of the normal attenuation images by revealing
clinical features which would be invisible in a conventional attenuation
image. Recently, spectral panoramic imaging specifically has also been
investigated 9 10 11. A previous work demonstrated that it is possible to
achieve material decomposition in panoramic imaging for as long as there is a
sufficient number of counts for calibration and decomposition 11.

Currently, no clinical studies on the feasibility of spectral panoramic
imaging have been published. While the results from a previous pre-clinical
study 11 seem plausible, the confidence in the results would increase if there
was a reference to compare these results to. Several approaches to generate
synthetic panoramic images from CT have been described 12 13 14 15. These
approaches aim at automatically or semi-automatically finding the dental arch,
e.g. by using the maximum intensity projection along the axial slices followed
by subsequent thresholding and morphological thinning. After the dental arch
has been extracted, the panoramic image is generated by extracting voxel
values along lines which are orthogonal to the dental arch. The length of
lines is set so that the lines run through the mandible, the maxilla, and the
teeth while the rest of the lower head is left out. This kind of process may
be suitable for generating synthetic panoramic images for reading in dental
practice. However, if verification is desired, then the whole image generation
process should be modeled. This includes forward projecting actual frames and
summing the frames to get the panoramic image. Additionally, the generation of
basis material images is required when modeling spectral imaging.

This work introduces a framework for generating spectral panoramic images from
CT volumes. We use a CT image of an anthropomorphic head phantom as input and
generate a set of spectral panoramic images using the framework. We then
cross-compare the output of the framework against experimental spectral
panoramic images.

##  II. Methods

###  II.A. Simulation framework

The framework was implemented in Python (Python Software Foundation,
Wilmington, Delaware, United States). A proprietary X-AID module (MITOS GmbH,
Munich, Germany) was used for the forward projection step. The main steps of
the workflow are shown in fig. 2 and the detailed explanation is below.

Figure 2: Main steps of the simulation framework. The dental arch used and the
geometry generated from the arch are shown overlaid with an axial slice of a
volume on the bottom left side. Every 60th projection is shown, with green
being the detector, dark blue the source, cyan the focal through, and red the
rotation center.

The first part of the framework deals with preparing the input CT volume for
the workflow. If the voxels are not isotropic, an interpolation is first
carried out so that the output is isotropic. Secondly, the user has to rotate
the volume along the three axes so that the volume is symmetric along the
coronal and axial views, and the dental arch is horizontal in the sagittal
view. The final step of the preparation is thresholding the volume to two
basis materials. There are two ways to carry out the thresholding, depending
on the type of the input volume.

A simple threshold is used in the segmentation if the imaging system’s
Hounsfield Unit (HU) values are not perfectly reliable due to scatter
contamination 16. A threshold of 1280 HU was used in this study. Both basis
material images are scaled with their respective mean signals so that the mean
value of the segmented volumes is 1.

If the HU values of the input volume are more reliable e.g., by reducing the
scatter contamination by using an anti-scatter grid 17, then a more
sophisticated thresholding is possible in the framework. The method uses a
sheetness score combined with a threshold in order to better segment the bone,
as bones in the skull mainly consist of sheet-like structures 18. The
sheetness score $S$ is calculated as in 18:

$S=\text{
exp}\left[\frac{-R_{\text{sheet}}^{2}}{2\alpha^{2}}\right]\left(1-\text{exp}\left[\frac{-R_{\text{tube}}^{2}}{2\beta^{2}}\right]\right)\left(1-\text{exp}\left[\frac{-R_{\text{noise}}^{2}}{2\gamma^{2}}\right]\right),$
(1)

where $R_{\text{sheet}}=|\lambda_{2}|/|\lambda_{3}|$,
$R_{\text{tube}}=|(2|\lambda_{3}|-|\lambda_{2}|)-|\lambda_{1}|)|/|\lambda_{3}|$,
and $R_{\text{noise}}=\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}}$
with $\lambda_{i}$ being the local Hessian eigenvalues so that
$|\lambda_{1}|\leq|\lambda_{2}|\leq|\lambda_{3}|$. The output HU values
$\Omega$ are then:

$\Omega=\begin{cases}-1000&\text{if }(\Omega_{i}S\cdot 5+\Omega_{i})\leq T\\\
\Omega_{i}S&\text{othwerwise}\\\ \end{cases},$ (2)

where $\Omega_{i}$ are the HU values of the input volume and $T$ is the
threshold value. $(\Omega_{i}S\cdot 5+\Omega_{i})$ was heuristically found to
yield good results in the thresholding.

The next step is to trace the dental arch of the volume. The dental arch needs
to be correctly registered in order for the panoramic reconstruction to
succeed (c.f. fig 1). This requires the user to manually provide the
coordinates of the dental arch. There are two ways to achieve this in the
framework: u-shaped (9 coordinates) and exact (11 coordinates) arches. The
u-shaped arch generates a smoother profile with a smaller number of rotations.
The u-shaped arch was used in this work (shown in fig. 2) as it corresponds
better to a realistic trajectory. It is also possible to restrict the overlap
of the mandible from the other side in the framework. This feature was not
used for the results presented in this paper, as it was desirable to have the
simulation results match the experimental results. The dental arch is created
by interpolating between the user-given points by using the piecewise cubic
hermite interpolating polynomial (PCHIP) algorithm. The algorithm is desirable
for creating the dental arch, as it is more resistant to ringing artifacts
than spline or regular polynomial-based interpolation algorithms.

After the dental arch is complete, the subsequent step is forward projecting
the frames. The X-AID module uses a ray-driven forward projector which traces
the path from the source to the target pixel and sums up everything along the
path. The panoramic reconstruction is carried out by summing the frames
according to the movement profile generated using the dental arch (c.f. fig.
2).

The output of the panoramic reconstruction is called a panoramic equivalent
thickness image (PETI) in this framework. In this case, the volume was
thresholded into soft tissue and cortical bone. Thus, the output here is two
PETIs of soft tissue and cortical bone with a unit of thickness (mm by
default).

For the generation of log images, cortical bone and soft tissue PETIs were
used as basis materials. The linear attenuation coefficients (LACs) are
obtained from 19. The pyXSFW framework 20 was used to model the tungsten
spectrum. The detector was modeled as having a sensor thickness of 0.75 mm of
cadmium telluride, with each absorbed photon inducing one count. With the
spectrum and the LAC(s) known for both basis materials and the detector, the
output signal can be calculated:

$I(E)=\sum_{E=1}^{\infty}I_{0}(E)D(E)F(E)e^{-t_{st}\mu_{st}(E)-t_{cb}\mu_{cb}(E)},$
(3)

where $I_{0}(E)$ is the input spectrum, $D(E)$ is the detector absorption
efficiency, $F(E)$ is the spectral response of the detector, $t_{st}$ and
$t_{cb}$ are the respective soft tissue and cortical bone line integral
lengths, while $\mu_{st}(E)$ and $\mu_{cb}(E)$ are the energy-dependent soft
tissue and cortical bone LACs, respectively.

No quantum noise was added to the images in this work, as the purpose was to
evaluate the feasibility of the results in the sense that the output is as
expected. While the experimental results are subject to quantum noise, it
would not make the comparison easier or better to add noise to the synthetic
images. If there was a need to add quantum noise, it should be added to the
forward projected frames before panoramic reconstruction. This would ensure
that the reconstructed image has a realistic noise distribution.

Virtual monoenergetic images (VMIs) are generated by multiplying the PETIs
with their respective LACs and then summing them together.

Anatomical noise in a panoramic image pixel is measured by recording all of
the pixel values from the frames which contribute to the pixel in the
panoramic image. In other words, this is an indication of how different the
individual line integrals paths making up the final panoramic image pixel are.
The mean of these values is the final PETI output value in the pixel.

###  II.B. Experimental measurements

####  II.B.1. CBCT acquisition

Planmeca Viso G7 (Planmeca Oy, Helsinki, Finland) was used for the cone beam
CT (CBCT) acquisition. The acquisition was carried out in offset mode, and two
rotations were conducted and stitched together to acquire the whole head
phantom. The reconstructed volume was used as input to the simulation
framework. The acquisition parameters are listed in table 1. The dose is
higher than what would be used in the typical patient acquisition, as high
resolution with lower noise was desired.

Table 1: CBCT acquisition parameters.

source to image distance | 70 cm
---|---
acceleration voltage | 120 kVp
tube current | 11 mA
exposure time | 11.4 s
dose-area product | 4100 mGy(cm)2
magnification | 1.5

The anthropomorphic head phantom RANDO SK150 (Radiation Analogue Dosimetry
System; The Phantom Laboratory, Salem, NY, USA) was used in this study. The
phantom is a skull encapsulated in plastic. The plastic has air cavities that
mimic the upper respiratory tract. A hole has been drilled into the phantom
from the bottom of the jaw so that there is also a cavity in the bottom of the
jaw in the plastic and in the mandible.

####  II.B.2. Spectral panoramic acquisition

For the experimental data, a low-order polynomial model of line integrals for
polychromatic X-ray spectrum is used 21. The expected number of counts
$\hat{y}_{i}^{s}$ is then

$\hat{y}_{i}^{s}=b_{i}^{s}e^{-\hat{l}_{i}^{s}},\
\hat{l}_{i}^{s}=P_{i}^{s}(A_{i};c_{is}),\
A_{i}={(A_{i}^{0},...,A_{i}^{n})}^{T},$ (4)

where $b_{i}^{s}$ is the flatfield signal for pixel $i$, spectral measurement
$s$ and $A_{i}$ are the basis material line integrals. For this paper, the
polynomial model $P_{i}^{s}$ was chosen to have two basis materials:

$\begin{split}P_{i}^{s}&=\frac{N}{D}\\\
N&=c_{0}+c_{1}A_{i}^{0}+c_{2}A_{i}^{1}+c_{3}(A_{i}^{0})^{2}+c_{4}(A_{i}^{1})^{2}+c_{5}A_{i}^{0}A_{i}^{1}\\\
D&=1+c_{6}A_{i}^{0}+c_{7}A_{i}^{1}\end{split},$ (5)

where $c$ are the fit coefficients. Here, the subscript is the number of the
fit coefficient, while pixel and spectral measurement subscripts are omitted.

Calibration measurements are least-squares fitted to establish the fit
parameters $c_{is}$:

$c_{is}=\text{arg
min}\sum_{k=1}^{K}w_{k}{\left(l_{ik}^{s}-P\left(A_{ik};c_{is}\right)\right)}^{2},$
(6)

where $w_{k}$ are the weights with $k$ being the index of the calibration
measurement, $l_{ik}^{s}$ the linearized calibration measurements, and
$A_{ik}$ the line integrals. The standard deviation of the particular
calibration measurement determines the weights.

After the forward model has been calibrated, a maximum likelihood estimation
is carried out to find the basis materials thicknesses $A_{i}$, which maximize
the likelihood given the spectral measurements. Assuming independent Poisson
statistics for the spectral measurements, minimizing the corresponding
negative log-likelihood is computationally more efficient:

$\mathcal{A}_{i}=\text{arg
min}-\mathcal{L}_{i}(A_{i})=\sum_{s=1}^{S}\hat{y}_{i}^{s}(A_{i})-y_{i}^{s}\text{ln}(\hat{y}_{i}^{s}(A_{i})).$
(7)

A dictionary learning-based denoising is also applied to the images, with
details provided in 22 11.

Panoramic imaging presents a new challenge to the decomposition process, as
individual frames in panoramic imaging have a very low number of counts.
Directly decomposing frames would induce severe statistical bias to the
decomposition 23. On the other hand, the count statistics in the final
panoramic image are sufficient for typical line integrals in the jaw area 11.
To decompose the whole panoramic image instead of individual frames, an
appropriate calibration needs to be carried out, however. One such way is
generating virtual panoramic images for calibration. This involves taking
still calibration exposures, taking the mean of each measurement, and
performing the panoramic reconstruction by using the mean frame as every
frame. The next step is calibrating the panoramic image using these virtual
panoramic images. Thus, the actual panoramic exposure will have the same
flatfield values and pixel weights as the calibrated panoramic image. A change
of basis from the calibration materials into soft tissue and dentin is carried
out after the decomposition. A detailed explanation of the whole process can
be found in 11.

An experimental spectral photon counting panoramic system was used to capture
an image of the same phantom. The Planmeca Promax Mid (Planmeca Oy, Helsinki,
Finland) system has a DC-Vela detector (Varex Imaging Corporation, Salt Lake
City, USA). DC-Vela is a photon counting detector with two adjustable
thresholds, a charge sharing correction, and a 0.75 mm cadmium telluride
sensor layer. The acquisition parameters for the experimental panoramic system
are shown in table 2.

Table 2: Panoramic acquisition parameters.

source to image distance | 60 cm
---|---
acceleration voltage | 70 kVp
tube current | 9 mA
exposure time | 18 s
dose-area product | 72.1 mGy(cm)2
magnification | 1.4

A decomposition process like the one described in this section was not modeled
in the simulation framework. If the same forward model is used to simulate the
transmission and to decompose, the output of the decomposition will be the
same as the input in the absence of noise. This holds as long as the input
volume is segmented into exactly two basis materials.

##  III. Results

The output of the simulation framework is shown in fig. 3. As expected, the
maxilla, the mandible and the teeth are the prominent features of the cortical
bone image (a). The spine is also visible to the sides. The soft tissue image
(b) has larger line integral thicknesses to the center. Outside of the center,
there are homogeneous areas with, darker air gaps with the upper respiratory
tract being the most prominent. The tract starts at the sides, next to the
spine, and extends upwards until splitting into the oral and nasal tracts and
arching towards the center. There are also darker areas in the soft tissue
image (b) in the same places where the bone features are in the cortical bone
image (a). For example, the spine, the mandible, and the teeth are visible as
darker areas in the soft tissue image (b). The drill hole (solid orange
ellipse) is visible in both images. This is expected, as the hole was indeed
drilled through the bottom of the jaw and partially through the mandible. The
air gaps highlighted in (b) by the larger and smaller dashed ellipses are not
part of the respiratory tracts. The gaps are roughly where the maxillary
sinuses would be. While darker areas are inside the ellipses, the areas are
only barely visible.

Figure 3: Synthetic pan images, generated from the CBCT scan of the RANDO
SK150 phantom, with window for cortical bone image (a) [0.00–2.28 cm] and for
soft tissue (b) [7.30–17.2 cm]. A solid orange ellipse highlights the drilled
hole in both images. The blue dashed lines highlight the boundary for the
overlapping mandible from the other side in the bone image (a). The two dashed
ellipses show small air gaps in the plastic in the soft tissue image (b). Many
dark features in the soft tissue image (b) are not air but the absence of bone
i.e., there is a strong negative correlation between the two images.

The basis material images from the experimental system are shown in fig. 4.
The maxilla, the mandible, and the teeth are the dominant features in the
dentin image (a). The soft tissue image (b) is quite homogeneous, with the air
gaps, especially the upper respiratory tract being the dominant features. The
air gaps not being part of the respiratory tract are again highlighted by the
dashed ellipses. In the experimental image 4, these gaps stand out more
prominently, although the smaller gap overlapping with the lower nasal cavity
limits its visibility.

Figure 4 is similar to what was described in a previous work 11. The most
notable difference is the presence of the upper respiratory tract in the soft
tissue image. Note that the images of this work are left-right flipped, unlike
the image in the previous work. Dentin was considered to be a good basis
material in the previous work, so it was kept the same in this work (instead
of using cortical bone as in the synthetic image).

Figure 4: The experimental panoramic images, the window for the dentin image
(a) is [0.00–1.92 cm] and for the soft tissue image (b) [4.22–11.1 cm]. The
highlighted details are the same as in fig. 3. The overlap from the mandibles
is not symmetric in (a), as the overlap on the right side of the image only
starts from a higher and more away from the center than the overlap on the
left side. Here, the darker features in (b) are mostly air gaps in the
phantom. The yellow and the red points are the anatomical noise measurement
pixels.

There are some differences in features between the bone component in figures 3
and 4. The mean line integral length of fig.3 (a) is lower than fig. 4 (a).
The overlap from the mandible from the other side has locally very large
values in fig. 3, however. This is why the windowing is wider in fig. 3 (a)
than in fig. 4 (a). In fig. 3 (a), the contrast between the body and the base
of the mandible is more prominent than in 4 (a).

The soft tissue images (b) have notable differences in 3 and 4. The soft
tissue image in 4 is much more homogeneous than 3. More specifically, the
bright features in the synthetic cortical bone image (a) in fig. 3 appear dark
in the soft tissue image (b) in fig. 3. In other words, the negative
correlation between the experimental basis material images in 4 is much
smaller than in the synthetic images in fig. 3. On the other hand, the upper
respiratory tract is clearly visible in both soft tissue images. Additionally,
the drill hole can be seen in both the dentin images and both soft tissue
images in figures 3 and 4.

The shape of one of the overlapping mandibles is highlighted by a blue dashed
line in figures 3 and 4. But the exact location and shape of the overlap is
different between the two figures. This is a result of different imaging
geometries. The experimental device uses a fixed geometry, whereas the
simulator uses a manually defined geometry, which is specifically tuned for
each volume.

The details in fig. 5 show slices from the CBCT acquisition. The slices
confirm that there are indeed air gaps in the positions highlighted in fig. 3
and fig. 4. The drill hole is inside the mandible but is visible in all of the
images in figures 3 and 4.

Figure 5: An axial and two coronal slices showing the same details in the CBCT
image as in the two pan images. From the CBCT image, it is confirmed that
there are indeed air gaps in the highlighted positions. The visible soft
tissue boundary in (b) and (c) around eye sockets’ upper part marks the two
rotations’ stitching boundary. The windowing is [0–2300 HU] for all images.

Figure 6 displays the anatomical noise for both the soft tissue and the
cortical bone PETIs of two pixels in fig. 3. More precisely, each histogram in
fig. 6 is generated by recording the individual projections contributing to
the panoramic image’s particular pixel. These two cases were included to
highlight the extremes in the anatomical noise. The upper two histograms are
near the boundary of the mandible overlap. The values close to 18 mm (b) are a
sum of the overlap and the actual signal, while the lower values are only from
the actual signal. The lower two histograms are from a region with little
anatomical noise. Thus, the values are very close to each other. Most of the
panoramic image has a noise magnitude close to pixel (900, 2390), while an
increase in the noise is observed for bone – soft tissue boundaries with pixel
(1131, 1180) as an extreme case with the mandible overlap.

Figure 6: Cortical bone and soft tissue anatomical noise histograms from two
locations with pixel indices given in (height, width). The yellow point in
fig. 3 corresponds to the top row (a) and (b) and the red point (c) and (d) to
the bottom row. Observe that the x-scales are very different for the two rows
(compare (a) vs (c) and (b) vs (d)).

##  IV. Discussion

The presented framework enables generation of panoramic images from CT
volumes. The PETIs generated this way resemble the experimental material
decomposition panoramic images. The framework also enables quantitative and
pixel-wise estimation of anatomical noise in panoramic images.

The thresholding approach for segmentation assumes that one voxel induces a
signal only on the soft tissue or bone image, but not both. The assumption may
be correct for voxels containing only plastic, as they should only contribute
to the soft tissue image. The skull, on the other hand, consists of more than
one type of tissue. For instance, cancellous bone 24, cortical bone 24, dentin
25, and enamel 25 have different effective atomic numbers and are to an extent
present in the skull. Furthermore, in reality, one voxel can have more than
one tissue type.

The overlap of the mandible from the other side appears different in figures 3
and 4 due to the different geometries. The difference has an impact on both
location and the intensity of the overlap. Ideally, the geometry of the
experimental setup should also be used in the simulation framework to mitigate
these differences. However, the exact geometry of the used imaging system is
not publicly available. The non-symmetrical overlap of the mandibles is likely
caused by a sub-optimal mounting of the detector in the experimental system.

The choice of PETI materials is arbitrary when segmenting CT volumes. Cortical
bone and soft tissue were used in this study as they are a good guess with
respect to the effective atomic numbers and densities of the underlying
materials. The experimental setup decomposes images initially into the
calibration materials, in this case, titanium and polyoxymethylene (POM). In
the subsequent change of basis it has to be decided which materials and
energies for low and high bins are used, which are again arbitrary to some
extent. In other words, there is no single correct value for the experimental
output images. Thus, it is not feasible to expect the synthetic and
experimental panoramic images to match exactly.

The input volumes should ideally have both high-resolution and reliable HU
values. This would enable reliable two-material segmentation while preserving
bone feature visibility. Reliable HU-values would also increase the volume of
segmented images, as with the system used in this study some of the cancellous
bone was segmented as soft tissue. Partially missing cancellous bone is a
plausible explanation as to why the contrast between the body and the base of
the mandible is stronger in fig. 3 than in fig. 4. An increase in the soft
tissue – bone threshold would result in soft tissue being segmented as bone.
With reliable HU values, the bone could be potentially segmented into more
than two basis materials, making the results even more realistic.

With the aforementioned differences in the generation of the experimental and
synthetic spectral panoramic images in mind, it was expected that the two
images are not exactly the same. Nevertheless, the close resemblance of the
two images increases the confidence that the output of the experimental
spectral panoramic imaging workflow is feasible and expected. Additionally,
the developed simulation framework will be useful in further exploring
spectral panoramic imaging.

The framework will be adapted for dual-energy CT (DECT) input in a future
work. DECT images can be decomposed and forward projected directly without
thresholding. This way, a voxel can consist of more than one material. DECT
decomposed images can also be used as a reference for estimating bone mineral
density (BMD) 26. Photon counting CT (PCCT) provides the same benefits as DECT
while having spatial resolution close to CBCT 27. Thus, PCCT would be the
ideal way to provide the ground truth for BMD estimation from spectral
panoramic images. Another future work will explore the clinical relevance of
PETIs.

##  V. Conclusion

A framework for generating synthetic spectral panoramic images from CT volumes
has been successfully implemented. The synthetic panoramic images obtained
using the framework support the results obtained from the experimental
spectral panoramic system. The framework also provides a basis for further
research.

## Acknowledgements

The authors acknowledge financial support through the EQAP project, which is
part of the One-Munich Strategy. The authors acknowledge Varex Imaging
Corporation for the hardware support with the photon counting detector and
Planmeca Oy for the support with the experimental setups. The authors
acknowledge Ari Hietanen and Juha Koivisto for helpful discussions.

## Conflicts of interest

Villeseveri Somerkivi is an employee of Planmeca Oy. York Hämisch, Tuomas
Pantsar and Henrik Lohman are employees of Varex Imaging Corp.

## References

##

  * 1 United Nations Scientific Committee on the Effects of Atomic Radiation,  Sources, effect and risks of ionizing radiation 2020/2021, volume I,  Unscear report, United Nations, New York, US, 2022.
  * 2 J. M. Ianucci and L. J. Howerton,  Dental Radiography: Principles and Techniques,  Elsevier-Saunders, St. Louis, MO, USA, 4 edition, 2012.
  * 3 W. D. McDavid, U. Welander, S. Brent Dove, and G. Tronjje,  Digital imaging in rotational panoramic radiography,  Dentomaxillofacial Radiology 24, 68–75 (1995).
  * 4 S. Perschbacher,  Interpretation of panoramic radiographs,  Australian dental journal 57, 40–45 (2012).
  * 5 J. Schlomka et al.,  Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,  Physics in medicine and biology 53, 4031 (2008).
  * 6 R. E. Alvarez,  Estimator for photon counting energy selective X-ray imaging with multibin pulse height analysis,  Medical physics 38, 2324–2334 (2011).
  * 7 S. Ehn, T. Sellerer, K. Mechlem, A. Fehringer, M. Epple, J. Herzen, F. Pfeiffer, and P. Noël,  Basis material decomposition in spectral CT using a semi-empirical, polychromatic adaption of the Beer–Lambert model,  Physics in Medicine and Biology 62, N1–N17 (2017).
  * 8 T. Flohr, M. Petersilka, A. Henning, S. Ulzheimer, J. Ferda, and B. Schmidt,  Photon-counting CT review,  Physica Medica 79, 126–136 (2020).
  * 9 T. Fukui, A. Katsumata, K. Ogawa, and S. Fujiwara,  Usefulness of an energy-binned photon-counting X-ray detector for dental panoramic radiographs,  in Medical Imaging 2015: Physics of Medical Imaging, volume 9412, page 94125O, International Society for Optics and Photonics, 2015.
  * 10 R. Langlais, A. Katsumata, S. Naidoo, K. Ogawa, T. Fukui, S. Shimoda, and K. Kobayashi,  The cadmium telluride photon counting sensor in panoramic radiology: gray value separation and its potential application for bone density evaluation,  Oral and maxillofacial radiology 120, 636–643 (2015).
  * 11 V. Somerkivi, T. Sellerer, T. Pantsar, H. Lohman, and F. Pfeiffer,  Spectral photon counting for panoramic dental imaging,  Biomedical Physics & Engineering Express 9, 035017 (2023).
  * 12 M. Bae, J.-W. Park, and N. Kim,  Semi-automatic and robust determination of dental arch form in dental cone-beam CT with B-spline approximation,  Computer Methods and Programs in Biomedicine 172, 95–101 (2019).
  * 13 T. Luo, C. Shi, X. Zhao, Y. Zhao, and J. Xu,  Automatic Synthesis of Panoramic Radiographs from Dental Cone Beam Computed Tomography Data,  PLOS ONE 11, 1–20 (2016).
  * 14 T. K. Papakosta, A. D. Savva, T. L. Economopoulos, G. K. Matsopoulos, and H. Gröhndal,  An automatic panoramic image reconstruction scheme from dental computed tomography images,  Dentomaxillofacial Radiology 46, 20160225 (2017).
  * 15 Z. Yun, S. Yang, E. Huang, L. Zhao, W. Yang, and Q. Feng,  Automatic reconstruction method for high-contrast panoramic image from dental cone-beam CT data,  Computer Methods and Programs in Biomedicine 175, 205–214 (2019).
  * 16 A. Sisniega, W. Zbijewski, A. Badal, I. S. Kyprianou, J. W. Stayman, J. J. Vaquero, and J. H. Siewerdsen,  Monte Carlo study of the effects of system geometry and antiscatter grids on cone-beam CT scatter distributions,  Medical Physics 40, 051915 (2013).
  * 17 C. Altunbas, Y. Park, Z. Yu, and A. Gopal,  A unified scatter rejection and correction method for cone beam computed tomography,  Medical physics 48, 1211–1225 (2021).
  * 18 O. Cuadros Linares, J. Bianchi, D. Raveli, J. Batista Neto, and B. Hamann,  Mandible and skull segmentation in cone beam computed tomography using super-voxels and graph clustering,  The Visual Computer 35, 1461–1474 (2019).
  * 19 J. H. Hubbell and S. M. Seltzer,  Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients 1 keV to 20 MeV for elements Z=1 to 92 and 48 additional substances of dosimetric interest,  Technical report, NIST, 1995.
  * 20 A. Malecki, G. Potdevin, and F. Pfeiffer,  Quantitative wave-optical numerical analysis of the dark-field signal in grating-based X-ray interferometry,  Europhysics Letters 99, 48001 (2012).
  * 21 H. N. Cardinal and A. Fenster,  An accurate method for direct dual-energy calibration and decomposition,  Medical physics 17, 327–341 (1990).
  * 22 K. Mechlem, T. Sellerer, S. Ehn, D. Münzel, E. Braig, J. Herzen, P. B. Noël, and F. Pfeiffer,  Spectral angiography material decomposition using an empirical forward model and a dictionary-based regularization,  IEEE Transactions on Medical Imaging 37, 2298–2309 (2018).
  * 23 K. Mechlem, S. Ehn, T. Sellerer, E. Braig, D. Münzel, F. Pfeiffer, and P. B. Noël,  Joint statistical iterative material image reconstruction for spectral computed tomography using a semi-empirical forward model,  IEEE transactions on medical imaging 37, 68–80 (2017).
  * 24 P. Hasgall, F. Di Gennaro, C. Baumgartner, E. Neufeld, B. Lloyd, M. Gosselin, D. Payne, A. Klingenböck, and N. Kuster,  IT’IS Database for thermal and electromagnetic parameters of biological tissues,  Technical Report Version 4.0, Foundation for Research on Information Technologies in Society, 2018.
  * 25 M. A. F. Zenóbio, M. S. N. Tavares, E. G. Zenóbio, and T. A. Silva,  Elemental composition of dental biologic tissues: study by means of different analytical techniques,  Journal of Radioanalytical and Nuclear Chemistry 289, 161–166 (2011).
  * 26 A. Laugerette et al.,  DXA-equivalent quantification of bone mineral density using dual-layer spectral CT scout scans,  European Radiology 29, 4624–4634 (2019).
  * 27 F. S. L. Thomsen, S. Horstmeier, J. H. Niehoff, J. A. Peña, and J. Borggrefe,  Effective spatial resolution of photon counting CT for imaging of trabecular structures is superior to conventional clinical CT and similar to high resolution peripheral CT,  Investigative Radiology 57, 620–626 (2022).


$\displaystyle\leq\frac{1}{2^{j}}\langle\xi\rangle_{k}^{\mu-j-|\alpha|}\Phi(x)^{\mu-j-|\beta|}\Bigg{[}\chi^{(1)}_{N}\langle\xi\rangle_{k}^{\tilde{m}_{1}}\omega(x)^{\tilde{m}_{2}}\tilde{\theta}(t)$
$\displaystyle\qquad+\chi^{(2)}_{N}\langle\xi\rangle_{k}^{m_{1}^{\prime}}\omega(x)^{m_{2}^{\prime}}\left(\frac{\theta(t)}{t}\right)^{m_{3}^{\prime}}\tilde{\theta}(t)^{m_{4}^{\prime}+m_{5}^{\prime}(|\alpha|+|\beta|+2j)}$
$\displaystyle\qquad+\chi^{(3)}_{N}\langle\xi\rangle_{k}^{m_{1}}\omega(x)^{m_{2}}\left(\frac{\theta(t)}{t}\right)^{m_{3}}e^{m_{4}\psi(t)}\tilde{\theta}(t)^{m_{5}+m_{6}(|\alpha|+|\beta|+2j)}\Bigg{]},$

where we have estimated $\frac{\Phi(x)^{\mu}}{2^{j}}\geq 1$and
$\frac{\langle\xi\rangle_{k}^{\mu}}{2^{j}}\geq 1$ (due to the support of cut-
off functions) once in each summand. In $Z_{mid}(N)$ and $Z_{ext}(N),$
$\tilde{\theta}(t)\leq\tilde{\theta}(h).$ For any $r\geq 0$ and $j\geq 1,$

$\tilde{\theta}(t)^{rj}\leq\tilde{\theta}(h)^{rj}\lesssim(\Phi(x)\langle\xi\rangle_{k})^{\varepsilon},$
(A.2)

as the singularity in our consideration is of logarithmic type. From (A.1) and
(A.2), we obtain

$\displaystyle|D_{x}^{\beta}$
$\displaystyle\partial_{\xi}^{\alpha}\gamma_{j}(x,\xi)a_{j}(t,x,\xi)|$ (A.3)
$\displaystyle\leq\frac{1}{2^{j}}\langle\xi\rangle_{k}^{\mu-j-|\alpha|}\Phi(x)^{\mu-j-|\beta|}\Bigg{[}\chi^{(1)}_{N}\langle\xi\rangle_{k}^{\tilde{m}_{1}}\omega(x)^{\tilde{m}_{2}}\tilde{\theta}(t)$
$\displaystyle\qquad+\chi^{(2)}_{N}\langle\xi\rangle_{k}^{m_{1}^{\prime}+\varepsilon}\Phi(x)^{\varepsilon}\omega(x)^{m_{2}^{\prime}}\left(\frac{\theta(t)}{t}\right)^{m_{3}^{\prime}}\tilde{\theta}(t)^{m_{4}^{\prime}+m_{5}^{\prime}(|\alpha|+|\beta|)}$
$\displaystyle\qquad+\chi^{(3)}_{N}\langle\xi\rangle_{k}^{m_{1}+\varepsilon}\Phi(x)^{\varepsilon}\omega(x)^{m_{2}}\left(\frac{\theta(t)}{t}\right)^{m_{3}}e^{m_{4}\psi(t)}\tilde{\theta}(t)^{m_{5}+m_{6}(|\alpha|+|\beta|)}\Bigg{]}$
$\displaystyle\leq\frac{1}{2^{j}}\langle\xi\rangle_{k}^{1-j-|\alpha|}\Phi(x)^{1-j-|\beta|}\Bigg{[}\chi^{(1)}_{N}\langle\xi\rangle_{k}^{\tilde{m}_{1}}\omega(x)^{\tilde{m}_{2}}\tilde{\theta}(t)$
$\displaystyle\qquad+\chi^{(2)}_{N}\langle\xi\rangle_{k}^{m_{1}^{\prime}}\omega(x)^{m_{2}^{\prime}}\left(\frac{\theta(t)}{t}\right)^{m_{3}^{\prime}}\tilde{\theta}(t)^{m_{4}^{\prime}+m_{5}^{\prime}(|\alpha|+|\beta|)}$
$\displaystyle\qquad+\chi^{(3)}_{N}\langle\xi\rangle_{k}^{m_{1}}\omega(x)^{m_{2}}\left(\frac{\theta(t)}{t}\right)^{m_{3}}e^{m_{4}\psi(t)}\tilde{\theta}(t)^{m_{5}+m_{6}(|\alpha|+|\beta|)}\Bigg{]},$

where $j\geq 1.$ Now

$\displaystyle|D_{x}^{\beta}$
$\displaystyle\partial_{\xi}^{\alpha}a(t,x,\xi)|$ $\displaystyle\leq$
$\displaystyle\left|D_{x}^{\beta}\partial_{\xi}^{\alpha}\left(\gamma_{0}(x,\xi)a_{0}(t,x,\xi)\right)\right|+\sum_{j=1}^{j_{0}-1}\left|D_{x}^{\beta}\partial_{\xi}^{\alpha}\left(\gamma_{j}(x,\xi)a_{j}(t,x,\xi)\right)\right|.$

Combining the symbol estimate of $a_{0}$ and the estimate (A.3), we readily
obtain

$\displaystyle a$ $\displaystyle\in
G^{\tilde{m}_{1},\tilde{m}_{2}}\\{\tilde{m}_{3}\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{m_{1}^{\prime},m_{2}^{\prime}}\\{m_{3}^{\prime},m_{4}^{\prime},m_{5}^{\prime}\\}(\omega,g_{\Phi,k})_{N}^{(2)}$
$\displaystyle\qquad\cap
G^{m_{1},m_{2}}\\{m_{3},m_{4},m_{5},m_{6}\\}(\omega,g_{\Phi,k})_{N}^{(3)}$

Arguing as above, we see that $\sum_{j=j_{0}}^{\infty}\gamma_{j}a_{j}$ belongs
to

$\displaystyle
G^{\tilde{m}_{1}-j_{0},1}\\{\tilde{m}_{3}\\}(\omega^{\tilde{m}_{2}}\Phi^{-j_{0}},g_{\Phi,k})^{(1)}_{N}\cap
G^{{m}_{1}^{\prime}-j_{0}+\varepsilon,1}\\{m_{3}^{\prime},m_{4}^{\prime},m_{5}^{\prime}\\}(\omega^{{m}_{2}^{\prime}}\Phi^{-j_{0}+\varepsilon},g_{\Phi,k})^{(2)}_{N}$
$\displaystyle\qquad\cap
G^{{m}_{1}-j_{0}+\varepsilon,1}\\{m_{3},m_{4},m_{5},m_{6}\\}(\omega^{{m}_{2}}\Phi^{-j_{0}+\varepsilon},g_{\Phi,k})^{(2)}_{N},$

and thus, $a(t,x,\xi)-\sum\limits_{j=0}^{j_{0}-1}a_{j}(t,x,\xi)$ belongs to

$\displaystyle
G^{\tilde{m}_{1}-j_{0},1}\\{\tilde{m}_{3}\\}(\omega^{\tilde{m}_{2}}\Phi^{-j_{0}},g_{\Phi,k})^{(1)}_{N}\cap
G^{{m}_{1}^{\prime}-j_{0}+\varepsilon,1}\\{m_{3}^{\prime},m_{4}^{\prime},m_{5}^{\prime}\\}(\omega^{{m}_{2}^{\prime}}\Phi^{-j_{0}+\varepsilon},g_{\Phi,k})^{(2)}_{N}$
$\displaystyle\qquad\cap
G^{{m}_{1}-j_{0}+\varepsilon,1}\\{m_{3},m_{4},m_{5},m_{6}\\}(\omega^{{m}_{2}}\Phi^{-j_{0}+\varepsilon},g_{\Phi,k})^{(2)}_{N}.$

Lastly, we use Propositions A.1 \- A.3 to conclude that

$t^{1-\varepsilon}a_{j}\in
C\left([0,T];G^{m_{1}^{*}-j,1}(\omega^{m_{2}^{*}}\Phi^{r^{*}-j},g_{\Phi,k})\right)$

for
$m_{1}^{*}=\max\\{\tilde{m}_{1},m_{1}^{\prime}+m_{3}^{\prime},m_{1}+m_{3}\\},$
$m_{1}^{*}=\max\\{\tilde{m}_{2},m_{2}^{\prime},m_{2}\\}$ and
$r^{*}=\max\\{m_{3}^{\prime},m_{3}\\}.$ As $j$ tends to $+\infty,$ the
intersection of all those spaces belongs to the space
$C\left((0,T];G^{-\infty}\right).$ This completes the proof. ∎

######  Lemma A.1.

Let $A$ and $B$ be pseudodifferential operators with symbols

$\displaystyle a=\sigma(A)$ $\displaystyle\in
G^{\tilde{m}_{1},\tilde{m}_{2}}\\{\tilde{m}_{3}\\}(\omega,g_{\Phi,k})_{2N}^{(1)}\cap
G^{m_{1}^{\prime},m_{2}^{\prime}}\\{m_{3}^{\prime},m_{4}^{\prime},m_{5}^{\prime}\\}(\omega,g_{\Phi,k})_{N}^{(2)}$
$\displaystyle\qquad\cap
G^{m_{1},m_{2}}\\{m_{3},m_{4},m_{5},m_{6}\\}(\omega,g_{\Phi,k})_{N}^{(3)}$

and

$\displaystyle b=\sigma(B)$ $\displaystyle\in
G^{\tilde{l}_{1},\tilde{l}_{2}}\\{\tilde{l}_{3}\\}(\omega,g_{\Phi,k})_{2N}^{(1)}\cap
G^{l_{1}^{\prime},l_{2}^{\prime}}\\{l_{3}^{\prime},l_{4}^{\prime},l_{5}^{\prime}\\}(\omega,g_{\Phi,k})_{N}^{(2)}$
$\displaystyle\qquad\cap
G^{l_{1},l_{2}}\\{l_{3},l_{4},l_{5},l_{6}\\}(\omega,g_{\Phi,k})_{N}^{(3)}.$

Then, the pseudodifferential operator $C=A\circ B$ has a symbol $c=\sigma(C)$
in

$\displaystyle
G^{\tilde{m}_{1}+\tilde{l}_{1},\tilde{m}_{2}+\tilde{l}_{2}}\\{\tilde{m}_{3}+\tilde{l}_{3}\\}(\omega,g_{\Phi,k})_{2N}^{(1)}\cap
G^{m_{1}^{\prime}+l_{1}^{\prime},m_{2}^{\prime}+l_{2}^{\prime}}\\{m_{3}^{\prime}+l_{3}^{\prime},m_{4}+l_{4}^{\prime},m_{5}^{\prime}+l^{\prime}_{5}\\}(\omega,g_{\Phi,k})_{N}^{(2)}$
$\displaystyle\qquad\cap
G^{m_{1}+l_{1},m_{2}+l_{2}}\\{m_{3}+l_{3},m_{4}+l_{4},m_{5}+l_{5},m_{6}+l_{6}\\}(\omega,g_{\Phi,k})_{N}^{(3)}$

and satisfies

$c(t,x,\xi)\sim\sum_{\alpha\in\mathbb{N}^{n}}\frac{1}{\alpha!}\partial_{\xi}^{\alpha}a(t,x,\xi)D_{x}^{\alpha}b(t,x,\xi).$
(A.4)

The operator $C$ is uniquely determined modulo an operator with symbol from
$C\left((0,T];G^{-\infty}\right).$

###### Proof.

In view of Propositions A.1 \- A.4, it is clear that the operator $C$ is a
well-defined pseudodifferential operator. Relation (A.4) is a direct
consequence of the standard composition rules (see [17, Section 1.2]). ∎

######  Lemma A.2.

Let $A$ be a pseudodifferential operator with an invertible symbol

$a=\sigma(A)\in G^{0,0}\\{0\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{0,0}\\{0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(2)}\cap
G^{0,0}\\{0,0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(3)}.$

Then, there exists a parametrix $A^{\\#}$ with symbol $a^{\\#}$ in

$G^{0,0}\\{0\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{0,0}\\{0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(2)}\cap
G^{0,0}\\{0,0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(3)}.$

###### Proof.

We use the existence of the inverse of $a$ and set

$\displaystyle a_{0}^{\\#}(t,x,\xi)=a(t,x,\xi)^{-1}$ $\displaystyle\in
G^{0,0}\\{0\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{0,0}\\{0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(2)}$ $\displaystyle\qquad\cap
G^{0,0}\\{0,0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(3)}.$

In view of Propositions A.1 \- A.3, one can define a sequence
$a_{j}^{\\#}(t,x,\xi)$ recursively by

$\sum_{1\leq|\alpha|\leq
j}\frac{1}{\alpha!}\partial_{\xi}^{\alpha}a(t,x,\xi)D_{x}^{\alpha}a_{j-|\alpha|}^{\\#}(t,x,\xi)=-a(t,x,\xi)a_{j}^{\\#}(t,x,\xi)$

with $a_{j}^{\\#}$ in

$G^{-j,1}\\{0\\}(\Phi^{-j},g_{\Phi,k})^{(1)}_{N}\cap
G^{-j,1}\\{0,0,0\\}(\Phi^{-j},g_{\Phi,k})^{(2)}_{N}\cap
G^{-j,1}\\{0,0,0,0\\}(\Phi^{-j},g_{\Phi,k})^{(3)}_{N}.$

Proposition A.4 then yields the existence of a symbol

$a_{R}^{\\#}\in G^{0,0}\\{0\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{0,0}\\{0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(2)}\cap
G^{0,0}\\{0,0,0,0\\}(\omega,g_{\Phi,k})_{N}^{(3)}$

and a right parametrix $A_{R}^{\\#}(t,x,\xi)$ with symbol
$\sigma\left(A_{R}^{\\#}\right)=a_{R}^{\\#}.$ We have

$AA_{R}^{\\#}-I\in C\left([0,T];G^{-\infty}\right).$

The existence of a left parametrix follows in similar lines. One can also
prove the existence of a parametrix $A^{\\#}$ by showing that right and left
parametrix coincide up to a regularizing operator. ∎

##  Appendix B Appendix II: Sharp Gårding Inequality For Parameter Dependent
Matrix

In this section we prove the sharp Gårding inequality for a matrix
pseudodifferential operator with symbol $a(t,x,\xi)$ in

$G^{1,1}\\{1\\}(\omega,g_{\Phi,k})_{N}^{(1)}\cap
G^{0,0}\\{1,0,1\\}(\omega,g_{\Phi,k})_{N}^{(2)}\cap
G^{-1,-1}\\{2,1,2,1\\}(\omega,g_{\Phi,k})_{N}^{(3)}.$ (B.1)

This is used Sections 5.4 and 5.5 to arrive at an energy estimate. By
Propositions A.1 and A.2, for any $\varepsilon\in(0,1),$

$\displaystyle a$ $\displaystyle\in
L^{\infty}\left([0,1];G^{1,1}(\Phi,g_{\Phi,k})\right),$
$\displaystyle\qquad\text{ if }\tilde{\theta}\text{ is bounded,}$ (B.2)
$\displaystyle t^{1-\varepsilon}a$ $\displaystyle\in
C\left([0,1];G^{1,1}(\Phi,g_{\Phi,k})\right),$ otherwise. (B.3)

We give below the sharp Gårding inequality when $\tilde{\theta}$ is unbounded.
The case for bounded $\tilde{\theta}$ follows in similar lines by taking
$\varepsilon=1$ in the proof.

######  Theorem B.1.

Let $a(t,x,\xi)$ be $2\times 2$ positive semi-definite matrix belonging to the
symbol class given in (B.1). Then, for each $t\in[0,T]$ and
$\varepsilon\in(0,1)$ we have

$\left(t^{1-\varepsilon}a^{w}(t,x,D)u,u\right)\geqq-C\|u\|^{2},\quad
u\in\mathcal{S}\left(\mathbb{R}^{n};{\mathbb{R}}^{2}\right).$ (B.4)

###### Proof.

We have the following region-wise seminorms for $a(t,x,\xi)$:

$\displaystyle\text{in }Z_{int}(N),\qquad$ $\displaystyle|a|_{j}^{g}(w)$
$\displaystyle\leq h^{-1}\tilde{\theta},$ $\displaystyle\text{in }Z_{mid}(N),$
$\displaystyle|a|_{j}^{g}(w)$
$\displaystyle\leq\frac{\theta(t)}{t}\tilde{\theta}^{j},$
$\displaystyle\text{in }Z_{ext}(N),$ $\displaystyle|a|_{j}^{g}(w)$
$\displaystyle\leq
h\left(\frac{\theta(t)}{t}\right)^{2}e^{\psi(t)}{\tilde{\theta}}^{2+j}.$

Here $w\in{\mathbb{R}}^{2n}$, $g=g_{\Phi,k}$ and the norm of the $j^{th}$
differential at $x$ given by

$|f|_{j}^{g}:=\sup_{y_{l}\in{\mathbb{R}}^{2n}}|f^{(j)}(x;y_{1},\dots,y_{j})|\Big{/}\prod_{1}^{j}g(y_{l})^{\frac{1}{2}}.$

Using an appropriate affine symplectic transformation one can assume that
$g_{\Phi,k}=he$ where $e$ is the Euclidean metric form. Then

$\displaystyle\text{in }Z_{int}(N),\qquad$ $\displaystyle|a|_{j}^{e}(w)$
$\displaystyle\leq h^{(j-2)/2}\tilde{\theta},$ (B.5) $\displaystyle\text{in
}Z_{mid}(N),$ $\displaystyle|a|_{j}^{e}(w)$ $\displaystyle\leq
h^{j/2}\frac{\omega(t)}{t}\tilde{\theta}^{j},$ (B.6) $\displaystyle\text{in
}Z_{ext}(N),$ $\displaystyle|a|_{k}^{e}(w)$ $\displaystyle\leq
h^{(j+2)/2}\left(\frac{\omega(t)}{t}\right)^{2}e^{\psi(t)}{\tilde{\theta}}^{2+j}.$
(B.7)

Let $a^{(0)}_{t}+a^{(1)}_{t}(x,\xi)$ be the first order Taylor expansion of
$a$ at $(x,\xi)=0.$ Let $v\in{\mathbb{R}}^{2}$. By the above semi-norms we
have

$\displaystyle\text{in }Z_{int}(N),\quad$
$\displaystyle\begin{cases}\begin{array}[]{l}\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(x,0)v,v\right)+\frac{|x|^{2}}{2}(\tilde{\theta}v,v)\geqq
0,\\\
\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(0,\xi)v,v\right)+\frac{|\xi|^{2}}{2}(\tilde{\theta}v,v)\geqq
0,\end{array}\end{cases}$ (B.8) $\displaystyle\text{in }Z_{mid}(N),\quad$
$\displaystyle\begin{cases}\begin{array}[]{l}\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(x,0)v,v\right)+\frac{|x|^{2}}{2}\left(\frac{\theta(t)}{t}\tilde{\theta}^{2}hv,v\right)\geqq
0,\\\
\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(0,\xi)v,v\right)+\frac{|\xi|^{2}}{2}\left(\frac{\theta(t)}{t}\tilde{\theta}^{2}hv,v\right)\geqq
0,\end{array}\end{cases}$ (B.9) $\displaystyle\text{in }Z_{ext}(N),\quad$
$\displaystyle\begin{cases}\begin{array}[]{l}\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(x,0)v,v\right)+\frac{|x|^{2}}{2}\left(\frac{\theta(t)^{2}}{t^{2}}e^{\psi(t)}\tilde{\theta}^{4}h^{2}v,v\right)\geqq
0,\\\
\left(a^{(0)}_{t}v,v\right)+\left(a^{(1)}_{t}(0,\xi)v,v\right)+\frac{|\xi|^{2}}{2}\left(\frac{\theta(t)^{2}}{t^{2}}e^{\psi(t)}\tilde{\theta}^{4}h^{2}v,v\right)\geqq
0.\end{array}\end{cases}$ (B.10)

Let us fix $0<\varepsilon<\varepsilon^{\prime}<1.$ Using the definition of the
regions, we see that in the whole of extended phase space

$\begin{cases}\begin{array}[]{l}\left(t^{1-\varepsilon^{\prime}}a_{t}^{(0)}v,v\right)+\left(t^{1-\varepsilon^{\prime}}a_{t}^{(1)}(x,0)v,v\right)+\frac{|x|^{2}}{2}(t^{1-\varepsilon^{\prime}}\tilde{\theta}(t)^{q}v,v)\geqq
0,\\\
\left(t^{1-\varepsilon^{\prime}}a_{t}^{(0)}v,v\right)+\left(t^{1-\varepsilon^{\prime}}a_{t}^{(1)}(0,\xi)v,v\right)+\frac{|\xi|^{2}}{2}(t^{1-\varepsilon^{\prime}}\tilde{\theta}(t)^{q}v,v)\geqq
0,\end{array}\end{cases}$

where $q=1$ in $Z_{int}(N)$ while in $Z_{mid}(N)$ and $Z_{ext}(N)$ $q=2.$
Since the function $\tilde{\theta}(t)$ is of logarithmic type, in the whole of
extended phase space, we have

$\begin{cases}\begin{array}[]{l}t^{1-\varepsilon}\left((a_{t}^{(0)}v,v)+(a_{t}^{(1)}(x,0)v,v)+\frac{|x|^{2}}{2}\|v\|^{2}\right)\geqq
0,\\\
t^{1-\varepsilon}\left((a_{t}^{(0)}v,v)+(a_{t}^{(1)}(0,\xi)v,v)+\frac{|\xi|^{2}}{2}\|v\|^{2}\right)\geqq
0.\end{array}\end{cases}$

From here on we proceed as in [12, Theorem 18.6.14] to obtain the result. ∎

## References

  * [1] Ascanelli, A. and Cappiello, M. (2006) Log-Lipschitz regularity for SG hyperbolic systems. J. Differ. Equ., 230, 556 – 578.
  * [2] Ascanelli, A. and Cappiello, M. (2008) Hölder continuity in time for SG hyperbolic systems. J. Differ. Equ., 244, 2091 – 2121.
  * [3] Beals, R. and Fefferman, C. (1974) Spatially inhomogeneous pseudodifferential operators I. Commun. Pur. Appl. Math., 27, 1–24.
  * [4] Bernardi, E., Bove, A., and Petkov, V. (2015) Cauchy problem for effectively hyperbolic operators with triple characteristics of variable multiplicity. J. Hyperbolic Differ. Equ., 12, 535–579.
  * [5] Cicognani, M. (2003) The Cauchy problem for strictly hyperbolic operators with non-absolutely continuous coefficients. Tsukuba J. Math., 27, 1–12.
  * [6] Cicognani, M. (2004) Coefficients with unbounded derivatives in hyperbolic equations. Math. Nachr., 276, 31–46.
  * [7] Cicognani, M. and Lorenz, D. (2018) Strictly hyperbolic equations with coefficients low-regular in time and smooth in space. J. Pseudo Differ. Oper. Appl., 9, 643–675.
  * [8] Colombini, F., Del Santo, D., and Kinoshita, T. (2002) Well-posedness of the Cauchy problem for a hyperbolic equation with non-Lipschitz coefficients. Ann. Scoula. Norm.-Sci., Ser. 5, 1, 327–358.
  * [9] Colombini, F., Del Santo, D., and Reissig, M. (2003) On the optimal regularity of coefficients in hyperbolic cauchy problems. Bulletin des Sciences Mathématiques, 127, 328–347.
  * [10] Ghisi, M. and Gobbino, M. (2021) Finite vs infinite derivative loss for abstract wave equations with singular time-dependent propagation speed. Bulletin des Sciences Mathématiques, 166, 102918.
  * [11] Hirosawa, F. (2003) Loss of regularity for the solutions to hyperbolic equations with non-regular coefficients—an application to kirchhoff equation. Mathematical Methods in the Applied Sciences, 26, 783–799.
  * [12] Hörmander, L. (1985) The Analysis of Linear Partial Differential Operators III. Springer-Verlag.
  * [13] Kajitani, K. and Nishitani, T. (1991) The Hyperbolic Cauchy Problem. Springer, Berlin, Heidelberg.
  * [14] Kubo, A. and Reissig, M. (2003) Construction of parametrix to strictly hyperbolic Cauchy problems with fast oscillations in nonLipschitz coefficients. Commun. Part. Diff. Eq., 28, 1471–1502.
  * [15] Lerner, N. (2010) Metrics on the Phase Space and Non-Selfadjoint Pseudo-Differential Operators. Birkhäuser Basel.
  * [16] Lorenz, D. and Reissig, M. (2019) A generalized levi condition for weakly hyperbolic cauchy problems with coefficients low regular in time and smooth in space. J. Pseudo Differ. Oper. Appl., 10, 199–240.
  * [17] Nicola, F. and Rodino, L. (2011) Global Pseudo-differential Calculus on Euclidean Spaces. Birkhäuser Basel.
  * [18] Pattar, R. R. and Uday Kiran, N. (2021), Global well-posedness of a class of strictly hyperbolic cauchy problems with coefficients non-absolutely continuous in time. Bulletin des Sciences Mathématiques , 171, 103037.
  * [19] Pattar, R. R., Kiran, N. U. (2021) Strictly hyperbolic Cauchy problems on ${\mathbb{R}}^{n}$ with unbounded and singular coefficients. Annali dell’Universita di Ferrara, https://doi.org/10.1007/s11565-021-00378-2.
  * [20] Rauch, J. (2012) Hyperbolic Partial Differential Equations and Geometric Optics. Graduate studies in mathematics, American Mathematical Society.
  * [21] Reissig, M. (2003) Hypebolic equations with non-Lipschitz coefficients. Rend. Sem. Mat. Univ. Pol. Torino, 61, 135–182.
  * [22] Uday Kiran, N., Coriasco, S., and Battisti, U. (2016) Hyperbolic operators with non-Lipschitz coefficients. Recent Advances in Theoretical & Computational Partial Differential Equations with Applications, Panjab University, Chandigarh.
  * [23] Yagdjian, K. (1997) The Cauchy Problem for Hyperbolic Operators: Multiple Characteristics. Micro-Local Approach. Wiley.